G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
21 namespace gismo {
22 template<short_t d, class T>
23 class gsC1SurfVertex
24 {
26  typedef memory::shared_ptr<gsC1SurfVertex> Ptr;
27 
29  typedef memory::unique_ptr<gsC1SurfVertex> uPtr;
30 
31 
32 public:
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 
1049 protected:
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::gsMultiPatch::addPatch
index_t addPatch(typename gsGeometry< T >::uPtr g)
Add a patch from a gsGeometry&lt;T&gt;::uPtr.
Definition: gsMultiPatch.hpp:210
gismo::patchSide
Struct which represents a certain side of a patch.
Definition: gsBoundary.h:231
gismo::gsTensorBasis::size
index_t size() const
Returns the number of elements in the basis.
Definition: gsTensorBasis.h:108
gismo::gsMultiPatch::basis
gsBasis< T > & basis(const size_t i) const
Return the basis of the i-th patch.
Definition: gsMultiPatch.hpp:171
gismo::gsBoxTopology::isBoundary
bool isBoundary(const patchSide &ps) const
Is the given patch side ps set to a boundary?
Definition: gsBoxTopology.h:223
index_t
#define index_t
Definition: gsConfig.h:32
gismo::gsTensorBSplineBasis
A tensor product B-spline basis.
Definition: gsTensorBSplineBasis.h:36
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
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::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:36
gismo::gsMultiPatch::patch
gsGeometry< T > & patch(size_t i) const
Return the i-th patch.
Definition: gsMultiPatch.h:226
gismo::gsMultiPatch::nPatches
size_t nPatches() const
Number of patches.
Definition: gsMultiPatch.h:208
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition: gsMultiPatch.h:33
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
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
gsG1AuxiliaryPatch.h
Reparametrize one Patch.
gismo::gsBoxTopology::isInterface
bool isInterface(const patchSide &ps) const
Is the given patch side ps set to an interface?
Definition: gsBoxTopology.cpp:103
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:366