G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsApproxC1Spline.hpp
Go to the documentation of this file.
1
14#pragma once
15
16#include<gsUnstructuredSplines/src/gsApproxC1Edge.h>
17#include<gsUnstructuredSplines/src/gsApproxC1Vertex.h>
18
19#include <gsUnstructuredSplines/src/gsApproxC1Utils.h>
20
21namespace gismo
22{
23
24template<short_t d,class T>
25void gsApproxC1Spline<d,T>::defaultOptions()
26{
27 /*
28 to do: general
29 */
30 gsTensorBSplineBasis<d, T> basis = dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(0));
31 index_t p = basis.degree(0);
32
33 // Check if the degree is the same for each patch
34 for (size_t np = 0; np < m_patches.nPatches(); np++) {
35 gsTensorBSplineBasis<d, T> basis_np = dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(np));
36 if (p != basis_np.degree(0))
37 gsWarn << "Not suitable for different degrees! \n";
38 }
39
40 // Check if the multipatch has no inner knots
41 for (size_t np = 0; np < m_patches.nPatches(); np++)
42 for (index_t dir = 0; dir < 2; dir++) {
43 const gsBSplineBasis<real_t> basis_geo =
44 dynamic_cast<const gsBSplineBasis<real_t> &>(m_patches[np].basis().component(dir));
45 gsKnotVector<real_t> kv_geo = dynamic_cast<const gsKnotVector<real_t> &>(basis_geo.knots());
46 std::vector<real_t> unique_geo = kv_geo.unique();
47
48 if (2 != unique_geo.size())
49 gsWarn << "Not suitable for geometries with inner knots! \n";
50 }
51
52 m_options.addSwitch("info","Print debug information", false );
53 m_options.addSwitch("plot","Print debug information", false );
54 m_options.addSwitch("interpolation","Compute the basis with interpolation", false );
55 m_options.addSwitch("second","Compute the second biharmonic problem", false );
56
57 m_options.addInt("gluingDataDegree","Print debug information", -1 );
58 m_options.addInt("gluingDataSmoothness","Print debug information", -1 );
59}
60
61template<short_t d,class T>
62void gsApproxC1Spline<d,T>::init()
63{
64 index_t row_dofs = 0;
65
66 // Fix orientation of each patch
67 m_patches.fixOrientation();
68
69 m_bases.clear();
70 m_bases.reserve(m_patches.nPatches()); // For each Patch
71 for (size_t np = 0; np < m_patches.nPatches(); np++)
72 {
73 // gsContainerBasisBase:
74 // # basisContainer
75 // - Interior space: [0] : inner,
76 // - Edge spaces: [1] : west, [2] : east, [3] : south, [4] : north,
77 // - Vertex spaces: [5] : southwest, [6] : southeast, [7] : northwest, [8] : northeast
78 // [9] For the initial space, maybe delete later
79 gsContainerBasis<d,T> containerBasis(9); // for 9 subspaces
80 m_bases.push_back(containerBasis);
81 }
82
83 // Create interior spline space
84 for (size_t np = 0; np < m_patches.nPatches(); np++)
85 {
86 gsTensorBSplineBasis<d, T> basis_inner = dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(np));
87
88 index_t dim_u = basis_inner.component(0).size();
89 index_t dim_v = basis_inner.component(1).size();
90
91 row_dofs += (dim_u - 4) * (dim_v - 4);
92
93 m_bases[np].setBasis(0, basis_inner); // Inner
94 }
95
96 // For loop over Interface to construct the spaces
97 for (size_t numInt = 0; numInt < m_patches.interfaces().size(); numInt++)
98 {
99 const boundaryInterface & item = m_patches.interfaces()[numInt];
100
101 const index_t dir_1 = item.first().side().index() > 2 ? 0 : 1;
102 const index_t dir_2 = item.second().side().index() > 2 ? 0 : 1;
103
104 // [!Basis space]
105 gsTensorBSplineBasis<d, T> basis_11 = dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(item.first().patch));
106 gsTensorBSplineBasis<d, T> basis_22 = dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(item.second().patch));
107
108 gsTensorBSplineBasis<d, T> basis;
109 index_t dir, patch;
110 if (basis_11.component(dir_1).numElements() > basis_22.component(dir_2).numElements())
111 {
112 basis = basis_22;
113 dir = dir_2;
114 patch = item.second().patch;
115 basis_11.component(dir_1) = basis_22.component(dir_2);
116 }
117 else
118 {
119 basis = basis_11;
120 dir = dir_1;
121 patch = item.first().patch;
122 basis_22.component(dir_2) = basis_11.component(dir_1);
123 }
124 // [!Basis space]
125
126 // [!Plus Minus space]
127 gsBSplineBasis<T> basis_plus;
128 gsBSplineBasis<T> basis_minus;
129 createPlusSpace(m_patches.patch(patch), basis, dir, basis_plus);
130 createMinusSpace(m_patches.patch(patch), basis, dir, basis_minus);
131 // [!Plus Minus space]
132
133 // [!Gluing data space]
134 gsBSplineBasis<T> basis_gluingData;
135 createGluingDataSpace(m_patches.patch(patch), basis, dir, basis_gluingData,
136 m_options.getInt("gluingDataDegree"), m_options.getInt("gluingDataSmoothness"));
137 // [!Gluing data space]
138
139 // [!Edge space]
140 gsTensorBSplineBasis<d, T> basis_edge_11, basis_edge_22;
141 createEdgeSpace(m_patches.patch(item.first().patch), basis_11, dir_1,
142 basis_plus, basis_minus, basis_gluingData, basis_edge_11);
143 createEdgeSpace(m_patches.patch(item.second().patch), basis_22, dir_2,
144 basis_plus, basis_minus, basis_gluingData, basis_edge_22);
145
146 // [!Edge space]
147 m_bases[item.first().patch].setBasis(item.first().side().index(), basis_edge_11);
148 m_bases[item.second().patch].setBasis(item.second().side().index(), basis_edge_22);
149
150 index_t numDofs = std::max(basis_plus.size() + basis_minus.size() - 10, (index_t)0);
151 row_dofs += numDofs; // The same as side_2
152 }
153
154 // For loop over the Edge to construct the spaces
155 for (size_t numBdy = 0; numBdy < m_patches.boundaries().size(); numBdy++)
156 {
157 const patchSide & bit = m_patches.boundaries()[numBdy];
158
159 index_t dir_1 = m_patches.boundaries()[numBdy].m_index < 3 ? 1 : 0;
160
161 gsTensorBSplineBasis<d, T> basis = dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(bit.patch));
162
163 // [!Plus Minus space]
164 gsBSplineBasis<T> basis_plus;
165 gsBSplineBasis<T> basis_minus;
166 createPlusSpace(m_patches.patch(bit.patch), basis, dir_1, basis_plus);
167 createMinusSpace(m_patches.patch(bit.patch), basis, dir_1, basis_minus);
168 // [!Plus Minus space]
169
170 gsTensorBSplineBasis<d, T> basis_edge_11;
171 createEdgeSpace(m_patches.patch(bit.patch), basis, dir_1, basis_plus, basis_minus, basis_edge_11);
172
173 m_bases[bit.patch].setBasis(bit.side().index(), basis_edge_11);
174
175 index_t numDofs = std::max(basis_plus.size() + basis_minus.size() - 10, (index_t)0);
176 row_dofs += numDofs;
177 }
178
179 // For loop over the Vertex to construct the spaces
180 for (size_t numVer = 0; numVer < m_patches.vertices().size(); numVer++)
181 {
182 std::vector<patchCorner> allcornerLists = m_patches.vertices()[numVer];
183 for (size_t j = 0; j < allcornerLists.size(); j++)
184 {
185 std::vector<patchSide> containingSides;
186 allcornerLists[j].getContainingSides(d, containingSides);
187 bool isInterface_1 = m_patches.isInterface(patchSide(allcornerLists[j].patch,
188 containingSides.at(0).side()));
189 bool isInterface_2 = m_patches.isInterface(patchSide(allcornerLists[j].patch,
190 containingSides.at(1).side()));
191
192 gsTensorBSplineBasis<d, T> basis =
193 dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(allcornerLists[j].patch));
194 if (isInterface_1)
195 {
196 patchSide result;
197 m_patches.getNeighbour(containingSides.at(0), result);
198
199 gsTensorBSplineBasis<d, T> basis2 =
200 dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(result.patch));
201 index_t dir_1 = containingSides.at(0).side() < 3 ? 1 : 0;
202 index_t dir_2 = result.side().index() < 3 ? 1 : 0;
203 if (basis.component(dir_1).numElements() > basis2.component(dir_2).numElements())
204 basis.component(dir_1) = basis2.component(dir_2);
205 }
206 if (isInterface_2)
207 {
208 patchSide result;
209 m_patches.getNeighbour(containingSides.at(1), result);
210
211 gsTensorBSplineBasis<d, T> basis2 =
212 dynamic_cast<gsTensorBSplineBasis<d, T> &>(m_multiBasis.basis(result.patch));
213 index_t dir_1 = containingSides.at(1).side() < 3 ? 1 : 0;
214 index_t dir_2 = result.side().index() < 3 ? 1 : 0;
215 if (basis.component(dir_1).numElements() > basis2.component(dir_2).numElements())
216 basis.component(dir_1) = basis2.component(dir_2);
217 }
218
219 if (containingSides.at(0).side() < 3) // If isInterface_1 == v, then switch
220 {
221 bool isInterface_temp = isInterface_1;
222 isInterface_1 = isInterface_2;
223 isInterface_2 = isInterface_temp;
224 }
225
226 gsTensorBSplineBasis<d, T> basis_vertex_11;
227 createVertexSpace(m_patches.patch(allcornerLists[j].patch), basis,
228 isInterface_1, isInterface_2, basis_vertex_11, m_options.getInt("gluingDataDegree"),
229 m_options.getInt("gluingDataSmoothness"));
230
231 m_bases[allcornerLists[j].patch].setBasis(allcornerLists[j].m_index + 4, basis_vertex_11);
232 }
233 row_dofs += 6;
234 }
235
236 // Initialise the matrix
237 m_matrix.clear();
238 index_t dim_col = 0;
239 for (size_t i = 0; i < m_bases.size(); i++)
240 {
241 dim_col += m_bases[i].size();
242 }
243
244 m_matrix.resize(row_dofs, dim_col);
245 const index_t nz = (m_multiBasis.basis(0).degree(0)*2 + 1)*row_dofs; // TODO
246 m_matrix.reserve(nz);
247
248
249 if (m_options.getSwitch("info"))
250 {
251 for (size_t i = 0; i < m_bases.size(); i++)
252 {
253 gsInfo << "Patch " << i << ": \n";
254 for (index_t j = 0; j < m_bases[i].nPieces(); j++)
255 {
256 gsInfo << (j == 0 ? "Interior Space: " : ( j < 5 ? "Edge Space: " : "Vertex Space: ") ) << "\n";
257 gsTensorBSplineBasis<d, T> basis =
258 dynamic_cast<const gsTensorBSplineBasis<d, T>&>(m_bases[i].piece(j));
259 std::vector<index_t> vec_1 = basis.knots(0).multiplicities();
260 std::vector<index_t> vec_2 = basis.knots(1).multiplicities();
261 gsAsVector<index_t> mult_1(vec_1);
262 gsAsVector<index_t> mult_2(vec_2);
263 gsInfo << mult_1.transpose() << "\n";
264 gsInfo << mult_2.transpose() << "\n";
265 gsInfo << "----------------------------------\n";
266 }
267 }
268 }
269}
270
271
272template<short_t d,class T>
273void gsApproxC1Spline<d,T>::compute()
274{
275 // Compute Inner Basis functions
276 index_t shift_row = 0, shift_col = 0;
277 for(size_t np = 0; np < m_patches.nPatches(); ++np)
278 {
279 index_t dim_u = m_bases[np].piece(0).component(0).size();
280 index_t dim_v = m_bases[np].piece(0).component(1).size();
281
282 index_t row_i = 0;
283 for (index_t j = 2; j < dim_v-2; ++j)
284 for (index_t i = 2; i < dim_u-2; ++i)
285 {
286 m_matrix.insert(shift_row + row_i, shift_col + j*dim_u+i) = 1.0;
287 ++row_i;
288 }
289
290 shift_row += row_i;
291 shift_col += m_bases[np].size();
292 }
293
294 // Interfaces
295 for (size_t numInt = 0; numInt < m_patches.interfaces().size(); numInt++)
296 {
297 const boundaryInterface & item = m_patches.interfaces()[numInt];
298 index_t side_1 = item.first().side().index();
299 index_t side_2 = item.second().side().index();
300
301 index_t patch_1 = item.first().patch;
302 index_t patch_2 = item.second().patch;
303
304 gsApproxC1Edge<d, T> approxC1Edge(m_patches, m_bases, item, numInt, m_options);
305 std::vector<gsMultiPatch<T>> basisEdge = approxC1Edge.getEdgeBasis();
306
307 index_t begin_col = 0, end_col = 0, shift_col = 0;
308 for (index_t np = 0; np < patch_1; ++np)
309 shift_col += m_bases[np].size();
310
311 for (index_t ns = 0; ns < side_1; ++ns)
312 begin_col += m_bases[patch_1].piece(ns).size();
313 for (index_t ns = 0; ns < side_1+1; ++ns)
314 end_col += m_bases[patch_1].piece(ns).size();
315
316 for (size_t ii = 0; ii < basisEdge[0].nPatches(); ++ii)
317 {
318 index_t jj = 0;
319 for (index_t j = begin_col; j < end_col; ++j, ++jj) {
320 if (basisEdge[0].patch(ii).coef(jj, 0) * basisEdge[0].patch(ii).coef(jj, 0) > 1e-25)
321 m_matrix.insert(shift_row + ii, shift_col + j) = basisEdge[0].patch(ii).coef(jj, 0);
322 }
323 }
324
325 begin_col = 0, end_col = 0, shift_col = 0;
326 for (index_t np = 0; np < patch_2; ++np)
327 shift_col += m_bases[np].size();
328
329 for (index_t ns = 0; ns < side_2; ++ns)
330 begin_col += m_bases[patch_2].piece(ns).size();
331 for (index_t ns = 0; ns < side_2+1; ++ns)
332 end_col += m_bases[patch_2].piece(ns).size();
333
334 for (size_t ii = 0; ii < basisEdge[1].nPatches(); ++ii)
335 {
336 index_t jj = 0;
337 for (index_t j = begin_col; j < end_col; ++j, ++jj) {
338 if (basisEdge[1].patch(ii).coef(jj, 0) * basisEdge[1].patch(ii).coef(jj, 0) > 1e-25)
339 m_matrix.insert(shift_row + ii, shift_col + j) = basisEdge[1].patch(ii).coef(jj, 0);
340 }
341 }
342
343 shift_row += basisEdge[0].nPatches();
344 }
345
346 // Compute Edge Basis functions
347 for (size_t numBdy = 0; numBdy < m_patches.boundaries().size(); numBdy++)
348 {
349 const patchSide & bit = m_patches.boundaries()[numBdy];
350 index_t side_1 = bit.side().index();
351 index_t patch_1 = bit.patch;
352
353 gsApproxC1Edge<d, T> approxC1Edge(m_patches, m_bases, bit, numBdy, m_options);
354 std::vector<gsMultiPatch<T>> basisEdge = approxC1Edge.getEdgeBasis();
355
356 index_t begin_col = 0, end_col = 0, shift_col = 0;
357 for (index_t np = 0; np < patch_1; ++np)
358 shift_col += m_bases[np].size();
359
360 for (index_t ns = 0; ns < side_1; ++ns)
361 begin_col += m_bases[patch_1].piece(ns).size();
362 for (index_t ns = 0; ns < side_1+1; ++ns)
363 end_col += m_bases[patch_1].piece(ns).size();
364
365 for (size_t ii = 0; ii < basisEdge[0].nPatches(); ++ii)
366 {
367 index_t jj = 0;
368 for (index_t j = begin_col; j < end_col; ++j, ++jj) {
369 if (basisEdge[0].patch(ii).coef(jj, 0) * basisEdge[0].patch(ii).coef(jj, 0) > 1e-20)
370 m_matrix.insert(shift_row + ii, shift_col + j) = basisEdge[0].patch(ii).coef(jj, 0);
371 }
372 }
373 shift_row += basisEdge[0].nPatches();
374 }
375
376 // Compute Vertex Basis functions
377 for (size_t numVer = 0; numVer < m_patches.vertices().size(); numVer++)
378 {
379 std::vector<patchCorner> allcornerLists = m_patches.vertices()[numVer];
380 std::vector<size_t> patchIndex;
381 std::vector<size_t> vertIndex;
382 for (size_t j = 0; j < allcornerLists.size(); j++)
383 {
384 patchIndex.push_back(allcornerLists[j].patch);
385 vertIndex.push_back(allcornerLists[j].m_index);
386 }
387
388 gsApproxC1Vertex<d, T> approxC1Vertex(m_patches, m_bases, patchIndex, vertIndex, numVer, m_options);
389 std::vector<gsMultiPatch<T>> basisVertex = approxC1Vertex.getVertexBasis();
390
391 for (size_t np = 0; np < patchIndex.size(); ++np)
392 {
393 index_t patch_1 = patchIndex[np];
394 index_t corner = vertIndex[np];
395
396 index_t begin_col = 0, end_col = 0, shift_col = 0;
397 for (index_t np = 0; np < patch_1; ++np)
398 shift_col += m_bases[np].size();
399
400 for (index_t ns = 0; ns < corner+4; ++ns)
401 begin_col += m_bases[patch_1].piece(ns).size();
402 for (index_t ns = 0; ns < corner+4+1; ++ns)
403 end_col += m_bases[patch_1].piece(ns).size();
404
405 for (size_t ii = 0; ii < basisVertex[np].nPatches(); ++ii)
406 {
407 index_t jj = 0;
408 for (index_t j = begin_col; j < end_col; ++j, ++jj) {
409 if (basisVertex[np].patch(ii).coef(jj, 0) * basisVertex[np].patch(ii).coef(jj, 0) > 1e-20)
410 {
411 m_matrix.insert(shift_row + ii, shift_col + j) = basisVertex[np].patch(ii).coef(jj, 0);
412 }
413 }
414 }
415 }
416 shift_row += basisVertex[0].nPatches(); // + 6
417 }
418 m_matrix.makeCompressed();
419}
420
421
422} // namespace gismo
gismo::boxSide::index
short_t index() const
Returns the index (as specified in boundary::side) of the box side.
Definition gsBoundary.h:140
gismo::gsApproxC1Spline::init
void init()
Initializes the method.
Definition gsApproxC1Spline.hpp:62
gismo::gsApproxC1Spline::defaultOptions
void defaultOptions()
Sets the default options.
Definition gsApproxC1Spline.hpp:25
gismo::gsApproxC1Spline::compute
void compute()
Computes the basis.
Definition gsApproxC1Spline.hpp:273
gismo::gsAsVector
Creates a mapped object or data pointer to a vector without copying data.
Definition gsAsMatrix.h:239
gismo::gsBSplineBasis
A univariate B-spline basis.
Definition gsBSplineBasis.h:700
gismo::gsKnotVector
Class for representing a knot vector.
Definition gsKnotVector.h:80
gismo::gsKnotVector::size
size_t size() const
Number of knots (including repetitions).
Definition gsKnotVector.h:242
gismo::gsTensorBSplineBasis< 1, T >::knots
const KnotVectorType & knots(int i=0) const
Returns the knot vector of the basis.
Definition gsBSplineBasis.h:371
gismo::gsTensorBSplineBasis< 1, T >::size
index_t size() const
Returns the number of elements in the basis.
Definition gsBSplineBasis.h:165
gismo::gsTensorBSplineBasis
A tensor product B-spline basis.
Definition gsTensorBSplineBasis.h:37
gismo::gsTensorBSplineBasis::component
const Basis_t & component(short_t dir) const
For a tensor product basis, return the (const) 1-d basis for the i-th parameter component.
Definition gsTensorBSplineBasis.h:202
gismo::gsTensorBasis::degree
short_t degree(short_t i) const
Returns the degree of the basis wrt variable i.
Definition gsTensorBasis.h:465
gismo::gsTensorBasis::size
index_t size() const
Returns the number of elements in the basis.
Definition gsTensorBasis.h:108
gsApproxC1Edge.h
Creates the (approx) C1 Edge space.
gsApproxC1Utils.h
Provides declaration of Basis abstract interface.
gsApproxC1Vertex.h
Creates the (approx.) C1 Vertex space.
index_t
#define index_t
Definition gsConfig.h:32
gsWarn
#define gsWarn
Definition gsDebug.h:50
gsInfo
#define gsInfo
Definition gsDebug.h:43
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::boundaryInterface
Struct which represents an interface between two patches.
Definition gsBoundary.h:650
gismo::boundaryInterface::first
patchSide & first()
first, returns the first patchSide of this interface
Definition gsBoundary.h:776
gismo::boundaryInterface::second
patchSide & second()
second, returns the second patchSide of this interface
Definition gsBoundary.h:782
gismo::patchSide
Struct which represents a certain side of a patch.
Definition gsBoundary.h:232
gismo::patchSide::patch
index_t patch
The index of the patch.
Definition gsBoundary.h:234