G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsC1SurfVertex.h
Go to the documentation of this file.
1
14#pragma once
15
16
17#include <gsUnstructuredSplines/src/gsG1AuxiliaryPatch.h>
18#include <gsUnstructuredSplines/src/gsC1SurfBasisVertex.h>
19
20
21namespace gismo {
22template<short_t d, class T>
23class gsC1SurfVertex
24{
26 typedef memory::shared_ptr<gsC1SurfVertex> Ptr;
27
29 typedef memory::unique_ptr<gsC1SurfVertex> uPtr;
30
31
32public:
34 ~gsC1SurfVertex() { }
35
36 gsC1SurfVertex(gsMultiPatch<T> & mp, const std::vector<index_t > patchesAroundVertex, const std::vector<index_t > vertexIndices) : m_mp(mp)
37 {
38 for(size_t i = 0; i < patchesAroundVertex.size(); i++)
39 {
40 auxGeom.push_back(gsG1AuxiliaryPatch<d,T>(mp.patch(patchesAroundVertex[i]), patchesAroundVertex[i]));
41 auxVertexIndices.push_back(vertexIndices[i]);
42 checkBoundary(mp, patchesAroundVertex[i], vertexIndices[i]);
43 }
44 sigma = 0.0;
45 }
46
47 gsMultiPatch<T> computeAuxTopology(){
48 gsMultiPatch<T> auxTop;
49 for(size_t i = 0; i < auxGeom.size(); i++)
50 {
51 auxTop.addPatch(auxGeom[i].getPatch());
52 }
53 auxTop.computeTopology();
54 return auxTop;
55 }
56
57 void reparametrizeG1Vertex()
58 {
59 for(size_t i = 0; i < auxGeom.size(); i++)
60 {
61 checkOrientation(i); // Check if the orientation is correct. If not, modifies vertex and edge vectors
62
63 switch (auxVertexIndices[i])
64 {
65 case 1:
66 break;
67 case 4:
68 auxGeom[i].rotateParamAntiClockTwice();
69 break;
70 case 2:
71 auxGeom[i].rotateParamAntiClock();
72 break;
73 case 3:
74 auxGeom[i].rotateParamClock();
75 break;
76 }
77 }
78 }
79
80 index_t kindOfVertex()
81 {
82 if(auxGeom.size() == 1)
83 return -1; // Boundary vertex
84
85 gsMultiPatch<T> top(computeAuxTopology());
86 index_t nInt = top.interfaces().size();
87 if((index_t)auxGeom.size() == nInt)
88 return 0; // Internal vertex
89 else
90 return 1; // Interface-Boundary vertex
91 }
92
93 void checkOrientation(index_t i)
94 {
95 if (auxGeom[i].getPatch().orientation() == -1)
96 {
97 auxGeom[i].swapAxis();
98 this->swapBdy(i); //Swap boundary edge bool-value
99
100 // Swap vertices index after swapping axis
101 if(auxVertexIndices[i] == 2)
102 auxVertexIndices[i] = 3;
103 else
104 if(auxVertexIndices[i] == 3)
105 auxVertexIndices[i] = 2;
106 }
107 }
108
109 void computeSigma()
110 {
111 T p = 0;
112 T h_geo = 0;
113 for(size_t i = 0; i < auxGeom.size(); i++)
114 {
115 gsTensorBSplineBasis<d, T> & bsp_temp = dynamic_cast<gsTensorBSplineBasis<d, T> & >(auxGeom[i].getPatch().basis());
116 T p_temp = bsp_temp.maxDegree();
117 p = (p < p_temp ? p_temp : p);
118
119 for(index_t j = 0; j < auxGeom[i].getPatch().parDim(); j++)
120 {
121 T h_geo_temp = bsp_temp.component(j).knots().at(p + 2);
122 h_geo = (h_geo < h_geo_temp ? h_geo_temp : h_geo);
123 }
124 }
125 T val = auxGeom.size();
126
127 gsMatrix<T> zero;
128 zero.setZero(2,1);
129 for (index_t i = 0; i < val; i++)
130 sigma += auxGeom[i].getPatch().deriv(zero).template lpNorm<gsEigen::Infinity>();
131 sigma *= h_geo/(val*p);
132 sigma = 1 / sigma;
133 }
134
135 void checkBoundary(gsMultiPatch<T> & mpTmp, index_t patchInd, index_t sideInd)
136 {
137 std::vector<bool> tmp;
138 switch (sideInd)
139 {
140 case 1: tmp.push_back(mpTmp.isBoundary(patchInd,3));
141 tmp.push_back(mpTmp.isBoundary(patchInd,1));
142 break;
143 case 2: tmp.push_back(mpTmp.isBoundary(patchInd, 2));
144 tmp.push_back(mpTmp.isBoundary(patchInd, 3));
145 break;
146 case 3: tmp.push_back(mpTmp.isBoundary(patchInd, 1));
147 tmp.push_back(mpTmp.isBoundary(patchInd, 4));
148 break;
149 case 4: tmp.push_back(mpTmp.isBoundary(patchInd, 4));
150 tmp.push_back(mpTmp.isBoundary(patchInd, 2));
151 break;
152 default:
153 break;
154 }
155 isBdy.push_back(tmp);
156 }
157
158 void swapBdy(index_t i)
159 {
160 bool tmp = isBdy[i][0];
161 isBdy[i][0] = isBdy[i][1];
162 isBdy[i][1] = tmp;
163 }
164
165 gsMatrix<T> computeBigSystemMatrix( index_t np)
166 {
167 gsMultiBasis<T> bas(auxGeom[np].getPatch());
168 gsTensorBSplineBasis<d, T> & temp_L = dynamic_cast<gsTensorBSplineBasis<d, T> &>(bas.basis(0));
169 index_t dimU = temp_L.size(0);
170 index_t dimV = temp_L.size(1);
171
172 gsMatrix<T> BigMatrix;
173 BigMatrix.setZero( 2 * (dimU + dimV - 2),auxGeom[np].getG1Basis().nPatches());
174
175 for(index_t bf = 0; bf < auxGeom[np].getG1Basis().nPatches(); bf++)
176 {
177 for (index_t i = 0; i < 2 * dimU; i++)
178 {
179 if (auxGeom[np].getG1BasisCoefs(bf).at(i) * auxGeom[np].getG1BasisCoefs(bf).at(i) > m_zero*m_zero)
180 BigMatrix(i, bf) = auxGeom[np].getG1BasisCoefs(bf).at(i);
181 }
182
183 for (index_t i = 1; i < dimV - 1; i++)
184 {
185 for(index_t j = i; j < i + 2; j++)
186 {
187 if (auxGeom[np].getG1BasisCoefs(bf).at((i + 1) * dimU + j - i) * auxGeom[np].getG1BasisCoefs(bf).at((i + 1) * dimU + j - i) > m_zero*m_zero)
188 BigMatrix(i + j + (2 * dimU ) - 2, bf) = auxGeom[np].getG1BasisCoefs(bf).at((i + 1) * dimU + j - i);
189 }
190 }
191 }
192 return BigMatrix;
193 }
194
195 gsMatrix<T> computeSmallSystemMatrix( index_t np)
196 {
197 gsMultiBasis<T> bas(auxGeom[np].getPatch());
198 gsTensorBSplineBasis<d, T> & temp_L = dynamic_cast<gsTensorBSplineBasis<d, T> &>(bas.basis(0));
199 index_t dimU = temp_L.size(0);
200 index_t dimV = temp_L.size(1);
201
202 gsMatrix<T> SmallMatrix;
203 SmallMatrix.setZero((dimU + dimV - 1),auxGeom[np].getG1Basis().nPatches());
204
205 for(index_t bf = 0; bf < auxGeom[np].getG1Basis().nPatches(); bf++)
206 {
207 for (index_t i = 0; i < dimU; i++)
208 {
209 if (auxGeom[np].getG1BasisCoefs(bf).at(i) * auxGeom[np].getG1BasisCoefs(bf).at(i) > m_zero*m_zero)
210 SmallMatrix(i, bf) = auxGeom[np].getG1BasisCoefs(bf).at(i);
211 }
212
213 for (index_t i = 1; i < dimV; i++)
214 {
215 if (auxGeom[np].getG1BasisCoefs(bf).at(i * dimU) * auxGeom[np].getG1BasisCoefs(bf).at(i * dimU) > m_zero*m_zero)
216 SmallMatrix(i + dimU -1, bf) = auxGeom[np].getG1BasisCoefs(bf).at(i * dimU);
217 }
218 }
219 return SmallMatrix;
220 }
221
222 gsMatrix<T> leftBoundaryBigSystem(index_t np)
223 {
224 gsMultiBasis<T> bas(auxGeom[np].getPatch());
225 gsTensorBSplineBasis<d, T> & temp_L = dynamic_cast<gsTensorBSplineBasis<d, T> &>(bas.basis(0));
226 index_t dimU = temp_L.size(0);
227 index_t dimV = temp_L.size(1);
228
229 gsMatrix<T> BigMatrix;
230 BigMatrix.setZero( 2 * dimV,6);
231
232 for(index_t bf = 0; bf < 6; bf++)
233 {
234 for (index_t i = 0; i < dimV ; i++)
235 {
236 for(index_t j = i; j < i + 2; j++)
237 {
238 if (auxGeom[np].getG1BasisCoefs(bf).at(i * dimU + j - i) * auxGeom[np].getG1BasisCoefs(bf).at(i * dimU + j - i) > m_zero*m_zero)
239 BigMatrix(i + j, bf) = auxGeom[np].getG1BasisCoefs(bf).at(i * dimU + j - i);
240 else
241 auxGeom[np].getG1BasisCoefs(bf).at(i) *= 0;
242 }
243 }
244 }
245 return BigMatrix;
246 }
247
248 gsMatrix<T> rightBoundaryBigSystem( index_t np)
249 {
250 gsMultiBasis<T> bas(auxGeom[np].getPatch());
251 gsTensorBSplineBasis<d, T> & temp_L = dynamic_cast<gsTensorBSplineBasis<d, T> &>(bas.basis(0));
252 index_t dimU = temp_L.size(0);
253
254 gsMatrix<T> BigMatrix;
255 BigMatrix.setZero( 2 * dimU ,6);
256
257 for(index_t bf = 0; bf < 6; bf++)
258 {
259 for (index_t i = 0; i < 2 * dimU; i++)
260 {
261 if (auxGeom[np].getG1BasisCoefs(bf).at(i) * auxGeom[np].getG1BasisCoefs(bf).at(i) > m_zero*m_zero)
262 BigMatrix(i, bf) = auxGeom[np].getG1BasisCoefs(bf).at(i);
263 else
264 auxGeom[np].getG1BasisCoefs(bf).at(i) *= 0;
265 }
266 }
267 return BigMatrix;
268 }
269
270 gsMatrix<T> leftBoundarySmallSystem( index_t np)
271 {
272 gsMultiBasis<T> bas(auxGeom[np].getPatch());
273 gsTensorBSplineBasis<d, T> & temp_L = dynamic_cast<gsTensorBSplineBasis<d, T> &>(bas.basis(0));
274 index_t dimU = temp_L.size(0);
275 index_t dimV = temp_L.size(1);
276
277 gsMatrix<T> SmallMatrix;
278 SmallMatrix.setZero(dimV,6);
279
280 for(index_t bf = 0; bf < 6; bf++)
281 {
282 for (index_t i = 0; i < dimV; i++)
283 {
284 if (auxGeom[np].getG1BasisCoefs(bf).at(i * dimU) * auxGeom[np].getG1BasisCoefs(bf).at(i * dimU) > m_zero*m_zero)
285 SmallMatrix(i, bf) = auxGeom[np].getG1BasisCoefs(bf).at(i * dimU);
286 else
287 auxGeom[np].getG1BasisCoefs(bf).at(i * dimU) *= 0;
288 }
289 }
290 return SmallMatrix;
291 }
292
293 gsMatrix<T> rightBoundarySmallSystem( index_t np)
294 {
295 gsMultiBasis<T> bas(auxGeom[np].getPatch());
296 gsTensorBSplineBasis<d, T> & temp_L = dynamic_cast<gsTensorBSplineBasis<d, T> &>(bas.basis(0));
297 index_t dimU = temp_L.size(0);
298
299 gsMatrix<T> SmallMatrix;
300 SmallMatrix.setZero( dimU, 6);
301
302 for(index_t bf = 0; bf < 6; bf++)
303 {
304 for (index_t i = 0; i < dimU; i++)
305 {
306 if (auxGeom[np].getG1BasisCoefs(bf).at(i ) * auxGeom[np].getG1BasisCoefs(bf).at(i ) > m_zero*m_zero)
307 SmallMatrix(i, bf) = auxGeom[np].getG1BasisCoefs(bf).at(i);
308 else
309 auxGeom[np].getG1BasisCoefs(bf).at(i) *= 0;
310 }
311 }
312 return SmallMatrix;
313 }
314
315 gsMatrix<T> bigInternalBoundaryPatchSystem( index_t np)
316 {
317 gsMultiBasis<T> bas(auxGeom[np].getPatch());
318 gsTensorBSplineBasis<d, T> & temp_L = dynamic_cast<gsTensorBSplineBasis<d, T> &>(bas.basis(0));
319 index_t dimU = temp_L.size(0);
320
321 gsMatrix<T> Matrix;
322 Matrix.setZero( 3 ,6);
323
324 for(index_t bf = 0; bf < 6; bf++)
325 {
326 Matrix(0, bf) = auxGeom[np].getG1BasisCoefs(bf).at(0);
327 Matrix(1, bf) = auxGeom[np].getG1BasisCoefs(bf).at(1);
328 Matrix(2, bf) = auxGeom[np].getG1BasisCoefs(bf).at(dimU);
329
330 }
331 return Matrix;
332 }
333
334 gsMatrix<T> smallInternalBoundaryPatchSystem( index_t np)
335 {
336 gsMatrix<T> Matrix;
337 Matrix.setZero( 1 ,6);
338
339 for(index_t bf = 0; bf < 6; bf++)
340 {
341 Matrix(0, bf) = auxGeom[np].getG1BasisCoefs(bf).at(0);
342 }
343 return Matrix;
344 }
345
346 std::pair<gsMatrix<T>, gsMatrix<T>> createSinglePatchSystem(index_t np)
347 {
348 if(isBdy[np][1] == 1)
349 return std::make_pair(leftBoundaryBigSystem(np), leftBoundarySmallSystem(np));
350 else
351 {
352 if (isBdy[np][0] == 1)
353 return std::make_pair(rightBoundaryBigSystem(np), rightBoundarySmallSystem(np));
354 else
355 return std::make_pair(bigInternalBoundaryPatchSystem(np), smallInternalBoundaryPatchSystem(np));
356 }
357 }
358
359 void checkValues(gsMatrix<T> & mat)
360 {
361 for(index_t bk = 0; bk < mat.cols(); bk++ )
362 {
363 for(index_t r=0; r < mat.rows(); r++)
364 {
365 if( abs( mat(r, bk) * mat(r, bk) ) < (m_zero*m_zero) )
366 mat(r, bk) = 0;
367 }
368 }
369 }
370
371 void addVertexBasis(gsMatrix<T> & basisV)
372 {
373 gsMatrix<T> vertBas;
374 vertBas.setIdentity(auxGeom[0].getG1Basis().nPatches(), auxGeom[0].getG1Basis().nPatches());
375 index_t count = 0;
376 index_t numBF = auxGeom[0].getG1Basis().nPatches();
377 while (basisV.cols() < numBF)
378 {
379 basisV.conservativeResize(basisV.rows(), basisV.cols() + 1);
380 basisV.col(basisV.cols() - 1) = vertBas.col(count);
381
382 gsEigen::FullPivLU<gsMatrix<T>> ker(basisV);
383 ker.setThreshold(1e-10);
384 if (ker.dimensionOfKernel() != 0)
385 {
386 basisV = basisV.block(0, 0, basisV.rows(), basisV.cols() - 1);
387 }
388 count++;
389 }
390 }
391
392 void addSmallKerBasis(gsMatrix<T> & basisV, gsMatrix<T> & smallK, index_t smallKDim)
393 {
394 for(index_t i=0; i < smallKDim; i++)
395 {
396 basisV.conservativeResize(basisV.rows(), basisV.cols() + 1);
397 basisV.col(basisV.cols()-1) = smallK.col(i);
398
399 gsEigen::FullPivLU<gsMatrix<T>> ker(basisV);
400 ker.setThreshold(1e-10);
401 if(ker.dimensionOfKernel() != 0)
402 {
403 basisV = basisV.block(0, 0, basisV.rows(), basisV.cols()-1);
404 }
405 }
406
407 }
408
409 std::pair<gsMatrix<T>, std::vector<index_t>> selectVertexBoundaryBasisFunction(gsMatrix<T> bigKernel, index_t bigKerDim, gsMatrix<T> smallKernel, index_t smallKerDim)
410 {
411 gsMatrix<T> basisVect;
412 std::vector<index_t> numberPerType;
413
414 numberPerType.push_back(bigKerDim); // Number of basis which has to be moved to the internal
415 numberPerType.push_back(smallKerDim - bigKerDim); // Number of basis which are boundary function of FIRST TYPE
416 numberPerType.push_back(bigKernel.cols() - smallKerDim); // Number of basis which are boundary function of SECOND TYPE
417
418 if(bigKerDim != 0)
419 {
420 checkValues(bigKernel);
421 checkValues(smallKernel);
422
423 basisVect = bigKernel;
424 addSmallKerBasis(basisVect, smallKernel, smallKerDim);
425 addVertexBasis(basisVect);
426 }
427 else
428 {
429 if (smallKerDim != 0)
430 {
431 checkValues(smallKernel);
432
433 basisVect = smallKernel;
434 addVertexBasis(basisVect);
435 }
436 else
437 {
438 gsMatrix<T> vertBas;
439 vertBas.setIdentity(auxGeom[0].getG1Basis().nPatches(), auxGeom[0].getG1Basis().nPatches());
440 basisVect = vertBas;
441 }
442 }
443 return std::make_pair(basisVect, numberPerType);
444 }
445
446 gsMatrix<T> selectGD(index_t i)
447 {
448 gsMatrix<T> coefs(4, 2);
449
450 if( kindOfVertex() == 1 ) // If the boundary it´s along u and along v there is an interface (Right Patch) or viceversa
451 {
452 gsMultiPatch<T> tmp(this->computeAuxTopology());
453 for(auto iter : tmp.interfaces())
454 {
455 if( i == iter.first().patch || i == iter.second().patch )
456 {
457 gsMultiPatch<T> aux;
458 if( iter.first().index() == 1 )
459 {
460 aux.addPatch(tmp.patch(iter.first().patch));
461 aux.addPatch(tmp.patch(iter.second().patch));
462 }
463 else
464 {
465 aux.addPatch(tmp.patch(iter.second().patch));
466 aux.addPatch(tmp.patch(iter.first().patch));
467 }
468
469 aux.computeTopology();
470 gsMultiBasis<T> auxB(aux);
471 gsC1SurfGluingData<T> ret(aux, auxB);
472 gsMatrix<T> sol = ret.getSol();
473 gsMatrix<T> solBeta = ret.getSolBeta();
474
475 if( (isBdy[i][0] == 0) && (isBdy[i][1] == 0))
476 {
477 if ( (i == iter.first().patch && iter.first().index() == 3)
478 || (i == iter.second().patch && iter.second().index() == 3) )
479 {
480 coefs(0, 0) = sol(2, 0);
481 coefs(1, 0) = sol(3, 0);
482 coefs(2, 0) = solBeta(2, 0);
483 coefs(3, 0) = solBeta(3, 0);
484 }
485 else
486 if ( (i == iter.first().patch && iter.first().index() == 1)
487 || (i == iter.second().patch && iter.second().index() == 1))
488 {
489 coefs(0, 1) = sol(0, 0);
490 coefs(1, 1) = sol(1, 0);
491 coefs(2, 1) = solBeta(0, 0);
492 coefs(3, 1) = solBeta(1, 0);
493 }
494 }
495 else
496 {
497 if ((i == iter.first().patch && iter.first().index() == 3)
498 || (i == iter.second().patch && iter.second().index() == 3))
499 {
500 coefs(0, 0) = sol(2, 0);
501 coefs(0, 1) = 1;
502 coefs(1, 0) = sol(3, 0);
503 coefs(1, 1) = 1;
504 coefs(2, 0) = solBeta(2, 0);
505 coefs(2, 1) = 0;
506 coefs(3, 0) = solBeta(3, 0);
507 coefs(3, 1) = 0;
508 }
509 else if ((i == iter.first().patch && iter.first().index() == 1)
510 || (i == iter.second().patch && iter.second().index() == 1))
511 {
512 coefs(0, 0) = 1;
513 coefs(0, 1) = sol(0, 0);
514 coefs(1, 0) = 1;
515 coefs(1, 1) = sol(1, 0);
516 coefs(2, 0) = 0;
517 coefs(2, 1) = solBeta(0, 0);
518 coefs(3, 0) = 0;
519 coefs(3, 1) = solBeta(1, 0);
520 }
521 }
522 }
523 }
524 }
525 else
526 if( kindOfVertex() == -1 ) // Single patch corner
527 {
528 coefs.setZero();
529
530 coefs(0, 0) = 1;
531 coefs(0, 1) = 1;
532 coefs(1, 0) = 1;
533 coefs(1, 1) = 1;
534 }
535 else
536 if( kindOfVertex() == 0 ) // Internal vertex -> Two interfaces
537 {
538 gsMultiPatch<T> tmp(this->computeAuxTopology());
539 for(auto iter : tmp.interfaces())
540 {
541 if( (i == iter.first().patch) || (i == iter.second().patch) )
542 {
543 gsMultiPatch<T> aux;
544 if( iter.first().index() == 1 )
545 {
546 aux.addPatch(tmp.patch(iter.first().patch));
547 aux.addPatch(tmp.patch(iter.second().patch));
548 }
549 else
550 {
551 aux.addPatch(tmp.patch(iter.second().patch));
552 aux.addPatch(tmp.patch(iter.first().patch));
553 }
554
555 aux.computeTopology();
556 gsMultiBasis<T> auxB(aux);
557 gsC1SurfGluingData<T> ret(aux, auxB);
558 gsMatrix<T> sol = ret.getSol();
559 gsMatrix<T> solBeta = ret.getSolBeta();
560
561 if ( (i == iter.first().patch && iter.first().index() == 3)
562 || (i == iter.second().patch && iter.second().index() == 3) )
563 {
564 coefs(0, 0) = sol(2, 0);
565 coefs(1, 0) = sol(3, 0);
566 coefs(2, 0) = solBeta(2, 0);
567 coefs(3, 0) = solBeta(3, 0);
568 }
569 else
570 if ( (i == iter.first().patch && iter.first().index() == 1)
571 || (i == iter.second().patch && iter.second().index() == 1))
572 {
573 coefs(0, 1) = sol(0, 0);
574 coefs(1, 1) = sol(1, 0);
575 coefs(2, 1) = solBeta(0, 0);
576 coefs(3, 1) = solBeta(1, 0);
577 }
578 }
579 }
580 }
581 return coefs;
582 }
583
584 void computeG1InternalVertexBasis()
585 {
586
587 m_zero = 1e-10;
588
589 this->reparametrizeG1Vertex();
590 this->computeSigma();
591
592 std::vector<gsMultiPatch<T>> g1BasisVector;
593 std::pair<gsMatrix<T>, std::vector<index_t>> vertexBoundaryBasis;
594
595 gsMatrix<T> Phi(6, 6);
596 Phi.setIdentity();
597
598 Phi.col(1) *= sigma;
599 Phi.col(2) *= sigma;
600 Phi.col(3) *= sigma * sigma;
601 Phi.col(4) *= sigma * sigma;
602 Phi.col(5) *= sigma * sigma;
603
604 gsMultiPatch<T> rotPatch;
605
606 if (auxGeom[0].getPatch().parDim() + 1 == auxGeom[0].getPatch().targetDim())
607 {
608 gsMatrix<T> zero;
609 zero.setZero(2, 1);
610 gsMatrix<T> Jk = auxGeom[0].getPatch().jacobian(zero);
611 gsMatrix<T> G = Jk.transpose() * Jk; // Symmetric
612 gsMatrix<T> G_inv = G.cramerInverse(); // Symmetric
613
614
615 gsMatrix<T> geoMapDeriv1 = auxGeom[0].getPatch()
616 .deriv(zero); // First derivative of the geometric mapping with respect to the parameter coordinates
617 gsMatrix<T> geoMapDeriv2 = auxGeom[0].getPatch()
618 .deriv2(zero); // Second derivative of the geometric mapping with respect to the parameter coordinates
619
620 //Computing the normal vector to the tangent plane along the boundary curve
621// gsVector<T,3> t1 = Jk.col(0);
622// gsVector<T,3> t2 = Jk.col(1);
623//
624// gsVector<T,3> n = t1.cross(t2);
625//
626// gsVector<T> normal = n.normalized();
627// n = n.normalized();
628// gsVector<T,3> z(0, 0, 1);
629//
630// gsVector<T,3> rotVec = n.cross(z);
631// rotVec = rotVec.normalized();
632//
633// T cos_t = n.dot(z) / (n.norm() * z.norm());
634// T sin_t = (n.cross(z)).norm() / (n.norm() * z.norm());
635
636 gsVector<T> n(3);
637 n.setZero();
638 n(0) = Jk(1,0)*Jk(2,1)-Jk(2,0)*Jk(1,1);
639 n(1) = Jk(2,0)*Jk(0,1)-Jk(0,0)*Jk(2,1);
640 n(2) = Jk(0,0)*Jk(1,1)-Jk(1,0)*Jk(0,1);
641
642 gsVector<T> z(3);
643 z.setZero();
644 z(2) = 1.0;
645
646 gsVector<T> rotVec(3);
647 rotVec.setZero(3);
648 rotVec(0) = n(1,0)*z(2,0)-n(2,0)*z(1,0);
649 rotVec(1) = n(2,0)*z(0,0)-n(0,0)*z(2,0);
650 rotVec(2) = n(0,0)*z(1,0)-n(1,0)*z(0,0);
651
652 T cos_t = (n.dot(z))/ (n.norm() * z.norm());
653 T sin_t = rotVec.norm() / (n.norm() * z.norm());
654
655// Rotation matrix
656 gsMatrix<T> R(3, 3);
657 R.setZero();
658// Row 0
659 R(0, 0) = cos_t + rotVec.x() * rotVec.x() * (1 - cos_t);
660 R(0, 1) = rotVec.x() * rotVec.y() * (1 - cos_t) - rotVec.z() * sin_t;
661 R(0, 2) = rotVec.x() * rotVec.z() * (1 - cos_t) + rotVec.y() * sin_t;
662// Row 1
663 R(1, 0) = rotVec.x() * rotVec.y() * (1 - cos_t) + rotVec.z() * sin_t;
664 R(1, 1) = cos_t + rotVec.y() * rotVec.y() * (1 - cos_t);
665 R(1, 2) = rotVec.y() * rotVec.z() * (1 - cos_t) - rotVec.x() * sin_t;
666// Row 2
667 R(2, 0) = rotVec.x() * rotVec.z() * (1 - cos_t) - rotVec.y() * sin_t;
668 R(2, 1) = rotVec.y() * rotVec.z() * (1 - cos_t) + rotVec.x() * sin_t;
669 R(2, 2) = cos_t + rotVec.z() * rotVec.z() * (1 - cos_t);
670
671 for (size_t np = 0; np < auxGeom.size(); np++)
672 {
673 gsMatrix<T> coeffPatch = auxGeom[np].getPatch().coefs();
674
675 for (index_t i = 0; i < coeffPatch.rows(); i++)
676 {
677 coeffPatch.row(i) =
678 (coeffPatch.row(i) - coeffPatch.row(0)) * R.transpose() + coeffPatch.row(0);
679 }
680
681 rotPatch.addPatch(auxGeom[np].getPatch());
682 rotPatch.patch(np).setCoefs(coeffPatch);
683 }
684
685 Phi.resize(13, 6);
686 Phi.setZero();
687 Phi(0, 0) = 1;
688 Phi(1, 1) = sigma;
689 Phi(2, 2) = sigma;
690 Phi(4, 3) = sigma * sigma;
691 Phi(5, 4) = sigma * sigma;
692 Phi(7, 4) = sigma * sigma;
693 Phi(8, 5) = sigma * sigma;
694
695 for (size_t i = 0; i < auxGeom.size(); i++)
696 {
697 gsMatrix<T> gdCoefs(selectGD(i));
698 gsMultiPatch<T> g1Basis;
699
700 gsC1SurfBasisVertex<T> g1BasisVertex_0(rotPatch.patch(i), rotPatch.patch(i).basis(), isBdy[i], Phi, gdCoefs);
701
702 g1BasisVertex_0.setG1BasisVertex(g1Basis);
703
704 g1BasisVector.push_back(g1Basis);
705 auxGeom[i].setG1Basis(g1Basis);
706 }
707 }
708 else
709 {
710 for (size_t i = 0; i < auxGeom.size(); i++)
711 {
712 gsMatrix<T> gdCoefs(selectGD(i));
713 gsMultiPatch<T> g1Basis;
714
715 gsC1SurfBasisVertex<T> g1BasisVertex_0(auxGeom[i].getPatch(), auxGeom[i].getPatch().basis(), isBdy[i], Phi, gdCoefs);
716
717 g1BasisVertex_0.setG1BasisVertex(g1Basis);
718
719 g1BasisVector.push_back(g1Basis);
720 auxGeom[i].setG1Basis(g1Basis);
721 }
722 }
723
724 // OLD
725// if (this->kindOfVertex() == 1) // Interface-Boundary vertex
726// {
727// gsMatrix<T> bigMatrix(0,0);
728// gsMatrix<T> smallMatrix(0,0);
729// for (index_t i = 0; i < auxGeom.size(); i++)
730// {
731// std::pair<gsMatrix<T>, gsMatrix<T>> tmp;
732// tmp = createSinglePatchSystem(i);
733// index_t row_bigMatrix = bigMatrix.rows();
734// index_t row_smallMatrix = smallMatrix.rows();
735//
736// bigMatrix.conservativeResize(bigMatrix.rows() + tmp.first.rows(), 6);
737// smallMatrix.conservativeResize(smallMatrix.rows() + tmp.second.rows(), 6);
738//
739// bigMatrix.block(row_bigMatrix, 0, tmp.first.rows(), 6) = tmp.first;
740// smallMatrix.block(row_smallMatrix, 0, tmp.second.rows(), 6) = tmp.second;
741// }
742//
743// gsEigen::FullPivLU<gsMatrix<T>> BigLU(bigMatrix);
744// gsEigen::FullPivLU<gsMatrix<T>> SmallLU(smallMatrix);
745// SmallLU.setThreshold(1e-10);
746// BigLU.setThreshold(1e-10);
747//
752//
753// vertexBoundaryBasis = selectVertexBoundaryBasisFunction(BigLU.kernel(), BigLU.dimensionOfKernel(), SmallLU.kernel(), SmallLU.dimensionOfKernel());
754//
755// }
756// else if(this->kindOfVertex() == -1) // Boundary vertex
757// {
758// gsEigen::FullPivLU<gsMatrix<T>> BigLU(computeBigSystemMatrix(0));
759// gsEigen::FullPivLU<gsMatrix<T>> SmallLU(computeSmallSystemMatrix(0));
760// SmallLU.setThreshold(1e-10);
761// BigLU.setThreshold(1e-10);
762//
767//
768// vertexBoundaryBasis = selectVertexBoundaryBasisFunction(BigLU.kernel(), BigLU.dimensionOfKernel(), SmallLU.kernel(), SmallLU.dimensionOfKernel());
769//
770// }
771//
772// if (this->kindOfVertex() != 0)
773// for (index_t i = 0; i < auxGeom.size(); i++)
774// {
775// gsMultiPatch<T> temp_mp_g1 = g1BasisVector[i];
776// for (index_t bf = 0; bf < temp_mp_g1.nPatches(); bf++)
777// {
779// gsMatrix<T> coef_bf;
780// coef_bf.setZero(temp_mp_g1.patch(bf).coefs().dim().first,1);
781// for (index_t lambda = 0; lambda < temp_mp_g1.nPatches(); lambda++)
782// coef_bf += temp_mp_g1.patch(lambda).coefs() * vertexBoundaryBasis.first(lambda,bf);
783//
784// g1BasisVector[i].patch(bf).setCoefs(coef_bf);
785// }
786// auxGeom[i].parametrizeBasisBack(g1BasisVector[i]);
787// }
788// else
789// for (index_t i = 0; i < auxGeom.size(); i++)
790// auxGeom[i].parametrizeBasisBack(g1BasisVector[i]);
791 // END
792
793 // NEW
794 if (this->kindOfVertex() != 0)
795 computeKernel(g1BasisVector);
796
797 for (size_t i = 0; i < auxGeom.size(); i++)
798 auxGeom[i].parametrizeBasisBack(g1BasisVector[i]);
799 // END
800
801 for (size_t i = 0; i < auxGeom.size(); i++)
802 {
803 for (size_t ii = 0; ii < auxGeom[i].getG1Basis().nPatches(); ii++)
804 {
805 gsMatrix<T> coefs_temp;
806 coefs_temp.setZero(auxGeom[i].getG1Basis().patch(ii).coefs().rows(),1);
807 for (index_t j = 0; j < auxGeom[i].getG1Basis().patch(ii).coefs().rows(); j++)
808 {
809 if (!math::almostEqual(auxGeom[i].getG1Basis().patch(ii).coefs().at(j)*auxGeom[i].getG1Basis().patch(ii).coefs().at(j),T(0)))
810 coefs_temp(j,0) = auxGeom[i].getG1Basis().patch(ii).coefs()(j,0);
811 }
812 auxGeom[i].getG1Basis().patch(ii).setCoefs(coefs_temp);
813 }
814 basisVertexResult.push_back(auxGeom[i].getG1Basis());
815 }
816
817
818 // just for plotting
819// std::string fileName;
820// std::string basename = "VerticesBasisFunctions" + util::to_string(auxGeom.size());
821// gsParaviewCollection collection(basename);
822//
823// for (index_t np = 0; np < auxGeom.size(); ++np)
824// {
825// if (basisVertexResult.size() != 0)
826// for (index_t i = 0; i < basisVertexResult[np].nPatches(); ++i)
827// {
828// fileName = basename + "_" + util::to_string(np) + "_" + util::to_string(i);
829// gsField<T> temp_field(m_mp.patch(auxGeom[np].getGlobalPatchIndex()), basisVertexResult[np].patch(i));
830// gsWriteParaview(temp_field, fileName, 5000);
831// collection.addTimestep(fileName, i, "0.vts");
832//
833// }
834// }
835// collection.save();
836 }
837
838 gsG1AuxiliaryPatch<d,T> & getSinglePatch(const index_t i){ return auxGeom[i]; }
839
840 std::vector<gsMultiPatch<T>> getBasis()
841 { return basisVertexResult; }
842
843
844
845
846 void computeKernel(std::vector<gsMultiPatch<T>> & g1BasisVector)
847 {
848
849 gsMultiPatch<T> mp_vertex;
850 for(size_t i = 0; i < auxGeom.size(); i++)
851 mp_vertex.addPatch(auxGeom[i].getPatch());
852
853 mp_vertex.computeTopology();
854
855 index_t dim_mat = 0;
856 std::vector<index_t> dim_u, dim_v, side;
857 std::vector<index_t> dim_u_iFace;
858 std::vector<size_t> patchID, patchID_iFace;
859 gsMatrix<T> matrix_det(m_mp.targetDim(), m_mp.targetDim()), points(m_mp.parDim(),1);
860 points.setZero();
861 for(size_t np = 0; np < mp_vertex.nPatches(); np++)
862 {
863 if (mp_vertex.isBoundary(np,3)) // u
864 {
865 side.push_back(3);
866 patchID.push_back(np);
867 dim_u.push_back(auxGeom[np].getG1Basis().basis(0).component(0).size());
868 dim_v.push_back(auxGeom[np].getG1Basis().basis(0).component(1).size());
869 dim_mat += auxGeom[np].getG1Basis().basis(0).component(0).size();
870 dim_mat += auxGeom[np].getG1Basis().basis(0).component(0).size();
871
872 if(m_mp.parDim() == m_mp.targetDim()) // Planar
873 matrix_det.col(0) = auxGeom[np].getPatch().jacobian(points).col(0); // u
874 else if(m_mp.parDim() + 1 == m_mp.targetDim()) // Surface
875 {
876 gsMatrix<T> N, ev;
877 auxGeom[np].getPatch().jacobian_into(points, ev);
878 N.setZero(3,1);
879 N(0,0) = ev(1,0)*ev(2,1)-ev(2,0)*ev(1,1);
880 N(1,0) = ev(2,0)*ev(0,1)-ev(0,0)*ev(2,1);
881 N(2,0) = ev(0,0)*ev(1,1)-ev(1,0)*ev(0,1);
882 matrix_det.col(0) = ev.col(0);
883 matrix_det.col(2) = N;
884 }
885 }
886 if (mp_vertex.isBoundary(np,1)) // v
887 {
888 side.push_back(1);
889 patchID.push_back(np);
890 dim_u.push_back(auxGeom[np].getG1Basis().basis(0).component(0).size());
891 dim_v.push_back(auxGeom[np].getG1Basis().basis(0).component(1).size());
892 dim_mat += auxGeom[np].getG1Basis().basis(0).component(1).size();
893 dim_mat += auxGeom[np].getG1Basis().basis(0).component(1).size();
894
895 if(m_mp.parDim() == m_mp.targetDim()) // Planar
896 matrix_det.col(1) = auxGeom[np].getPatch().jacobian(points).col(1); // u
897 else if(m_mp.parDim() + 1 == m_mp.targetDim()) // Surface
898 {
899 gsMatrix<T> N, ev;
900 auxGeom[np].getPatch().jacobian_into(points, ev);
901 N.setZero(3,1);
902 N(0,0) = ev(1,0)*ev(2,1)-ev(2,0)*ev(1,1);
903 N(1,0) = ev(2,0)*ev(0,1)-ev(0,0)*ev(2,1);
904 N(2,0) = ev(0,0)*ev(1,1)-ev(1,0)*ev(0,1);
905 matrix_det.col(1) = ev.col(1);
906 }
907 }
908 if(mp_vertex.isInterface(patchSide(np,1)) && mp_vertex.isInterface(patchSide(np,3)))
909 {
910 dim_mat += 4;
911
912 patchID_iFace.push_back(np);
913 dim_u_iFace.push_back(auxGeom[np].getG1Basis().basis(0).component(0).size());
914 }
915
916 }
917 GISMO_ASSERT(patchID.size()==2,"Something went wrong");
918
919 index_t dofsCorner = 1;
920 if (matrix_det.determinant()*matrix_det.determinant() > 1e-15) // There is (numerically) a kink
921 {
922 dofsCorner = 0; // With Neumann
923 }
924
925 // gsDebug << "Det: " << matrix_det.determinant() << "\n";
926
927 gsMatrix<T> coefs_corner(dim_mat, 6);
928 coefs_corner.setZero();
929
930 index_t shift_row = 0;
931 for (size_t bdy_index = 0; bdy_index < patchID.size(); ++bdy_index)
932 {
933 if (side[bdy_index] < 3) // v
934 {
935 index_t shift_row_neumann = dim_v[bdy_index];
936 for (index_t i = 0; i < dim_v[bdy_index]; ++i)
937 {
938 for (index_t j = 0; j < 6; ++j)
939 {
940 T coef_temp = auxGeom[patchID[bdy_index]].getG1Basis().patch(j).coef(i*dim_u[bdy_index], 0); // v = 0
941 if (!math::almostEqual(coef_temp * coef_temp,T(0)))
942 coefs_corner(shift_row+i, j) = coef_temp;
943
944 coef_temp = auxGeom[patchID[bdy_index]].getG1Basis().patch(j).coef(i*dim_u[bdy_index] +1, 0); // v = 0
945 if (!math::almostEqual(coef_temp * coef_temp,T(0)))
946 coefs_corner(shift_row_neumann + shift_row+i, j) = coef_temp;
947
948 }
949 }
950 shift_row += dim_v[bdy_index];
951 shift_row += dim_v[bdy_index];
952 }
953 else // u
954 {
955 index_t shift_row_neumann = dim_u[bdy_index];
956 for (index_t i = 0; i < dim_u[bdy_index]; ++i)
957 {
958 for (index_t j = 0; j < 6; ++j)
959 {
960 T coef_temp = auxGeom[patchID[bdy_index]].getG1Basis().patch(j).coef(i, 0); // v = 0
961 if (!math::almostEqual(coef_temp * coef_temp,T(0)))
962 coefs_corner(shift_row+i, j) = coef_temp;
963
964 coef_temp = auxGeom[patchID[bdy_index]].getG1Basis().patch(j).coef(i+dim_u[bdy_index], 0); // v = 0
965 if (!math::almostEqual(coef_temp * coef_temp,T(0)))
966 coefs_corner(shift_row_neumann + shift_row+i, j) = coef_temp;
967
968 }
969 }
970 shift_row += dim_u[bdy_index];
971 shift_row += dim_u[bdy_index];
972 }
973 }
974 for (size_t iFace_index = 0; iFace_index < patchID_iFace.size(); ++iFace_index)
975 {
976 for (index_t i = 0; i < 2; ++i) // Only the first two
977 {
978 for (index_t j = 0; j < 6; ++j)
979 {
980 T coef_temp = auxGeom[patchID_iFace[iFace_index]].getG1Basis().patch(j).coef(i, 0); // v = 0
981 if (!math::almostEqual(coef_temp * coef_temp,T(0)))
982 coefs_corner(shift_row + i, j) = coef_temp;
983 coef_temp = auxGeom[patchID_iFace[iFace_index]].getG1Basis().patch(j).coef(i + dim_u_iFace[iFace_index],
984 0); // v = 0
985 if (!math::almostEqual(coef_temp * coef_temp,T(0)))
986 coefs_corner(shift_row + 2 + i, j) = coef_temp; // +2 bcs of the previous adding
987 }
988 }
989 shift_row += 4;
990 }
991
992 gsMatrix<T> kernel;
993 if (dofsCorner > 0)
994 {
995 T threshold = 1e-10;
996 gsEigen::FullPivLU<gsMatrix<T>> KernelCorner(coefs_corner);
997 KernelCorner.setThreshold(threshold);
998 //gsDebug << "Coefs: " << coefs_corner << "\n";
999 while (KernelCorner.dimensionOfKernel() < dofsCorner) {
1000 threshold += 1e-8;
1001 KernelCorner.setThreshold(threshold);
1002 }
1003 // gsDebug << "Dimension of Kernel: " << KernelCorner.dimensionOfKernel() << " With " << threshold << "\n";
1004
1005 gsMatrix<T> vertBas;
1006 vertBas.setIdentity(6, 6);
1007
1008 kernel = KernelCorner.kernel();
1009
1010 index_t count = 0;
1011 while (kernel.cols() < 6) {
1012 kernel.conservativeResize(kernel.rows(), kernel.cols() + 1);
1013 kernel.col(kernel.cols() - 1) = vertBas.col(count);
1014
1015 gsEigen::FullPivLU<gsMatrix<T>> ker_temp(kernel);
1016 ker_temp.setThreshold(1e-6);
1017 if (ker_temp.dimensionOfKernel() != 0) {
1018 kernel = kernel.block(0, 0, kernel.rows(), kernel.cols() - 1);
1019 }
1020 count++;
1021 }
1022 }
1023 else
1024 kernel.setIdentity(6, 6);
1025
1026 // gsDebug << "NumDofs: " << dofsCorner << " with Kernel: \n" << kernel << "\n";
1027
1028 for(size_t np = 0; np < auxGeom.size(); np++)
1029 {
1030 gsMultiPatch<T> temp_result_0 = auxGeom[np].getG1Basis();
1031
1032 for (index_t j = 0; j < 6; ++j)
1033 {
1034 index_t dim_uv = temp_result_0.basis(j).size();
1035 gsMatrix<T> coef_bf;
1036 coef_bf.setZero(dim_uv, 1);
1037 for (index_t i = 0; i < 6; ++i)
1038 if (!math::almostEqual(kernel(i, j) * kernel(i, j),T(0)))
1039 coef_bf += temp_result_0.patch(i).coefs() * kernel(i, j);
1040
1041
1042 g1BasisVector[np].patch(j).setCoefs(coef_bf);
1043 }
1044 }
1045
1046 }
1047
1048
1049protected:
1050
1051 std::vector<gsG1AuxiliaryPatch<d,T>> auxGeom;
1052 std::vector<index_t > auxVertexIndices;
1053 std::vector< std::vector<bool>> isBdy;
1054
1055 T sigma;
1056 index_t dim_kernel;
1057
1058 T m_zero;
1059
1060 // Store temp solution
1061 std::vector<gsMultiPatch<T>> basisVertexResult;
1062
1063 // only for plotting
1064 gsMultiPatch<T> & m_mp;
1065
1066}; // gsC1SurfVertex
1067
1068} // namespace gismo
gismo::gsBoxTopology::isInterface
bool isInterface(const patchSide &ps) const
Is the given patch side ps set to an interface?
Definition gsBoxTopology.cpp:103
gismo::gsBoxTopology::isBoundary
bool isBoundary(const patchSide &ps) const
Is the given patch side ps set to a boundary?
Definition gsBoxTopology.h:223
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition gsMatrix.h:211
gismo::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition gsMultiBasis.h:37
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition gsMultiPatch.h:100
gismo::gsMultiPatch::computeTopology
bool computeTopology(T tol=1e-4, bool cornersOnly=false, bool tjunctions=false)
Attempt to compute interfaces and boundaries automatically.
Definition gsMultiPatch.hpp:377
gismo::gsMultiPatch::patch
gsGeometry< T > & patch(size_t i) const
Return the i-th patch.
Definition gsMultiPatch.h:292
gismo::gsMultiPatch::addPatch
index_t addPatch(typename gsGeometry< T >::uPtr g)
Add a patch from a gsGeometry<T>::uPtr.
Definition gsMultiPatch.hpp:211
gismo::gsMultiPatch::basis
gsBasis< T > & basis(const size_t i) const
Return the basis of the i-th patch.
Definition gsMultiPatch.hpp:172
gismo::gsMultiPatch::nPatches
size_t nPatches() const
Number of patches.
Definition gsMultiPatch.h:274
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::maxDegree
short_t maxDegree() const
If the basis is of polynomial or piecewise polynomial type, then this function returns the maximum po...
Definition gsTensorBasis.h:470
gismo::gsTensorBasis::size
index_t size() const
Returns the number of elements in the basis.
Definition gsTensorBasis.h:108
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
index_t
#define index_t
Definition gsConfig.h:32
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsG1AuxiliaryPatch.h
Reparametrize one Patch.
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition gsExpressions.h:4486
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::patchSide
Struct which represents a certain side of a patch.
Definition gsBoundary.h:232