G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsDPatchBase.hpp
Go to the documentation of this file.
1
14#include <gsNurbs/gsTensorNurbsBasis.h>
15#include <gsHSplines/gsHTensorBasis.h>
16#include <gsHSplines/gsTHBSpline.h>
17
18#include <gsAssembler/gsExprHelper.h>
19#include <gsAssembler/gsExprEvaluator.h>
20#include <gsAssembler/gsAssembler.h>
21
22namespace gismo
23{
24
25 template<short_t d,class T>
26 void gsDPatchBase<d,T>::compute()
27 {
28 _initialize();
29 _computeMapper();
30 _computeSmoothMatrix();
31 GISMO_ASSERT(this->_checkMatrix(m_matrix),"Mapper does not have column sum equal to 1");
32 _makeTHB();
33 _computeEVs();
34 GISMO_ASSERT(this->_checkMatrix(m_matrix),"Mapper does not have column sum equal to 1");
35 m_computed=true;
36 }
37
38 template<short_t d,class T>
39 void gsDPatchBase<d,T>::defaultOptions()
40 {
41 m_options.addSwitch("SharpCorners","Reproduce sharp corners",true);
42 m_options.addReal("SharpCornerTolerance","Sharp corner tolerance",1e-2);
43 m_options.addSwitch("Verbose","Verbose output",false);
44 }
45
46 template<short_t d,class T>
47 std::vector<bool> gsDPatchBase<d,T>::getSharpCorners(T tol) const
48 {
49 std::fill(m_vertCheck.begin(), m_vertCheck.end(), false);
50 std::vector<bool> result(m_vertCheck.size());
51
52 if (m_patches.empty())
53 {
54 gsWarn<<"Sharp corners should be detected, but no multi-patch object is assigned to this class\n";
55 return result;
56 }
57
58 std::fill(result.begin(), result.end(), false);
59 index_t cidx;
60 patchCorner pcorner;
61 for (size_t p=0; p<m_Bbases.nBases(); p++)
62 {
63 for (index_t c=1; c<5; c++)
64 {
65 cidx = _vertIndex(p,c);
66 if (m_vertCheck.at(cidx))
67 continue;
68
69 // Get the patchCorners which are on the boundary
70 pcorner = patchCorner(p,c);
71 std::vector<patchCorner> pcorners;
72 m_topology.getCornerList(pcorner,pcorners);
73 std::vector<std::pair<patchCorner,patchSide>> boundaries;
74 for (std::vector<patchCorner>::iterator it=pcorners.begin(); it!=pcorners.end(); it++)
75 {
76 std::vector<patchSide> psides;
77 it->getContainingSides(d,psides);
78 for (size_t s=0; s!=psides.size(); s++)
79 {
80 // the 0,k (k=0,1) DoF should be eliminated
81 if (!m_topology.isInterface(psides[s]))
82 {
83 boundaries.push_back(std::make_pair(*it,psides[s]));
84 continue;
85 }
86 }
87
88 // Mark the corner as passed
89 m_vertCheck[ _vertIndex(it->patch, it->corner()) ] = true;
90 }
91
92 if (boundaries.size()==0)
93 continue;
94 else if (boundaries.size()==2)
95 {
96 std::vector<gsVector<T>> normals;
97 for (std::vector<std::pair<patchCorner,patchSide>>::iterator it=boundaries.begin(); it!=boundaries.end(); it++)
98 {
99 gsMapData<T> md;
100 md.flags = NEED_OUTER_NORMAL;
101
102 // Get the coordinates of the corner
103 gsVector<bool> pars;
104 it->first.parameters_into(m_Bbases.domainDim(),pars); // get the parametric coordinates of the corner
105 gsMatrix<T> supp = m_Bbases.basis(it->first.patch).support();
106 gsVector<T> vec(supp.rows());
107 for (index_t r = 0; r!=supp.rows(); r++)
108 vec(r) = supp(r,pars(r));
109 // Assign the corner coordinates to the mapper
110 md.points = vec;
111 md.side = it->second;
112 m_patches.patch(it->first.patch).computeMap(md);
113
114 normals.push_back(md.outNormal(0).normalized());
115 }
116
117 if ( (std::abs(normals[0].transpose() * normals[1] - 1)) > tol )
118 for (std::vector<patchCorner>::iterator it=pcorners.begin(); it!=pcorners.end(); it++)
119 result[ _vertIndex(it->patch, it->corner()) ] = true;
120 }
121 else
122 GISMO_ERROR("Size of the stored boundary corners and sides must be 0 or 2, but is "<<boundaries.size());
123 }
124
125 }
126 return result;
127 }
128
129 /*=====================================================================================
130 Information functions
131 =====================================================================================*/
132
133 template<short_t d,class T>
134 void gsDPatchBase<d,T>::mapperInfo() const
135 {
136 for (size_t p=0; p!=m_bases.nBases(); p++)
137 {
138 gsInfo<<"----------------------------------\n";
139 gsInfo<<"patch "<<p<<"\n";
140 index_t size = m_mapModified.patchSize(p);
141 for (index_t k=0; k!=size; k++)
142 {
143 bool free = m_mapModified.is_free(k,p);
144 std::string str = free ? "free" : "eliminated";
145 gsInfo<<"DoF "<<k<<" is "<<str<<"\n";
146 }
147 }
148 }
149
150 template<short_t d,class T>
151 void gsDPatchBase<d,T>::vertexInfo(patchCorner corner) const
152 {
153 std::pair<index_t,bool> data = this->_vertexData(corner);
154 gsInfo<<"Patch "<<corner.patch<<", corner "<<corner<<" has valence "<<data.first<<" and is "<<(data.second ? "an interior vertex" : "a boundary vertex")<<"\n";
155
156 }
157
158 template<short_t d,class T>
159 void gsDPatchBase<d,T>::sideInfo(patchSide side) const
160 {
161 gsInfo<<"Patch "<<side.patch<<", side "<<side<<" is "<<(m_topology.isBoundary(side) ? "a boundary side" : "an interface")<<"\n";
162 }
163
164 template<short_t d,class T>
165 void gsDPatchBase<d,T>::sideInfo() const
166 {
167 gsInfo<<"**D-Patch Side info**\n";
168 for(size_t i = 0;i<m_bases.nBases();++i)
169 for(index_t j=1;j<=4;++j)
170 sideInfo(patchSide(i,j));
171 }
172
173 template<short_t d,class T>
174 void gsDPatchBase<d,T>::cornerInfo() const
175 {
176 gsInfo<<"**D-Patch Corner info**\n";
177 for(size_t i = 0;i<m_bases.nBases();++i)
178 for(index_t j=1;j<=4;++j)
179 vertexInfo(patchCorner(i,j));
180 }
181
182 /*=====================================================================================
183 Output functions
184 =====================================================================================*/
185 template<short_t d,class T>
186 gsGeometry<T>* gsDPatchBase<d,T>::exportPatch(index_t patch, bool computeCoefs)
187 {
188 GISMO_ASSERT(m_computed,"The method has not been computed! Call compute().");
190 // This can be done more efficient!!
191 // Do it once instead of for every patch
193 if (computeCoefs || m_coefs.rows()==0)
194 {
195 m_coefs = this->_preCoefficients(); // gets coefficients of the modified size
196 m_coefs = m_matrix.transpose() * m_coefs; // maps to local size
197 }
198
200 index_t size,offset = 0;
201 for (index_t p=0; p!=patch; p++)
202 offset += m_mapOriginal.patchSize(p);
203
204 size = m_mapOriginal.patchSize(patch);
205 gsMatrix<T> local = m_coefs.block(offset,0,size,m_coefs.cols());
206 return m_bases[patch].makeGeometry( give(local) ).release();
207 }
208
209 // template<short_t d,class T>
210 // gsMultiPatch<T> gsDPatchBase<d,T>::exportToPatches()
211 // {
212 // m_coefs = this->_preCoefficients();
213 // m_coefs = m_matrix.transpose() * m_coefs;
214
215 // std::vector<gsGeometry<T> *> patches(m_bases.nBases());
216 // for (size_t p=0; p!=m_bases.nBases(); p++)
217 // patches[p]= this->exportPatch(p,false);
218
219 // return gsMultiPatch<T>(patches,m_topology.boundaries(),m_topology.interfaces());
220 // }
221
222 /*=====================================================================================
223 Output functions
224 =====================================================================================*/
225
226 template<short_t d,class T>
227 index_t gsDPatchBase<d,T>::_indexFromSides(index_t index1, const patchSide side1, index_t index2, const patchSide side2)
228 {
229 /*
230 Finds the index index1 away from side1 and index2 from side2
231 index 1 is the index parallel to side 1
232 index 2 is the index parallel to side 2
233 */
234 GISMO_ASSERT(side1.patch==side2.patch,"Sides must be from the same patch");
235 GISMO_ASSERT(side1.side().direction()!=side2.side().direction(),"Sides must have different direction");
236 index_t index;
237
238 gsBasis<T> * basis = &m_bases.basis(side1.patch);
239
240 gsVector<index_t> indices1 = static_cast<gsVector<index_t>>(basis->boundaryOffset(side1.side(),index2));
241
242 index_t n = indices1.rows();
243 if (side1.side()==1) //west
244 if (side2.side().parameter()==0) //south
245 index = indices1(index1);
246 else
247 index = indices1(n-1-index1); //north
248 else if (side1.side()==2) //east
249 if (side2.side().parameter()==0) //south
250 index = indices1(index1);
251 else
252 index = indices1(n-1-index1); //north
253 else if (side1.side()==3) //south
254 if (side2.side().parameter()==0) //west
255 index = indices1(index1);
256 else
257 index = indices1(n-1-index1); //east
258 else if (side1.side()==4) //north
259 if (side2.side().parameter()==0) //west
260 index = indices1(index1);
261 else
262 index = indices1(n-1-index1); //east
263 else
264 GISMO_ERROR("Side unknown. index = "<<side1.side());
265 return index;
266 }
267
268 template<short_t d,class T>
269 index_t gsDPatchBase<d,T>::_indexFromVert(const index_t index, const patchCorner corner, const patchSide side, const index_t offset) const
270 {
271 return _indexFromVert(&m_bases.basis(corner.patch),index, corner, side, offset);
272 }
273
274 template<short_t d,class T>
275 index_t gsDPatchBase<d,T>::_indexFromVert(const gsMultiBasis<T> & bases, const index_t index, const patchCorner corner, const patchSide side, const index_t offset) const
276 {
277 return _indexFromVert(&bases.basis(corner.patch),index, corner, side, offset);
278 }
279
280 template<short_t d,class T>
281 index_t gsDPatchBase<d,T>::_indexFromVert(const gsBasis<T> * basis, const index_t index, const patchCorner corner, const patchSide side, const index_t offset) const
282 {
283 const gsTensorBSplineBasis<d,T> * tbbasis;
284 const gsTensorNurbsBasis<d,T> * tnbasis;
285 const gsHTensorBasis<d,T> * thbasis;
286
287 if ((index==0) && (offset==0))
288 {
289 // gsHTensorBasis<d,T> *basis = dynamic_cast<gsHTensorBasis<d,T>*>(&bases.basis(corner.patch));
290 index_t idx = basis->functionAtCorner(corner.corner());
291 return idx;
292 }
293 else
294 {
295 if ((tbbasis = dynamic_cast<const gsTensorBSplineBasis<d,T> * >(basis)))
296 return _indexFromVert_impl(tbbasis,index,corner,side,offset);
297 else if ((tnbasis = dynamic_cast<const gsTensorNurbsBasis<d,T> * >(basis)))
298 return _indexFromVert_impl(&tnbasis->source(),index,corner,side,offset);
299 else if ((thbasis = dynamic_cast<const gsHTensorBasis<d,T> * >(basis)))
300 return _indexFromVert_impl(thbasis,index,corner,side,offset);
301 else
302 GISMO_ERROR("Basis type unknown");
303 }
304 }
305
306 template<short_t d,class T>
307 template<class U>
308 typename util::enable_if<util::is_same<U, const gsTensorBSplineBasis<d,T> *>::value,
309 const index_t>::type
310 gsDPatchBase<d,T>::_indexFromVert_impl(U basis, const index_t index, const patchCorner corner, const patchSide side, const index_t offset) const
311 {
312 /*
313 Finds indices i in the direction of side away from the vertex
314 if index = 0, the corner index is requested
315 */
316
317 index_t result = -1;
318
319 gsVector<index_t,2> sizes;
320 basis->size_cwise(sizes);
321 if (side.side()==1) //west
322 {
323 GISMO_ASSERT(offset<sizes[0],"Offset is out of bounds");
324 GISMO_ASSERT(index<sizes[1],"Index is out of bounds");
325 if (corner.corner()==1)//southwest
326 result = basis->index(offset,index);
327 else if (corner.corner()==3) //northwest
328 result = basis->index(offset,sizes[1]-1-index);
329 else
330 GISMO_ERROR(corner.corner() << " is not adjacent to side "<<side.side()<<"!");
331 }
332 else if (side.side()==2) //east
333 {
334 GISMO_ASSERT(offset<sizes[0],"Offset is out of bounds");
335 GISMO_ASSERT(index<sizes[1],"Index is out of bounds");
336 if (corner.corner()==2)//southeast
337 result = basis->index(sizes[0]-1-offset,index);
338 else if (corner.corner()==4) //northeast
339 result = basis->index(sizes[0]-1-offset,sizes[1]-1-index);
340 else
341 GISMO_ERROR(corner.corner() << " is not adjacent to side "<<side.side()<<"!");
342 }
343 else if (side.side()==3) //south
344 {
345 GISMO_ASSERT(offset<sizes[1],"Offset is out of bounds");
346 GISMO_ASSERT(index<sizes[0],"Index is out of bounds");
347 if (corner.corner()==1)//southwest
348 result = basis->index(index,offset);
349 else if (corner.corner()==2) //southeast
350 result = basis->index(sizes[0]-1-index,offset);
351 else
352 GISMO_ERROR(corner.corner() << " is not adjacent to side "<<side.side()<<"!");
353 }
354 else if (side.side()==4) //north
355 {
356 GISMO_ASSERT(offset<sizes[1],"Offset is out of bounds");
357 GISMO_ASSERT(index<sizes[0],"Index is out of bounds");
358 if (corner.corner()==3)//northwest
359 result = basis->index(index,sizes[1]-1-offset);
360 else if (corner.corner()==4) //northeast
361 result = basis->index(sizes[0]-1-index,sizes[1]-1-offset);
362 else
363 GISMO_ERROR(corner.corner() << " is not adjacent to side "<<side.side()<<"!");
364
365 }
366
367 return result;
368 }
369
370 template<short_t d,class T>
371 template<class U>
372 typename util::enable_if<util::is_same<U, const gsHTensorBasis<d,T> *>::value,
373 const index_t>::type
374 gsDPatchBase<d,T>::_indexFromVert_impl(U basis, const index_t index, const patchCorner corner, const patchSide side, const index_t offset) const
375 {
376 /*
377 Finds indices i in the direction of side away from the vertex
378 if index = 0, the corner index is requested
379 */
380
381 index_t result = -1;
382 index_t level = basis->levelAtCorner(corner.corner());
383 gsTensorBSplineBasis<d,T> tbasis;
384 tbasis = basis->tensorLevel(level);
385 result = _indexFromVert(&tbasis,index,corner,side,offset);
386 result = basis->flatTensorIndexToHierachicalIndex(result,level);
387 if (result!=-1) return result;
388
389 GISMO_ASSERT(basis->maxLevel()!=0,"Basis must have more than one level!");
390 GISMO_ASSERT(level!=0,"The level in the corner should not be 0!");
391
392 // Find last active index
393 index_t nActive_f = 0;
394 result = 0;
395 while (result!=-1)
396 {
397 result = _indexFromVert(&tbasis,nActive_f++,corner,side,offset);
398 result = basis->flatTensorIndexToHierachicalIndex(result,level);
399 }
400 nActive_f--;
401 // nActive_f is the number of functions from the corner that are active in the finest level
402
403 // Now, we find the number of functions that are not active in the basis one level coarser
404 index_t nInactive_c = 0;
405 level--;
406 result = -1;
407 tbasis = basis->tensorLevel(level);
408 while (result==-1) //&& nInactive_c < MAX!!
409 {
410 result = _indexFromVert(&tbasis,nInactive_c++,corner,side,offset);
411 result = basis->flatTensorIndexToHierachicalIndex(result,level);
412 }
413 nInactive_c--;
414 // nInactive_c number of functions from the corner that are inactive one level coarser
415
416 // Now, the distance of the function with nInactive_c is with respect to the corner, is nActive_f
417 index_t tmpIdx = index - nActive_f + nInactive_c;
418 result = _indexFromVert(&tbasis,tmpIdx,corner,side,offset);
419 result = basis->flatTensorIndexToHierachicalIndex(result,level);
420 GISMO_ASSERT(result!=-1,"Something went wrong");
421 return result;
422 }
423
424 template<short_t d,class T>
425 const std::pair<index_t,bool> gsDPatchBase<d,T>::_vertexData(const patchCorner corner) const
426 {
427 std::vector<patchCorner> corners;
428 std::pair<index_t,bool> output;
429 output.second = m_topology.getCornerList(corner,corners); // bool is true if interior vertex
430 output.first = corners.size();
431 return output;
432 }
433
434 template<short_t d,class T>
435 void gsDPatchBase<d,T>::_getLowestCorners(std::vector<patchCorner> & pcorners, index_t n) const
436 {
437 struct {
438 bool operator()(patchCorner a, patchCorner b) const { return a.patch < b.patch; }
439 } customLess;
440
441 // Sort
442 std::sort(pcorners.begin(), pcorners.end(),customLess);
443 // Resize; this vector are the indices we want to keep the DoFs of
444 pcorners.resize(n);
445 }
446
447 template<short_t d,class T>
448 void gsDPatchBase<d,T>::_removeLowestCorners(std::vector<patchCorner> & pcorners, index_t n) const
449 {
450 index_t size = pcorners.size();
451 GISMO_ASSERT(n<=size,"You cannot remove more corners than there are actually stored in the container. Container size = "<<size<<" no. corners to be removed = "<<n);
452 struct {
453 bool operator()(patchCorner a, patchCorner b) const { return a.patch > b.patch; }
454 } customGreater;
455
456 // Sort
457 std::sort(pcorners.begin(), pcorners.end(),customGreater);
458 // Resize; this vector are the indices we want to keep the DoFs of
459 pcorners.resize(size-n);
460 }
461
462 template<short_t d,class T>
463 void gsDPatchBase<d,T>::_getLowestIndices(std::vector<std::pair<index_t,index_t>> & indices, index_t n) const
464 {
465 // indices are pairs with first=patchID, second=localID
466 struct {
467 bool operator()(std::pair<index_t,index_t> a, std::pair<index_t,index_t> b) const
468 {
469 return
470 (a.first < b.first)
471 ||
472 ((a.first == b.first) &&
473 (a.second < b.second) );
474 }
475 } customLess;
476
477 // Sort
478 std::sort(indices.begin(), indices.end(),customLess);
479 // Resize; this vector are the indices we want to keep the DoFs of
480 indices.resize(n);
481 }
482
483 template<short_t d,class T>
484 void gsDPatchBase<d,T>::_removeLowestIndices(std::vector<std::pair<index_t,index_t>> & indices, index_t n) const
485 {
486 index_t size = indices.size();
487 GISMO_ASSERT(n<=size,"You cannot remove more corners than there are actually stored in the container. Container size = "<<size<<" no. corners to be removed = "<<n);
488 // indices are pairs with first=patchID, second=localID
489 struct {
490 bool operator()(std::pair<index_t,index_t> a, std::pair<index_t,index_t> b) const
491 {
492 return
493 (a.first > b.first)
494 ||
495 ((a.first == b.first) &&
496 (a.second > b.second) );
497 }
498 } customGreater;
499
500 // Sort
501 std::sort(indices.begin(), indices.end(),customGreater);
502 // Resize; this vector are the indices we want to keep the DoFs of
503 indices.resize(size-n);
504 }
505
506 // gets all interface DoFs on \a depth from the corner
507 template<short_t d,class T>
508 std::vector<std::pair<index_t,index_t>> gsDPatchBase<d,T>::_getInterfaceIndices(patchCorner pcorner, index_t depth, const gsMultiBasis<T> & mbasis) const
509 {
510 std::vector<std::pair<index_t,index_t>> result;
511 std::vector<patchSide> csides(d);
512 index_t index, patch;
513
514 pcorner.getContainingSides(d,csides);
515 patch = pcorner.patch;
516 const gsBasis<T> * basis = &mbasis.basis(patch);
517
518 if (depth==0)
519 {
520 // add the 0,0 one
521 index = basis->functionAtCorner(pcorner);
522 result.push_back(std::make_pair(patch,index));
523 }
524 else
525 {
526 for (index_t i = 0; i!=d; i++)
527 {
528 if ( m_topology.isBoundary(csides[i]) ) continue;
529 index = _indexFromVert(mbasis,depth,pcorner,csides[i],0);
530 result.push_back(std::make_pair(patch,index));
531 }
532 }
533 return result;
534 }
535
536 // gets all interface DoFs on \a depth from the all corners surrounding pcorner
537 template<short_t d,class T>
538 std::vector<std::pair<index_t,index_t>> gsDPatchBase<d,T>::_getAllInterfaceIndices(patchCorner pcorner, index_t depth, const gsMultiBasis<T> & mbasis) const
539 {
540 std::vector<std::pair<index_t,index_t>> result, tmp;
541 std::vector<patchCorner> pcorners;
542
543 m_topology.getCornerList(pcorner,pcorners); // bool is true if interior vertex
544 for (std::vector<patchCorner>::const_iterator it=pcorners.begin(); it!=pcorners.end(); it++)
545 {
546 tmp = this->_getInterfaceIndices(*it,depth,mbasis);
547 result.insert(result.end(), tmp.begin(), tmp.end());
548 }
549 return result;
550 }
551
552 template<short_t d,class T>
553 bool gsDPatchBase<d,T>::_checkMatrix(const gsSparseMatrix<T> & matrix) const // ! makes a deep copy (otherwise the contents of m_matrix get destroyed somehow...)
554 {
555 GISMO_ASSERT(matrix.cols()==matrix.outerSize(),"is the matrix ColMajor?");
556 gsVector<T> colSums(matrix.cols());
557 colSums.setZero();
558 for (index_t i = 0; i<matrix.outerSize(); ++i)
559 for (typename gsSparseMatrix<T>::iterator it = matrix.begin(i); it; ++it)
560 colSums.at(i) += it.value();
561
562 // TODO: Make tolerance depending on real_t
563 colSums.array() -= 1;
564 return (colSums.array() < 1e-8).all() && (colSums.array() > -1e-8).all(); // Lower?
565 }
566
567// /*=====================================================================================
568// Construction functions
569// =====================================================================================*/
570
571
572// template<short_t d,class T>
573// void gsDPatchBase<d,T>::_makeTHB()
574// {
575// // prepare the geometry
576// std::vector<std::vector<patchCorner> > cornerLists;
577// // m_RefPatches.getEVs(cornerLists,true);
578
579// // get the corners that need refinement
580// std::vector<patchCorner> cornerList;
581// patchCorner pcorner;
582// index_t cidx;
583// for(size_t p = 0;p<m_RefPatches.nPatches();++p)
584// {
585// for(index_t c=1;c<=4;++c)
586// {
587// pcorner=patchCorner(p,c);
588// cidx = _vertIndex(p,c);
589// bool C0 = m_C0s[cidx];
590// bool isCycle = m_RefPatches.getCornerList(pcorner,cornerList);
591// bool alreadyReached = false;
592// for(size_t k = 0;k<cornerList.size();++k)
593// if((size_t)cornerList[k].patch<p)
594// alreadyReached = true;
595
596// // add if
597// // interior vertex with valence!=4
598// // or
599// // boundary vertex with valence > 2 (unless C0, then valence > 1)
600// if(((isCycle&&cornerList.size()!=4)||((!isCycle)&&cornerList.size()>2-C0))&&!alreadyReached)
601// cornerLists.push_back(cornerList);
602// }
603// }
604
605// if (cornerLists.size()!=0)
606// {
607// /// Change the coefficients
608// gsMatrix<T> coefs = this->freeCoefficients(); // gets coefficients of the modified size
609// coefs = m_matrix.transpose() * coefs; // maps to local size
610
611// this->setCoefficients(coefs,m_RefPatches);
612
613// /// Handle the EVs
614// std::vector< std::vector<index_t> > elVec(m_RefPatches.nPatches());
615// for (size_t v =0; v!=cornerLists.size(); v++)
616// for (size_t c = 0; c!=cornerLists[v].size(); c++)
617// {
618// patchCorner corner = cornerLists[v].at(c);
619// gsVector<bool> pars;
620// corner.corner().parameters_into(m_RefPatches.parDim(),pars); // get the parametric coordinates of the corner
621// gsMatrix<T> mat = pars.template cast<T>(); // cast to real coordinates
622
623// gsMatrix<T> boxes(m_RefPatches.parDim(),2);
624// boxes.col(0) << mat;
625// boxes.col(1) << mat;
626
627// gsHTensorBasis<d,T> *basis = dynamic_cast<gsHTensorBasis<d,T>*>(&m_RefPatches.basis(corner.patch));
628// std::vector<index_t> elements = basis->asElements(boxes,0); // 0-ring
629
630// elVec.at(corner.patch).insert(elVec.at(corner.patch).end(), elements.begin(), elements.end());
631
632// // gsHTensorBasis<d,T> *basis = dynamic_cast<gsHTensorBasis<d,T>*>(&m_RefPatches.basis(corner.patch));
633
634// // basis->refineElements(elements, m_tMatrix);
635// }
636
637// gsSparseMatrix<T> tmp;
638// index_t rows = 0, cols = 0;
639// std::vector<gsEigen::Triplet<T,index_t>> tripletList;
640// for (size_t p=0; p!=m_RefPatches.nPatches(); p++)
641// {
642// gsHTensorBasis<d,T> *basis = dynamic_cast<gsHTensorBasis<d,T>*>(&m_RefPatches.basis(p));
643// std::vector< gsSortedVector< index_t > > xmat = basis->getXmatrix();
644
645// m_RefPatches.patch(p).refineElements(elVec[p]);
646
647// basis->transfer(xmat,tmp);
648
649// for (index_t i = 0; i<tmp.outerSize(); ++i)
650// for (typename gsSparseMatrix<T>::iterator it(tmp,i); it; ++it)
651// tripletList.push_back(gsEigen::Triplet<T,index_t>(it.row()+rows,it.col()+cols,it.value()));
652
653// rows += tmp.rows();
654// cols += tmp.cols();
655// }
656
657// m_tMatrix.resize(rows,cols);
658// m_tMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
659
660// m_tMatrix.makeCompressed();
661// m_bases = gsMultiBasis<T>(m_RefPatches);
662// }
663
664// // redefine the mappers
665// m_mapOriginal = gsDofMapper(m_bases);
666// m_mapOriginal.finalize();
667
668// // gsWriteParaview<T>(m_RefPatches,"mp_ref",1000,true);
669// }
670
671// template<short_t d,class T>
672// void gsDPatchBase<d,T>::_computeEVs()
673// {
674// /*
675// Our goal is to create three vectors c11, c12, c21 which all contain the
676// c11, c12 and c21 coefficients of the patches around the EV in the right order
677// (counter)-clockwise.
678// */
679
680// std::vector<std::vector<patchCorner> > cornerLists;
681// // m_topology.getEVs(cornerLists);
682
683// // get the corners that need refinement
684// std::vector<patchCorner> cornerList;
685// patchCorner pcorner;
686// index_t cidx;
687// for(size_t p = 0;p<m_RefPatches.nPatches();++p)
688// {
689// for(index_t c=1;c<=4;++c)
690// {
691// pcorner=patchCorner(p,c);
692// cidx = _vertIndex(p,c);
693// bool C0 = m_C0s[cidx];
694// bool isCycle = m_RefPatches.getCornerList(pcorner,cornerList);
695// bool alreadyReached = false;
696// for(size_t k = 0;k<cornerList.size();++k)
697// if((size_t)cornerList[k].patch<p)
698// alreadyReached = true;
699
700// // add if
701// // interior vertex with valence!=4
702// // or
703// // boundary vertex with valence > 2 (unless C0, then valence > 1)
704// if(((isCycle&&cornerList.size()!=4)||((!isCycle)&&cornerList.size()>2-C0))&&!alreadyReached)
705// cornerLists.push_back(cornerList);
706// }
707// }
708
709// std::fill(m_vertCheck.begin(), m_vertCheck.end(), false);
710// if (cornerLists.size()!=0)
711// {
712// m_matrix = m_matrix * m_tMatrix.transpose();
713
714// std::vector<patchCorner> pcorners;
715// patchCorner pcorner;
716// gsMatrix<T> Cg; // coefficients
717// gsMatrix<T> ub; // baricentric coordinates
718// gsMatrix<index_t> uind; // column indices of baricentric coordinates
719// index_t cidx;
720
721// for (std::vector<std::vector<patchCorner> >::iterator it=cornerLists.begin(); it!=cornerLists.end(); it++)
722// {
723
724// std::vector<patchCorner> pcorners = *it;
725// pcorner = it->at(0);
726// cidx = _vertIndex(pcorner.patch,pcorner.corner());
727// if (m_vertCheck.at(cidx))
728// continue;
729
730// std::pair<index_t,bool> vdata = _vertexData(pcorner); // corner c
731
732// // get the triangle
733// gsMatrix<T> Cg;
734// std::tie(Cg,ub,uind) = _makeTriangle(pcorner);
735
736// // The barycentric coordinates will be attached to the matrix rows corresponding to the 0,0 corners (and the three lowest patch index corners whenever valence > 3)
737// // We use _getLowestIndices such that the corners are assigned to increasing patch corners
738// // We need the index on the old basis (the unrefined basis), because we plug it into the mapModified (which maps the local DoFs to the global ones)
739// std::vector<std::pair<index_t,index_t>> indices, tmp;
740// if (vdata.first==2)
741// {
742// // These are two indices
743// indices = _getAllInterfaceIndices(pcorner,0,m_Bbases);
744// tmp = _getAllInterfaceIndices(pcorner,1,m_Bbases);
745// _getLowestIndices(tmp,1);
746// indices.push_back(tmp[0]);
747// }
748// else
749// {
750// indices = _getAllInterfaceIndices(pcorner,0,m_Bbases);
751// _getLowestIndices(indices,3);
752// }
753
754
755// std::vector<index_t> rowIndices;
756// rowIndices.reserve(3);
757// for (std::vector<std::pair<index_t,index_t>>::iterator it=indices.begin(); it!=indices.end(); it++)
758// {
759// // We need the index on the old basis (the unrefined basis), because we plug it into the mapModified (which maps the local DoFs to the global ones)
760// GISMO_ASSERT(m_mapModified.is_free(it->second,it->first),"This DoF must be free! patch = "<<it->first<<"; index = "<<it->first);
761// rowIndices.push_back(m_mapModified.index(it->second,it->first));
762// }
763
764// index_t rowIdx,colIdx;
765// // set the colums related to the barycentric columns equal to zero
766// for (index_t j=0; j!=ub.cols(); j++)
767// {
768// colIdx = uind(0,j);
769// m_matrix.prune(
770// [&colIdx](index_t i, index_t j, T)
771// { return j!=colIdx; }
772// );
773// }
774
775// for (index_t i=0; i!=ub.rows(); i++)
776// for (index_t j=0; j!=ub.cols(); j++)
777// {
778// rowIdx = rowIndices[i];
779// colIdx = uind(0,j);
780// m_matrix(rowIdx,colIdx) = ub(i,j);
781// }
782
783// for (size_t k = 0; k!=pcorners.size(); k++)
784// m_vertCheck[ _vertIndex(pcorners[k].patch, pcorners[k].corner()) ] = true;
785// }
786// m_matrix.makeCompressed();
787// }
788// }
789//
790
791 template<short_t d,class T>
792 void gsDPatchBase<d,T>::_initialize() // also initialize the mappers!
793 {
794 _initChecks();
795
796 m_C0s.resize(m_nVerts);
797 if (m_options.getSwitch("SharpCorners"))
798 m_C0s = getSharpCorners(m_options.getReal("SharpCornerTolerance"));
799 else
800 std::fill(m_C0s.begin(), m_C0s.end(), false);
801
802 _initTHB();
803 _initBasis();
804 _countDoFs();
805 _initMappers();
806 _initMatrix();
807 }
808
809 template<short_t d,class T>
810 void gsDPatchBase<d,T>::_initChecks()
811 {
812 m_topology.checkConsistency();
813 m_nSides = 2*m_topology.nInterfaces() + m_topology.nBoundary();
814 m_nVerts = 4*m_Bbases.nBases();
815
816 m_sideCheck.resize(m_nSides);
817 std::fill(m_sideCheck.begin(), m_sideCheck.end(), false);
818 m_vertCheck.resize(m_nVerts);
819 std::fill(m_vertCheck.begin(), m_vertCheck.end(), false);
820 }
821
822 template<short_t d,class T>
823 void gsDPatchBase<d,T>::_initTHB()
824 {
825 // Cast all patches of the mp object to THB splines
826 for (size_t k=0; k!=m_Bbases.nBases(); ++k)
827 {
828 if ( (dynamic_cast<const gsTensorBSplineBasis<d,T> * > (&m_Bbases.basis(k))) )
829 m_bases.addBasis(new gsTHBSplineBasis<d,T>(m_Bbases.basis(k)));
830 else if ((dynamic_cast<const gsTHBSplineBasis<d,T> * > (&m_Bbases.basis(k))))
831 m_bases.addBasis(m_Bbases.basis(k).clone());
832 else
833 gsWarn<<"No THB basis was constructed";
834 }
835 m_bases.setTopology(m_topology);
836 }
837
838 template<short_t d,class T>
842
843 template<short_t d,class T>
844 void gsDPatchBase<d,T>::_initMappers()
845 {
846 m_mapModified = gsDofMapper(m_bases);
847 m_mapOriginal = gsDofMapper(m_bases);
848 }
849
850 template<short_t d,class T>
851 void gsDPatchBase<d,T>::_initMatrix()
852 {
853 m_matrix.resize(m_size,m_bases.totalSize());
854 }
855
856 template<short_t d,class T>
857 void gsDPatchBase<d,T>::_whichHandled()
858 {
859 for (size_t p=0; p!=m_bases.nBases(); p++)
860 {
861 gsInfo<<"----------------------------------\n";
862 gsInfo<<"patch "<<p<<"\n";
863 for (size_t b=0; b!=m_mapOriginal.patchSize(p); b++)
864 {
865 index_t idx = m_mapModified.index(b,p);
866 if (m_mapModified.is_free_index(idx))
867 gsInfo<<"basis function "<<b<<" (check="<<m_basisCheck[idx]<<") on patch "<<p<<" is "<<(m_basisCheck[idx] ? "":"not ")<<"handled\n";
868 else
869 gsInfo<<"basis function "<<b<<" on patch "<<p<<" is "<<"eliminated\n";
870 }
871 }
872 }
873
874 template<short_t d,class T>
875 void gsDPatchBase<d,T>::_performChecks(bool basis)
876 {
877 bool checkSides = std::all_of(m_sideCheck.begin(), m_sideCheck.end(), [](bool m_sideCheck) { return m_sideCheck; });
878 GISMO_ENSURE(checkSides,"Not all sides are checked");
879 bool checkVerts = std::all_of(m_vertCheck.begin(), m_vertCheck.end(), [](bool m_vertCheck) { return m_vertCheck; });
880 GISMO_ENSURE(checkVerts,"Not all vertices are checked");
881
882 if (!basis)
883 return;
884
885 bool checkBasis = std::all_of(m_basisCheck.begin(), m_basisCheck.end(), [](bool m_basisCheck) { return m_basisCheck; });
886 GISMO_ENSURE(checkBasis,"Not all basis functions are checked");
887 }
888
889 template<short_t d,class T>
890 void gsDPatchBase<d,T>::_resetChecks(bool basis)
891 {
892 m_sideCheck.resize(m_nSides);
893 std::fill(m_sideCheck.begin(), m_sideCheck.end(), false);
894 m_vertCheck.resize(m_nVerts);
895 std::fill(m_vertCheck.begin(), m_vertCheck.end(), false);
896
897 if (!basis)
898 return;
899
900 m_basisCheck.resize(m_size);
901 std::fill(m_basisCheck.begin(), m_basisCheck.end(), false);
902
903 }
904
905 template<short_t d,class T>
906 void gsDPatchBase<d,T>::_handleInteriorVertex(patchCorner pcorner, index_t valence)
907 {
908 std::vector<patchSide> psides;
909 std::vector<patchCorner> corners;
910 sparseEntry_t entries;
911
912 boundaryInterface iface;
913 patchSide otherSide;
914
915 index_t colIdx, rowIdx;
916 T weight;
917
918 pcorner.getContainingSides(d,psides);
919
920 index_t b11_p1 = _indexFromVert(1,pcorner,psides[0],1); // point 1,1 (does not matter which reference side is taken)
921 rowIdx = m_mapModified.index(b11_p1,pcorner.patch);
922 colIdx = m_mapOriginal.index(b11_p1,pcorner.patch);
923 // Influence of 1,1 to itself
924 weight = 1.;
925 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
926
927 m_topology.getCornerList(pcorner,corners);
928
929 // The 1,1 corners of each patch will be given 0.5 weight in the interface handling, but in addition they will receive a 1/v weight from the 0,0 DoFs on each patch
930 gsBasis<T> * tmpBasis;
931 index_t index;
932 for (std::vector<patchCorner>::iterator corn = corners.begin(); corn != corners.end(); ++corn)
933 {
934 tmpBasis = &m_bases.basis(corn->patch);
935 index = tmpBasis->functionAtCorner(corn->corner());
936 colIdx = m_mapOriginal.index(index,corn->patch);
937 weight = 1./valence;
938 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
939
940 // m_matrix(rowIdx,colIdx) = 1./valence;
941 }
942
943 // Interface contributions
944 for (std::vector<patchSide>::iterator side = psides.begin(); side != psides.end(); ++side)
945 {
946 GISMO_ENSURE(m_topology.getInterface(*side,iface),"Side must be an interface!");
947 m_topology.getNeighbour(*side,otherSide);
948 patchCorner otherCorner = iface.mapCorner(pcorner);
949
950 index_t b10_p1 = _indexFromVert(1,pcorner,*side,0); // index from vertex pcorners[c] along side psides[0] with offset 0.
951 index_t b10_p2 = _indexFromVert(1,otherCorner,otherSide,0); // point 0,1
952
953 weight = 0.5;
954 colIdx = m_mapOriginal.index(b10_p1,side->patch);
955 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
956 colIdx = m_mapOriginal.index(b10_p2,otherSide.patch);
957 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
958 }
959
960
961 #pragma omp critical (handle_interior_vertex)
962 {
963 _pushAndCheck(entries);
964 m_vertCheck[ _vertIndex(pcorner.patch, pcorner.corner()) ] = true;
965 }
966 }
967
968
969 template<short_t d,class T>
970 void gsDPatchBase<d,T>::_handleVertex(patchCorner pcorner)
971 {
972
973 if (m_vertCheck[ _vertIndex(pcorner.patch, pcorner.corner()) ])
974 {
975 // gsDebug<<"corner "<<pcorner.corner()<<" ("<<pcorner.patch<<") skipped!\n";
976 return;
977 }
978
979 std::pair<index_t,bool> vdata = _vertexData(pcorner); // corner c
980
981 bool C0 = m_C0s[_vertIndex(pcorner.patch,pcorner.corner())];
982 if (!vdata.second) // boundary vertices
983 {
984 if (vdata.first==1)
985 _handleRegularCorner(pcorner);
986 else if (vdata.first==2 && !C0)
987 _handleRegularBoundaryVertexSmooth(pcorner,vdata.first);
988 else if (vdata.first==2 && C0)
989 _handleRegularBoundaryVertexNonSmooth(pcorner,vdata.first);
990 else if (vdata.first> 2 && !C0)
991 _handleIrregularBoundaryVertexSmooth(pcorner,vdata.first);
992 else if (vdata.first> 2 && C0)
993 _handleIrregularBoundaryVertexNonSmooth(pcorner,vdata.first);
994 else
995 GISMO_ERROR("Something went wrong");
996 } // end boundary vertices
997 else // interior vertices
998 _handleInteriorVertex(pcorner,vdata.first);
999 }
1000
1001 template<short_t d,class T>
1002 void gsDPatchBase<d,T>::_handleInterface(boundaryInterface iface)
1003 {
1004 if (m_sideCheck[ _sideIndex(iface.first().patch, iface.first().side()) ] || m_sideCheck[ _sideIndex(iface.second().patch, iface.second().side()) ])
1005 {
1006 // gsDebug<<"sides "<<iface.first().side()<<" ("<<iface.first().patch<<") and "<<iface.second().side()<<" ("<<iface.second().patch<<") skipped!\n";
1007 return;
1008 }
1009
1010 std::vector<patchCorner> pcorners;
1011
1012 std::vector<std::vector<index_t>> selectedIndices(2);
1013 std::vector<std::vector<index_t>> selectedOIndices(2);
1014
1015 std::vector<gsBasis<T> *> basis(2);
1016 std::vector<gsMatrix<index_t>> indices(2); // interface indices
1017 std::vector<gsMatrix<index_t>> oindices(2); // interface indices
1018
1019 sparseEntry_t entries;
1020
1021 // get the bases belonging to both patches
1022 basis[0] = &m_bases.basis(iface.first().patch);
1023 basis[1] = &m_bases.basis(iface.second().patch);
1024
1025 // this assumes the directions are handled correctly in matchWith (indices has the same direction as oindices)
1026 basis[0]->matchWith(iface,*basis[1],indices[0],indices[1],0);
1027 basis[0]->matchWith(iface,*basis[1],oindices[0],oindices[1],1);
1028
1029 index_t np;
1030 for (index_t p =0; p!=2; p++)
1031 {
1032 patchSide ifacep = (p==0) ? iface.first() : iface.second();
1033 ifacep.getContainedCorners(d,pcorners);
1034 for (index_t c =0; c!=2; c++)
1035 {
1036 selectedIndices[p].push_back(_indexFromVert(0,pcorners[c],ifacep,0)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1037 selectedIndices[p].push_back(_indexFromVert(1,pcorners[c],ifacep,0)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1038
1039 selectedOIndices[p].push_back(_indexFromVert(0,pcorners[c],ifacep,1)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1040 selectedOIndices[p].push_back(_indexFromVert(1,pcorners[c],ifacep,1)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1041 }
1042 std::vector<index_t> allIndices(indices[p].data(), indices[p].data() + indices[p].rows() * indices[p].cols());
1043 std::vector<index_t> result;
1044 std::copy_if(allIndices.begin(), allIndices.end(), std::back_inserter(result),
1045 [&selectedIndices,&p] (index_t entry)
1046 {
1047 std::vector<index_t>::const_iterator res = std::find(selectedIndices[p].begin(), selectedIndices[p].end(), entry);
1048 return (res == selectedIndices[p].end());
1049 });
1050 indices[p] = gsAsMatrix<index_t>(result);
1051
1052 std::vector<index_t> allOIndices(oindices[p].data(), oindices[p].data() + oindices[p].rows() * oindices[p].cols());
1053 result.clear();
1054 std::copy_if(allOIndices.begin(), allOIndices.end(), std::back_inserter(result),
1055 [&selectedOIndices,&p] (index_t entry)
1056 {
1057 std::vector<index_t>::const_iterator res = std::find(selectedOIndices[p].begin(), selectedOIndices[p].end(), entry);
1058 return (res == selectedOIndices[p].end());
1059 });
1060 oindices[p] = gsAsMatrix<index_t>(result);
1061 }
1062
1063 GISMO_ASSERT(indices[0].size()==indices[1].size(),"Indices do not have the right size, indices[0].size()="<<indices[0].size()<<",indices[1].size()="<<indices[1].size());
1064 GISMO_ASSERT(oindices[0].size()==oindices[1].size(),"Offset indices do not have the right size, oindices[0].size()="<<oindices[0].size()<<",oindices[1].size()="<<oindices[1].size());
1065
1066 index_t rowIdx,colIdx;
1067 T weight;
1068 // loop over adjacent patches and couple the DoFs.
1069 for (index_t p =0; p!= 2; p++)
1070 {
1071 np = 1-p; // not index p;
1072 patchSide ifacep = (p==0) ? iface.first() : iface.second();
1073 patchSide ifacenp= (np==0) ? iface.first() : iface.second();
1074 for (index_t k=0; k!= indices[p].size(); k++ )
1075 {
1076 GISMO_ASSERT(m_mapModified.is_free(oindices[p].at(k),ifacep.patch),"Index "<<oindices[p].at(k)<<" on patch "<<ifacep.patch<<" is eliminated. Something went wrong?");
1077 rowIdx = m_mapModified.index(oindices[p].at(k),ifacep.patch);
1078 // rowIdx1 = m_mapOriginal.index(oindices[p].at(k),patches[p]);
1079 GISMO_ASSERT(m_mapOriginal.is_free(oindices[p].at(k),ifacep.patch),"Index is eliminated. Something went wrong?");
1080 colIdx = m_mapOriginal.index(oindices[p].at(k),ifacep.patch);
1081 // m_matrix(rowIdx,colIdx) = 1.0;
1082 weight = 1.0;
1083 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1084
1085 GISMO_ASSERT(m_mapOriginal.is_free(indices[p].at(k),ifacep.patch),"Index is eliminated. Something went wrong?");
1086 colIdx = m_mapOriginal.index(indices[p].at(k),ifacep.patch);
1087 // m_matrix(rowIdx,colIdx) = 0.5;
1088 weight = 0.5;
1089 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1090
1091 GISMO_ASSERT(m_mapOriginal.is_free(indices[np].at(k),ifacenp.patch),"Index is eliminated. Something went wrong?");
1092 colIdx = m_mapOriginal.index(indices[np].at(k),ifacenp.patch);
1093 weight = 0.5;
1094 // m_matrix(rowIdx,colIdx) = 0.5;
1095 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1096
1097 // m_basisCheck[rowIdx] = true;
1098 }
1099 // m_sideCheck[ _sideIndex(iface[p].patch, iface[p].side()) ] = true; // side finished
1100 }
1101
1102 #pragma omp critical (handle_interface)
1103 {
1104 _pushAndCheck(entries);
1105 m_sideCheck[ _sideIndex(iface.first().patch, iface.first().side()) ] = true; // side finished
1106 m_sideCheck[ _sideIndex(iface.second().patch, iface.second().side()) ] = true; // side finished
1107 }
1108
1109 }
1110
1111 template<short_t d,class T>
1112 void gsDPatchBase<d,T>::_handleBoundary(patchSide side)
1113 {
1114 std::vector<patchCorner> pcorners;
1115 std::vector<index_t> selectedIndices;
1116 gsBasis<T> * basis = &m_bases.basis(side.patch);
1117 sparseEntry_t entries;
1118
1119 gsMatrix<index_t> indices = basis->boundaryOffset(side.side(),0);
1120 side.getContainedCorners(d,pcorners);
1121 for (index_t c =0; c!=2; c++)
1122 {
1123 selectedIndices.push_back(_indexFromVert(0,pcorners[c],side,0)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1124 selectedIndices.push_back(_indexFromVert(1,pcorners[c],side,0)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1125 }
1126
1127 std::sort(selectedIndices.begin(),selectedIndices.end());
1128 std::vector<index_t> allIndices(indices.data(), indices.data() + indices.rows() * indices.cols());
1129 std::vector<index_t> result(allIndices.size());
1130 std::vector<index_t>::iterator it=std::set_difference (allIndices.begin(), allIndices.end(), selectedIndices.begin(), selectedIndices.end(), result.begin());
1131 result.resize(it-result.begin());
1132
1133 index_t rowIdx,colIdx;
1134 T weight;
1135 for (std::vector<index_t>::iterator it = result.begin(); it!=result.end(); ++it)
1136 {
1137 rowIdx = m_mapModified.index(*it,side.patch);
1138 colIdx = m_mapOriginal.index(*it,side.patch);
1139 weight = 1.0;
1140 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1141
1142 // m_matrix(rowIdx,colIdx) = 1.0;
1143 // m_basisCheck[rowIdx] = true;
1144 }
1145
1146 #pragma omp critical (handle_boundary)
1147 {
1148 _pushAndCheck(entries);
1149
1150 m_sideCheck.at( _sideIndex(side.patch,side.side()) ) = true;
1151 }
1152 }
1153
1154 template<short_t d,class T>
1155 void gsDPatchBase<d,T>::_handleInterior()
1156 {
1157 #pragma omp critical (handle_interior)
1158 {
1159 index_t rowIdx,colIdx;
1160 for(size_t p=0; p!=m_bases.nBases(); p++)
1161 for (index_t b=0; b!=m_bases.basis(p).size(); b++)
1162 {
1163 rowIdx = m_mapModified.index(b,p);
1164 // rowIdx = m_mapOriginal.index(b,p);
1165 if ( (!m_mapModified.is_free(b,p)) || (m_basisCheck[rowIdx]) )
1166 // if ( (m_basisCheck[rowIdx]) )
1167 continue;
1168 colIdx = m_mapOriginal.index(b,p);
1169 m_matrix(rowIdx,colIdx) = 1;
1170 m_basisCheck[rowIdx] = true;
1171 // gsInfo<<"Basis function "<<rowIdx<<"(patch: "<<p<<"; fun: "<<b<<") is "<< (m_basisCheck[rowIdx] ? "" : "not ")<<"processed\n";
1172 }
1173 }
1174 }
1175
1176 template<short_t d,class T>
1177 void gsDPatchBase<d,T>::_computeSmoothMatrix()
1178 {
1179 GISMO_ASSERT((size_t)m_mapModified.freeSize()==m_size,"Size does not match predicted size, m_mapModified.freeSize()="<<m_mapModified.freeSize()<<"; m_size="<<m_size);
1180
1181 _resetChecks(true);
1182
1183 // iterate over the vertices
1184// #pragma omp parallel
1185// {
1186 // #pragma omp parallel for collapse(2)
1187 for (size_t p=0; p<m_bases.nBases(); p++)
1188 for (index_t c=1; c<5; c++)
1189 _handleVertex(patchCorner(p,c));
1190
1191 // #pragma omp parallel for
1192 for(gsBoxTopology::const_iiterator iit = m_topology.iBegin(); iit< m_topology.iEnd(); iit++)
1193 _handleInterface(*iit);
1194
1195 // boundaries
1196 // #pragma omp parallel for
1197 for(gsBoxTopology::const_biterator bit = m_topology.bBegin(); bit< m_topology.bEnd(); bit++)
1198 _handleBoundary(*bit);
1199
1200 _handleInterior();
1201// }
1202 if (m_options.getSwitch("Verbose")) { _whichHandled(); }
1203
1204 _performChecks(true);
1205 m_matrix.makeCompressed();
1206 }
1207
1208 template<short_t d,class T>
1209 void gsDPatchBase<d,T>::_computeMapper() // also initialize the mappers!
1210 {
1211 // interfaces
1212 _resetChecks(false);
1213
1214 patchCorner pcorner;
1215
1216 // For the interfaces, we eliminate all DoFs located on the interface, except the ones coinciding with the end vertices
1217// #pragma omp parallel
1218// {
1219// #pragma omp parallel for
1220 for(gsBoxTopology::const_iiterator iit = m_topology.iBegin(); iit!= m_topology.iEnd(); iit++)
1221 _computeInterfaceMapper(*iit);
1222// }
1223 // On the boundaries, we don't do anything
1224// #pragma omp parallelz
1225// {
1226 // #pragma omp parallel for
1227 for(gsBoxTopology::const_biterator bit = m_topology.bBegin(); bit!= m_topology.bEnd(); bit++)
1228 _computeBoundaryMapper(*bit);
1229// }
1230
1231// m_mapModified.finalize();
1232// m_mapOriginal.finalize();
1233
1234
1235 // For the vertices, we eliminate as follows (where v is the valence):
1236 // - No elimination when v==1
1237 // - One on each side when v==2
1238 // - All but three when v>2
1239 for (size_t p=0; p!=m_bases.nBases(); p++)
1240 {
1241 for (index_t c=1; c<5; c++)
1242 {
1243 pcorner = patchCorner(p,c);
1245 _computeVertexMapper(pcorner);
1246 }
1247 }
1248 m_mapModified.finalize();
1249 m_mapOriginal.finalize();
1250
1251 GISMO_ASSERT((size_t)m_mapModified.freeSize()==m_size,"Size does not match predicted size, m_mapModified.freeSize()="<<m_mapModified.freeSize()<<"; m_size="<<m_size);
1252 m_matrix.resize( m_size, m_mapOriginal.freeSize() );
1253
1254 // gsDebugVar(m_mapModified.coupledSize());
1255 // gsDebugVar(m_mapModified.boundarySize());
1256
1257 }
1258
1259 // Same for DPatch and AlmostC1
1260 template<short_t d,class T>
1261 void gsDPatchBase<d,T>::_computeInterfaceMapper(boundaryInterface iface)
1262 {
1263 index_t sidx = _sideIndex( iface.second().patch,iface.second().side());
1264 if (m_sideCheck.at(sidx))
1265 return;
1266
1267 gsBasis<T> * basis;
1268 std::pair<index_t,bool> vdata1, vdata2;
1269 std::vector<index_t> patches(2);
1270 std::vector<patchSide> psides(2);
1271 gsVector<index_t> indices;
1272 std::vector<patchCorner> pcorners;
1273
1274 patches[0] = iface.first().patch;
1275 patches[1] = iface.second().patch;
1276 psides[0] = patchSide(iface.first().patch,iface.first().side()); // the interface on the first patch
1277 psides[1] = patchSide(iface.second().patch,iface.second().side()); // the interface on the second patch
1278
1279 for (index_t p = 0; p != 2; p++)
1280 {
1281 sidx = _sideIndex( patches[p] ,psides[p] );
1282 if (m_sideCheck.at(sidx))
1283 continue;
1284
1285 /*
1286 Eliminates the interior nodes on the interfaces
1287
1288 o o o @ X | X @ o o o X: eliminated DoFs by interface/boundary rule
1289 o o o @ X | X @ o o o @: modified DoFs by interface rule
1290 o o o @ X | X @ o o o o: preserved DoFs (interior)
1291 o o o x x | x x @ @ @ ?: Depends on the vertex (X and @ if not (interior vertex & valence = 4))
1292 o o o x x | x x X X X x: handled in vertex rule
1293 ----------|----------
1294 | x x X X X
1295 | x x @ @ @
1296 | o o o o o
1297 | o o o o o
1298 | o o o o o
1299
1300 */
1301
1302 basis = &m_bases.basis(patches[p]);
1303 indices = static_cast<gsVector<index_t>>( basis->boundary(psides[p]) );
1304
1305 patchSide(patches[p],psides[p]).getContainedCorners(d,pcorners);
1306 vdata1 = this->_vertexData(pcorners[0]);
1307 vdata2 = this->_vertexData(pcorners[1]);
1308
1309 // cast indices to an std::vector
1310 std::vector<index_t> allIndices(indices.data(), indices.data() + indices.rows() * indices.cols());
1311 // for(size_t i=0; i < allIndices.size(); i++)
1312 // std::cout << allIndices.at(i) << ' ';
1313
1314 std::vector<index_t> selectedIndices;
1315 // for both vertices of the side, add the indices at the vertex and one inside
1316 for (index_t c =0; c!=2; c++)
1317 {
1318 selectedIndices.push_back(_indexFromVert(0,pcorners[c],psides[p],0)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1319 selectedIndices.push_back(_indexFromVert(1,pcorners[c],psides[p],0)); // index from vertex pcorners[c] along side psides[0] with offset 0.
1320 }
1321
1322 std::sort(selectedIndices.begin(),selectedIndices.end());
1323 // for(size_t i=0; i < selectedIndices.size(); i++)
1324 // std::cout << selectedIndices.at(i) << ' ';
1325
1326 std::vector<index_t> result(allIndices.size());
1327 std::vector<index_t>::iterator it=std::set_difference (allIndices.begin(), allIndices.end(), selectedIndices.begin(), selectedIndices.end(), result.begin());
1328 result.resize(it-result.begin());
1329
1330 gsAsMatrix<index_t> indices(result,result.size(),1);
1331
1332 #pragma omp critical (side_interface)
1333 {
1334 m_mapModified.markBoundary(patches[p], indices);
1335 m_sideCheck.at(sidx) = true;
1336 }
1337 }
1338 }
1339
1340 template<short_t d,class T>
1341 void gsDPatchBase<d,T>::_computeBoundaryMapper(patchSide boundary)
1342 {
1343 index_t sidx = _sideIndex(boundary.patch,boundary.side());
1344 m_sideCheck.at(sidx) = true;
1345 }
1346
1347 template<short_t d,class T>
1348 void gsDPatchBase<d,T>::_computeVertexMapper(patchCorner pcorner)
1349 {
1350 index_t cidx = _vertIndex(pcorner.patch,pcorner.corner());
1351 if (m_vertCheck.at(cidx))
1352 return;
1353
1354 bool C0 = m_C0s[cidx];
1355 bool interior;
1356 index_t valence;
1357
1358 std::tie(valence,interior) = _vertexData(pcorner); // corner c
1359 if (!interior && valence==1) //valence = 1
1360 _computeMapperRegularCorner_v1(pcorner,valence);
1361 else if (!interior && valence==2 && C0)
1362 _computeMapperRegularBoundaryVertexNonSmooth_v2(pcorner,valence);
1363 else if (!interior && valence==2 && !C0)
1364 _computeMapperRegularBoundaryVertexSmooth_v2(pcorner,valence);
1365 else if (!interior && valence >2 && C0)
1366 _computeMapperIrregularBoundaryVertexNonSmooth_v(pcorner,valence);
1367 else if (!interior && valence >2 && !C0)
1368 _computeMapperIrregularBoundaryVertexSmooth_v(pcorner,valence);
1369 else if (interior)
1370 _computeMapperInteriorVertex_v(pcorner,valence);
1371 else
1372 GISMO_ERROR("Something went terribly wrong, interior="<<interior<<"; valence="<<valence);
1373
1374 // label vertex as processed
1375 m_vertCheck[ cidx ] = true;
1376 }
1377
1378 template<short_t d,class T>
1379 void gsDPatchBase<d,T>::_computeMapperRegularCorner_v1(patchCorner pcorner, index_t valence)
1380 {
1381 GISMO_UNUSED(pcorner);
1382 GISMO_UNUSED(valence);
1383
1384 // do nothing,
1385 }
1386
1387 // Same for DPatch and AlmostC1
1388 template<short_t d,class T>
1389 void gsDPatchBase<d,T>::_computeMapperRegularBoundaryVertexSmooth_v2(patchCorner pcorner, index_t valence)
1390 {
1391 GISMO_UNUSED(valence);
1392
1393 std::vector<patchSide> psides(2);
1394
1395 /*
1396 o o o @ X |i| X @ o o o x: eliminated DoFs by vertex rule
1397 o o o @ X |n| X @ o o o X: eliminated DoFs by interface/boundary rule
1398 o o o @ X |t| X @ o o o o: preserved DoFs (interior)
1399 o o o * x |e| x * o o o @: modified DoFs by interface rule
1400 o o o * x |r| x * o o o *: modified DoFs by vertex rule (unique DoFs)
1401 ----------------------- %: modified DoFs by vertex rule (matched DoFs)
1402 -boundary-| | -boundary-
1403 -----------------------
1404
1405 v = 4 (almost C1)
1406 o o o @ X |i| X @ o o o x: eliminated DoFs by vertex rule
1407 o o o @ X |n| X @ o o o X: eliminated DoFs by interface/boundary rule
1408 o o o @ X |t| X @ o o o o: preserved DoFs (interior)
1409 @ @ @ * x |e| x * @ @ @ @: modified DoFs by interface rule
1410 X X X x x |r| x x X X X *: modified DoFs by vertex rule (unique DoFs)
1411 -----------------------
1412 interface-| |-interface
1413 -----------------------
1414 X X X x x |i| x x X X X
1415 @ @ @ * x |n| x * @ @ @
1416 o o o @ X |t| X @ o o o
1417 o o o @ X |e| X @ o o o
1418 o o o @ X |r| X @ o o o
1419
1420 */
1421
1422 // we mark the nodes belonging to the interface
1423 pcorner.getContainingSides(d,psides);
1424 for (size_t p=0; p!=psides.size(); p++)
1425 {
1426 // the 0,k (k=0,1) DoF should be eliminated
1427 if (m_topology.isInterface(psides[p]))
1428 {
1429 for (index_t k=0; k!=2; k++)
1430 m_mapModified.eliminateDof(_indexFromVert(k,pcorner,psides[p],0),pcorner.patch);
1431 }
1432 }
1433
1434 // // we mark the nodes belonging to the interface
1435 // pcorner.getContainingSides(d,psides);
1436 // for (size_t p=0; p!=psides.size(); p++)
1437 // {
1438 // if (m_topology.isInterface(psides[p]))
1439 // {
1440 // // the 0,0 vertex should be eliminated
1441 // m_mapModified.eliminateDof(basis->functionAtCorner(pcorner),pcorner.patch);
1442 // }
1443 // }
1444 }
1445
1446 // Different for DPatch and AlmostC1
1447 template<short_t d,class T>
1448 void gsDPatchBase<d,T>::_computeMapperRegularBoundaryVertexNonSmooth_v2(patchCorner pcorner, index_t valence)
1449 {
1450 GISMO_UNUSED(valence);
1451
1452 std::vector<patchSide> psides(2);
1453
1454 /*
1455 o o o @ X |i| X @ o o o x: eliminated DoFs by vertex rule
1456 o o o @ X |n| X @ o o o X: eliminated DoFs by interface/boundary rule
1457 o o o @ X |t| X @ o o o o: preserved DoFs (interior)
1458 o o o * x |e| x * o o o @: modified DoFs by interface rule
1459 o o o o o |r| o o o o o *: modified DoFs by vertex rule (unique DoFs)
1460 ----------------------- %: modified DoFs by vertex rule (matched DoFs)
1461 -boundary-| | -boundary-
1462 -----------------------
1463 */
1464 // we mark the nodes belonging to the interface
1465 pcorner.getContainingSides(d,psides);
1466 for (size_t p=0; p!=psides.size(); p++)
1467 {
1468 if (m_topology.isInterface(psides[p]))
1469 {
1470 m_mapModified.eliminateDof(this->_indexFromVert(1,pcorner,psides[p],0),pcorner.patch);
1471 }
1472 }
1473 }
1474
1475 template<short_t d,class T>
1476 void gsDPatchBase<d,T>::_computeMapperIrregularBoundaryVertexSmooth_v(patchCorner pcorner, index_t valence)
1477 {
1478 // std::vector<std::pair<index_t,index_t>> indices0 = _getAllInterfaceIndices(pcorner,0,m_bases);
1479 // std::vector<std::pair<index_t,index_t>> indices1 = _getAllInterfaceIndices(pcorner,1,m_bases);
1480 // std::vector<patchCorner> pcorners;
1481 // m_topology.getCornerList(pcorner,pcorners);
1482 // for (std::vector<patchCorner>::iterator it=pcorners.begin(); it!=pcorners.end(); it++)
1483 // {
1484 // // mark the vertex as passed
1485 // m_vertCheck[ _vertIndex(it->patch, it->corner()) ] = true;
1486 // }
1487
1488 // // Eliminate the 1,0 and 0,1s
1489 // for (std::vector<std::pair<index_t,index_t>>::iterator it=indices1.begin(); it!=indices1.end(); it++)
1490 // m_mapModified.eliminateDof(it->second,it->first);
1491
1492 // _removeLowestIndices(indices0,3);
1493 // for (std::vector<std::pair<index_t,index_t>>::iterator it=indices0.begin(); it!=indices0.end(); it++)
1494 // m_mapModified.eliminateDof(it->second,it->first);
1495 gsWarn<<"C0 handling for boundary corners with valence >2 has not yet been implemented. Using the default approach\n";
1496 this->_computeMapperIrregularBoundaryVertexNonSmooth_v(pcorner,valence);
1497 }
1498
1499 template<short_t d,class T>
1500 void gsDPatchBase<d,T>::_computeMapperIrregularBoundaryVertexNonSmooth_v(patchCorner pcorner, index_t valence)
1501 {
1502 GISMO_UNUSED(valence);
1503
1504 // Get all the 0,1 or 1,0 DoFs on the interfaces
1505 // The 0,0 DoFs are kept but deleted later from the matrix
1506 std::vector<patchSide> psides(2);
1507 pcorner.getContainingSides(d,psides);
1508 for (size_t p=0; p!=psides.size(); p++)
1509 {
1510 // the 0,k (k=0,1) DoF should be eliminated
1511 if (m_topology.isInterface(psides[p]))
1512 m_mapModified.eliminateDof(_indexFromVert(1,pcorner,psides[p],0),pcorner.patch);
1513 }
1514 }
1515
1516 template<short_t d,class T>
1517 void gsDPatchBase<d,T>::_computeMapperInteriorVertex_v(patchCorner pcorner, index_t valence)
1518 {
1519 GISMO_UNUSED(valence);
1520
1521 std::vector<patchSide> psides(2);
1522 /*
1523 o o o @ X |i| X @ o o o x: eliminated DoFs by vertex rule
1524 o o o @ X |n| X @ o o o X: eliminated DoFs by interface/boundary rule
1525 o o o @ X |t| X @ o o o o: preserved DoFs (interior)
1526 @ @ @ * x |e| x * @ @ @ @: modified DoFs by interface rule
1527 X X X x x |r| x x X X X *: modified DoFs by vertex rule (unique DoFs)
1528 -----------------------
1529 interface-| |-interface
1530 -----------------------
1531 X X X x x |i| x x X X X
1532 @ @ @ * x |n| x * @ @ @
1533 o o o @ X |t| X @ o o o
1534 o o o @ X |e| X @ o o o
1535 o o o @ X |r| X @ o o o
1536
1537 */
1538 // we mark the nodes belonging to the interfaces (both sides bordering the vertex)
1539 pcorner.getContainingSides(d,psides);
1540 for (size_t p=0; p!=psides.size(); p++)
1541 {
1542 m_mapModified.eliminateDof(this->_indexFromVert(0,pcorner,psides[p],0),pcorner.patch);
1543 m_mapModified.eliminateDof(this->_indexFromVert(1,pcorner,psides[p],0),pcorner.patch);
1544 }
1545 }
1546
1547 template<short_t d,class T>
1548 void gsDPatchBase<d,T>::_handleRegularCorner(patchCorner pcorner)
1549 {
1550 std::vector<patchSide> psides;
1551 std::vector<index_t> indices(4);
1552 sparseEntry_t entries;
1553
1554 pcorner.getContainingSides(d,psides);
1555 indices[0] = _indexFromVert(0,pcorner,psides[0],0); // b00
1556 indices[1] = _indexFromVert(1,pcorner,psides[0],0); // b01
1557 indices[2] = _indexFromVert(1,pcorner,psides[1],0); // b10
1558 indices[3] = _indexFromVert(1,pcorner,psides[1],1); // b11
1559
1560 T weight = 1.0;
1561 index_t colIdx, rowIdx;
1562 for (std::vector<index_t>::iterator it = indices.begin(); it!=indices.end(); ++it)
1563 {
1564 rowIdx = m_mapModified.index(*it,pcorner.patch);
1565 colIdx = m_mapOriginal.index(*it,pcorner.patch);
1566 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1567 }
1568
1569 #pragma omp critical (handle_boundary_vertex_tt)
1570 {
1571 _pushAndCheck(entries);
1572 m_vertCheck[ _vertIndex(pcorner.patch, pcorner.corner()) ] = true;
1573 }
1574 // gsInfo<<"patch = "<<pcorner.patch<<", corner = "<<pcorner.corner()<<"\n";
1575 return;
1576 }
1577
1578 template<short_t d,class T>
1579 void gsDPatchBase<d,T>::_handleRegularBoundaryVertexSmooth(patchCorner pcorner, index_t valence)
1580 {
1581 GISMO_UNUSED(valence);
1582
1583 std::vector<patchSide> psides;
1584 std::vector<index_t> indices(3);
1585
1586 sparseEntry_t entries;
1587
1588 index_t colIdx, rowIdx;
1589 T weight;
1590 boundaryInterface iface;
1591
1592 // Get the sides joining at the corner.
1593 pcorner.getContainingSides(d,psides);
1594
1595 // 1. find the interface
1596 index_t iindex = m_topology.isInterface(psides[0]) ? 0 : 1;
1597
1598 GISMO_ENSURE(m_topology.getInterface(psides[iindex],iface),"Must be an interface");
1599
1600 // 2. collect indices
1601 // If we want C0 at this vertex, we only handle the row k=1.
1602 patchSide otherSide = iface.other(psides[iindex]);
1603 patchCorner otherCorner = iface.mapCorner(pcorner);
1604 for (index_t k = 0; k!=2; k++) // index of point over the interface
1605 {
1606 indices[0] = _indexFromVert(k,pcorner,psides[iindex],1); // bk1 on patch of iface
1607 indices[1] = _indexFromVert(k,pcorner,psides[iindex],0); // bk0 on patch of iface
1608 indices[2] = _indexFromVert(k,otherCorner,otherSide,0); // bk0 on other patch
1609
1610 rowIdx = m_mapModified.index(indices[0],pcorner.patch);
1611 colIdx = m_mapOriginal.index(indices[0],pcorner.patch);
1612 weight = 1.0;
1613 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1614 colIdx = m_mapOriginal.index(indices[1],psides[iindex].patch);
1615 weight = 0.5;
1616 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1617 colIdx = m_mapOriginal.index(indices[2],otherSide.patch);
1618 weight = 0.5;
1619 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1620 }
1621
1622 #pragma omp critical (handle_boundary_vertex_tt)
1623 {
1624 _pushAndCheck(entries);
1625 m_vertCheck[ _vertIndex(pcorner.patch, pcorner.corner()) ] = true;
1626 }
1627 }
1628
1629 template<short_t d,class T>
1630 void gsDPatchBase<d,T>::_handleRegularBoundaryVertexNonSmooth(patchCorner pcorner, index_t valence)
1631 {
1632 GISMO_UNUSED(valence);
1633
1634 std::vector<patchSide> psides;
1635 std::vector<index_t> indices(3);
1636
1637 sparseEntry_t entries;
1638
1639 index_t colIdx, rowIdx;
1640 T weight;
1641 boundaryInterface iface;
1642
1643 // Get the sides joining at the corner.
1644 pcorner.getContainingSides(d,psides);
1645
1646 // 1. find the interface
1647 index_t iindex = m_topology.isInterface(psides[0]) ? 0 : 1;
1648
1649 GISMO_ENSURE(m_topology.getInterface(psides[iindex],iface),"Must be an interface");
1650
1651 // 2. collect indices
1652 // If we want C0 at this vertex, we only handle the row k=1.
1653 patchSide otherSide = iface.other(psides[iindex]);
1654 patchCorner otherCorner = iface.mapCorner(pcorner);
1655 for (index_t k = 1; k!=2; k++) // index of point over the interface
1656 {
1657 indices[0] = _indexFromVert(k,pcorner,psides[iindex],1); // bk1 on patch of iface
1658 indices[1] = _indexFromVert(k,pcorner,psides[iindex],0); // bk0 on patch of iface
1659 indices[2] = _indexFromVert(k,otherCorner,otherSide,0); // bk0 on other patch
1660
1661 rowIdx = m_mapModified.index(indices[0],pcorner.patch);
1662 colIdx = m_mapOriginal.index(indices[0],pcorner.patch);
1663 weight = 1.0;
1664 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1665 colIdx = m_mapOriginal.index(indices[1],psides[iindex].patch);
1666 weight = 0.5;
1667 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1668 colIdx = m_mapOriginal.index(indices[2],otherSide.patch);
1669 weight = 0.5;
1670 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1671 }
1672
1673 #pragma omp critical (handle_boundary_vertex_tt)
1674 {
1675 _pushAndCheck(entries);
1676 m_vertCheck[ _vertIndex(pcorner.patch, pcorner.corner()) ] = true;
1677 }
1678 }
1679
1680 template<short_t d,class T>
1681 void gsDPatchBase<d,T>::_handleIrregularBoundaryVertexSmooth(patchCorner pcorner, index_t valence)
1682 {
1683 GISMO_UNUSED(valence);
1684
1685 gsWarn<<"C1 handling for boundary corners with valence >2 has not yet been implemented. Using the default approach\n";
1686 this->_handleIrregularBoundaryVertexNonSmooth(pcorner,valence);
1687 }
1688
1689 template<short_t d,class T>
1690 void gsDPatchBase<d,T>::_handleIrregularBoundaryVertexNonSmooth(patchCorner pcorner, index_t valence)
1691 {
1692 GISMO_UNUSED(valence);
1693
1694 std::vector<patchSide> psides;
1695 std::vector<patchCorner> corners;
1696 std::vector<index_t> indices;
1697 sparseEntry_t entries;
1698
1699 boundaryInterface iface;
1700 patchSide otherSide;
1701
1702 index_t colIdx, rowIdx;
1703 T weight;
1704
1705 pcorner.getContainingSides(d,psides);
1706
1707 std::vector<index_t> rowIndices, colIndices, patchIndices;
1708
1709 // pcorner is the current corner
1710 m_topology.getCornerList(pcorner,corners);
1711
1713 // Influence of 1,1 to itself
1714 index_t b11_p1 = _indexFromVert(1,pcorner,psides[0],1); // point 1,1 (does not matter which reference side is taken)
1715 rowIdx = m_mapModified.index(b11_p1,pcorner.patch);
1716 colIdx = m_mapOriginal.index(b11_p1,pcorner.patch);
1717
1718 weight = 1.;
1719 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1720
1722 index_t idx;
1723 for (index_t k = 0; k!=2; k++)
1724 {
1725 // Check if one of the adjacent interface is a boundary;
1726 // if so, add weight 1.0 to itself
1727 if (!m_topology.getInterface(psides[k],iface)) // check if the side is NOT an interface
1728 {
1729 idx = _indexFromVert(1,pcorner,psides[k],0);
1730 rowIdx = m_mapModified.index(idx,pcorner.patch); //1,0 corner (on the boundary)
1731 colIdx = m_mapOriginal.index(idx,pcorner.patch); //1,0 corner (on the boundary)
1732 weight = 1.0;
1733 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1734 }
1735 // else, add weight of 0.5 from the 0,1 or 1,0 vertices across the interface
1736 else
1737 {
1738 weight = 0.5;
1739 rowIdx = m_mapModified.index(b11_p1,pcorner.patch);
1740
1741 patchSide otherSide = iface.other(psides[k]);
1742 patchCorner otherCorner = iface.mapCorner(pcorner);
1743
1744 idx = _indexFromVert(1,pcorner,psides[k],0); // bk0 on patch of iface
1745 colIdx = m_mapOriginal.index(idx,pcorner.patch);
1746 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1747
1748 idx = _indexFromVert(1,otherCorner,otherSide,0); // bk0 on other patch
1749 colIdx = m_mapOriginal.index(idx,otherCorner.patch);
1750 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1751 }
1752 }
1753
1754 // Lastly, give the 0,0 a weight 1 to itself
1755 index_t b00_p1 = _indexFromVert(0,pcorner,psides[0],0); // point 0,0 (does not matter which reference side is taken)
1756 rowIdx = m_mapModified.index(b00_p1,pcorner.patch);
1757 colIdx = m_mapOriginal.index(b00_p1,pcorner.patch);
1758
1759 weight = 1.;
1760 entries.push_back(std::make_tuple(rowIdx,colIdx,weight));
1761
1762 #pragma omp critical (handle_boundary_vertex_ff)
1763 {
1764 _pushAndCheck(entries);
1765 m_vertCheck[ _vertIndex(pcorner.patch, pcorner.corner()) ] = true;
1766 }
1767 }
1768
1769} // namespace gismo
gismo::boxSide::parameter
bool parameter() const
Returns the parameter value (false=0=start, true=1=end) that corresponds to this side.
Definition gsBoundary.h:128
gismo::boxSide::direction
short_t direction() const
Returns the parametric direction orthogonal to this side.
Definition gsBoundary.h:113
gismo::gsAsMatrix
Creates a mapped object or data pointer to a matrix without copying data.
Definition gsAsMatrix.h:32
gismo::gsBasis
A basis represents a family of scalar basis functions defined over a common parameter domain.
Definition gsBasis.h:79
gismo::gsDPatchBase
Constructs the D-Patch, from which the transformation matrix can be called.
Definition gsDPatchBase.h:37
gismo::gsDPatchBase::_handleVertex
virtual void _handleVertex(patchCorner pcorner)
Handles a vertex in the global matrix.
Definition gsDPatchBase.hpp:970
gismo::gsDPatchBase::mapperInfo
virtual void mapperInfo() const
Returns for each basis function if it is free or eliminated.
Definition gsDPatchBase.hpp:134
gismo::gsDPatchBase::exportPatch
virtual gsGeometry< T > * exportPatch(index_t patch, bool computeCoefs=true)
Exports a single modified patch with index patch.
Definition gsDPatchBase.hpp:186
gismo::gsDPatchBase::defaultOptions
virtual void defaultOptions()
Sets the default options.
Definition gsDPatchBase.hpp:39
gismo::gsDPatchBase::_indexFromSides
virtual index_t _indexFromSides(index_t index1, const patchSide side1, index_t index2, const patchSide side2)
Computes the index of a basis function using sides as reference.
Definition gsDPatchBase.hpp:227
gismo::gsDPatchBase::_removeLowestCorners
virtual void _removeLowestCorners(std::vector< patchCorner > &pcorners, index_t n=3) const
From a list of patchCorners pcorners, remove all but the lowest n corners.
Definition gsDPatchBase.hpp:448
gismo::gsDPatchBase::vertexInfo
virtual void vertexInfo(patchCorner corner) const
Returns information about a vertex.
Definition gsDPatchBase.hpp:151
gismo::gsDPatchBase::sideInfo
virtual void sideInfo() const
Returns information for all the sides in the topology.
Definition gsDPatchBase.hpp:165
gismo::gsDPatchBase::_initialize
virtual void _initialize()
Initializes the class:-.
Definition gsDPatchBase.hpp:792
gismo::gsDPatchBase::_initTHB
virtual void _initTHB()
Initializes the matrix, the basis and the mappers.
Definition gsDPatchBase.hpp:823
gismo::gsDPatchBase::compute
virtual void compute()
Computes the construction.
Definition gsDPatchBase.hpp:26
gismo::gsDPatchBase::_initChecks
virtual void _initChecks()
Initializes the matrix, the basis and the mappers.
Definition gsDPatchBase.hpp:810
gismo::gsDPatchBase::_handleRegularCorner
virtual void _handleRegularCorner(patchCorner pcorner)
Handles a regular corner.
Definition gsDPatchBase.hpp:1548
gismo::gsDPatchBase::_handleInterior
virtual void _handleInterior()
Handles the interior in the global matrix.
Definition gsDPatchBase.hpp:1155
gismo::gsDPatchBase::_computeMapper
virtual void _computeMapper()
Calculates the mapper.
Definition gsDPatchBase.hpp:1209
gismo::gsDPatchBase::_initBasis
virtual void _initBasis()
Initializes the basis.
Definition gsDPatchBase.hpp:839
gismo::gsDPatchBase::_computeSmoothMatrix
virtual void _computeSmoothMatrix()
Calculates the smooth matrix.
Definition gsDPatchBase.hpp:1177
gismo::gsDPatchBase::_removeLowestIndices
virtual void _removeLowestIndices(std::vector< std::pair< index_t, index_t > > &indices, index_t n=3) const
From a list of tuples (patch,index), remove all but the lowest n tuples.
Definition gsDPatchBase.hpp:484
gismo::gsDPatchBase::_initMatrix
virtual void _initMatrix()
Initializes the matrix.
Definition gsDPatchBase.hpp:851
gismo::gsDPatchBase::_getLowestCorners
virtual void _getLowestCorners(std::vector< patchCorner > &pcorners, index_t n=3) const
From a list of patchCorners pcorners, get the lowest n corners.
Definition gsDPatchBase.hpp:435
gismo::gsDPatchBase::_indexFromVert
virtual index_t _indexFromVert(const index_t index, const patchCorner corner, const patchSide side, const index_t offsets=0) const
Computes the index of a basis function taking one corner and one side as reference.
Definition gsDPatchBase.hpp:269
gismo::gsDPatchBase::_getLowestIndices
virtual void _getLowestIndices(std::vector< std::pair< index_t, index_t > > &indices, index_t n=3) const
From a list of tuples (patch,index), get the lowest n tuples.
Definition gsDPatchBase.hpp:463
gismo::gsDPatchBase::cornerInfo
virtual void cornerInfo() const
Returns information for all the corners in the topology.
Definition gsDPatchBase.hpp:174
gismo::gsDPatchBase::_initMappers
virtual void _initMappers()
Initializes the matrix, the basis and the mappers.
Definition gsDPatchBase.hpp:844
gismo::gsDPatchBase::_performChecks
virtual void _performChecks(bool basis)
Performs checks on sides, vertices and bases.
Definition gsDPatchBase.hpp:875
gismo::gsDPatchBase::_whichHandled
virtual void _whichHandled()
Prints which DoFs have been handled and which have been eliminated.
Definition gsDPatchBase.hpp:857
gismo::gsDPatchBase::getSharpCorners
virtual std::vector< bool > getSharpCorners(T tol=1e-2) const
Checks if corners are sharp or not.
Definition gsDPatchBase.hpp:47
gismo::gsDPatchBase::_handleBoundary
virtual void _handleBoundary(patchSide side)
Handles a boundary in the global matrix.
Definition gsDPatchBase.hpp:1112
gismo::gsDPatchBase::_handleInterface
virtual void _handleInterface(boundaryInterface iface)
Handles an interface in the global matrix.
Definition gsDPatchBase.hpp:1002
gismo::gsDPatchBase::_resetChecks
virtual void _resetChecks(bool basis)
Resets checks on sides, vertices and bases.
Definition gsDPatchBase.hpp:890
gismo::gsDPatchBase::_vertexData
virtual const std::pair< index_t, bool > _vertexData(const patchCorner corner) const
Returns the valence and whether a corner is interior or boundary.
Definition gsDPatchBase.hpp:425
gismo::gsDofMapper
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition gsDofMapper.h:69
gismo::gsFunctionSet::basis
const gsBasis< T > & basis(const index_t k) const
Helper which casts and returns the k-th piece of this function set as a gsBasis.
Definition gsFunctionSet.hpp:33
gismo::gsGeometry
Abstract base class representing a geometry map.
Definition gsGeometry.h:93
gismo::gsMapData
the gsMapData is a cache of pre-computed function (map) values.
Definition gsFuncData.h:349
gismo::gsMapData::points
gsMatrix< T > points
input (parametric) points
Definition gsFuncData.h:372
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition gsMultiBasis.h:37
gismo::gsMultiBasis::basis
const gsBasis< T > & basis(const size_t i) const
Return the i-th basis block.
Definition gsMultiBasis.h:267
gismo::gsTHBSplineBasis
Truncated hierarchical B-spline basis.
Definition gsTHBSplineBasis.h:36
gismo::gsTensorBSplineBasis
A tensor product B-spline basis.
Definition gsTensorBSplineBasis.h:37
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gsAssembler.h
Provides generic assembler routines.
index_t
#define index_t
Definition gsConfig.h:32
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
gsWarn
#define gsWarn
Definition gsDebug.h:50
GISMO_UNUSED
#define GISMO_UNUSED(x)
Definition gsDebug.h:112
GISMO_ENSURE
#define GISMO_ENSURE(cond, message)
Definition gsDebug.h:102
gsInfo
#define gsInfo
Definition gsDebug.h:43
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsExprEvaluator.h
Generic expressions evaluator.
gsExprHelper.h
Generic expressions helper.
gsHTensorBasis.h
Provides definition of HTensorBasis abstract interface.
gsTHBSpline.h
Provides declaration of THBSplineBasis class.
gsTensorNurbsBasis.h
Provides declaration of TensorNurbsBasis abstract interface.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::NEED_OUTER_NORMAL
@ NEED_OUTER_NORMAL
Outward normal on the boundary.
Definition gsForwardDeclarations.h:86
gismo::give
S give(S &x)
Definition gsMemory.h:266
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::other
patchSide & other(const patchSide &ps)
Returns the second side if ps is the first side, otherwise returns the second side.
Definition gsBoundary.h:789
gismo::boundaryInterface::second
patchSide & second()
second, returns the second patchSide of this interface
Definition gsBoundary.h:782
gismo::patchCorner
Struct which represents a certain corner of a patch.
Definition gsBoundary.h:393
gismo::patchCorner::getContainingSides
void getContainingSides(short_t dim, std::vector< patchSide > &sides) const
returns a vector of patchSides that contain this corner
Definition gsBoundary.h:416
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
gismo::patchSide::getContainedCorners
void getContainedCorners(short_t dim, std::vector< patchCorner > &corners) const
returns the vector of the corners contained in the side
Definition gsBoundary.cpp:35