G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsApproxC1Vertex.hpp
Go to the documentation of this file.
1 
14 #pragma once
15 
16 #include <gsUnstructuredSplines/src/gsApproxC1Vertex.h>
17 
18 namespace gismo
19 {
20 
21 template<short_t d,class T>
22 void gsApproxC1Vertex<d, T>::reparametrizeVertexPatches()
23 {
24  for(size_t i = 0; i < m_auxPatches.size(); i++)
25  {
26  checkOrientation(i); // Check if the orientation is correct. If not, modifies vertex and edge vectors
27 
28  index_t vertex_1 = m_vertexIndices[i];
29  if(m_auxPatches[i].getOrient() == 0) // not switched
30  switch (vertex_1) // == vertex
31  {
32  case 1:
33  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " not rotated\n";
34  break;
35  case 4:
36  m_auxPatches[i].rotateParamAntiClockTwice();
37  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " rotated twice anticlockwise\n";
38  break;
39  case 2:
40  m_auxPatches[i].rotateParamAntiClock();
41  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " rotated anticlockwise\n";
42  break;
43  case 3:
44  m_auxPatches[i].rotateParamClock();
45  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " rotated clockwise\n";
46  break;
47  }
48  else if (m_auxPatches[i].getOrient() == 1) // switched
49  switch (vertex_1) // == vertex
50  {
51  case 1:
52  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " not rotated\n";
53  break;
54  case 4:
55  m_auxPatches[i].rotateParamAntiClockTwice();
56  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " rotated twice anticlockwise\n";
57  break;
58  case 3:
59  m_auxPatches[i].rotateParamAntiClock();
60  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " rotated anticlockwise\n";
61  break;
62  case 2:
63  m_auxPatches[i].rotateParamClock();
64  //gsInfo << "Patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << " rotated clockwise\n";
65  break;
66  }
67  }
68 }
69 
70 template<short_t d,class T>
71 T gsApproxC1Vertex<d, T>::computeSigma(const std::vector<size_t> &vertexIndices)
72 {
73  T sigma = 0;
74 
75  T p = 0;
76  T h_geo = 1;
77  for(size_t i = 0; i < m_auxPatches.size(); i++)
78  {
79  const gsTensorBSplineBasis<2, T> bsp_temp = dynamic_cast<const gsTensorBSplineBasis<d, T>&>(m_auxPatches[0].getBasisRotated().piece(vertexIndices[i]+4));
80 
81  T p_temp = math::max(bsp_temp.degree(0), bsp_temp.degree(1));
82 
83  p = (p < p_temp ? p_temp : p);
84 
85  for(index_t j = 0; j < m_auxPatches[i].getPatchRotated().parDim(); j++)
86  {
87  T h_geo_temp = bsp_temp.knot(j,p + 1);
88  h_geo = (h_geo > h_geo_temp ? h_geo_temp : h_geo);
89  }
90  }
91 
92  gsMatrix<T> zero;
93  gsMatrix<T> deriv;
94  zero.setZero(2,1);
95  for (size_t i = 0; i < m_auxPatches.size(); i++)
96  {
97  deriv = m_auxPatches[i].getPatchRotated().deriv(zero);
98  sigma += deriv.template lpNorm<gsEigen::Infinity>();
99  }
100  sigma *= h_geo/(m_auxPatches.size()*p);
101 
102  return (1.0 / sigma);
103 }
104 
105 template<short_t d,class T>
106 void gsApproxC1Vertex<d, T>::checkOrientation(size_t i) {
107  if (m_auxPatches[i].getPatchRotated().orientation() == -1)
108  {
109  m_auxPatches[i].swapAxis();
110  //gsInfo << "Changed axis on patch: " << m_patchesAroundVertex[i] << " with side " << m_vertexIndices[i] << "\n";
111  }
112 }
113 
114 template<short_t d,class T>
115 void gsApproxC1Vertex<d, T>::computeKernel()
116 {
117 
118  gsMultiPatch<T> mp_vertex;
119  for(size_t i = 0; i < m_patchesAroundVertex.size(); i++)
120  mp_vertex.addPatch(m_auxPatches[i].getPatchRotated());
121 
122  mp_vertex.computeTopology();
123 
124  index_t dim_mat = 0;
125  std::vector<index_t> dim_u, dim_v, side, patchID;
126  std::vector<index_t> dim_u_iFace, patchID_iFace;
127  gsMatrix<T> matrix_det(m_mp.targetDim(), m_mp.targetDim()), points(m_mp.parDim(),1);
128  points.setZero();
129  for(size_t np = 0; np < mp_vertex.nPatches(); np++)
130  {
131  if (mp_vertex.isBoundary(np,3)) // u
132  {
133  side.push_back(3);
134  patchID.push_back(np);
135  dim_u.push_back(basisVertexResult[np].basis(0).component(0).size());
136  dim_v.push_back(basisVertexResult[np].basis(0).component(1).size());
137  dim_mat += basisVertexResult[np].basis(0).component(0).size();
138  if (!m_optionList.getSwitch("second"))
139  dim_mat += basisVertexResult[np].basis(0).component(0).size();
140 
141  if(m_mp.parDim() == m_mp.targetDim()) // Planar
142  matrix_det.col(0) = m_auxPatches[np].getPatchRotated().jacobian(points).col(0); // u
143  else if(m_mp.parDim() + 1 == m_mp.targetDim()) // Surface
144  {
145  gsMatrix<T> N, ev;
146  m_auxPatches[np].getPatchRotated().jacobian_into(points, ev);
147  N.setZero(3,1);
148  N(0,0) = ev(1,0)*ev(2,1)-ev(2,0)*ev(1,1);
149  N(1,0) = ev(2,0)*ev(0,1)-ev(0,0)*ev(2,1);
150  N(2,0) = ev(0,0)*ev(1,1)-ev(1,0)*ev(0,1);
151  matrix_det.col(0) = ev.col(0);
152  matrix_det.col(2) = N;
153  }
154  }
155  if (mp_vertex.isBoundary(np,1)) // v
156  {
157  side.push_back(1);
158  patchID.push_back(np);
159  dim_u.push_back(basisVertexResult[np].basis(0).component(0).size());
160  dim_v.push_back(basisVertexResult[np].basis(0).component(1).size());
161  dim_mat += basisVertexResult[np].basis(0).component(1).size();
162  if (!m_optionList.getSwitch("second"))
163  dim_mat += basisVertexResult[np].basis(0).component(1).size();
164 
165  if(m_mp.parDim() == m_mp.targetDim()) // Planar
166  matrix_det.col(1) = m_auxPatches[np].getPatchRotated().jacobian(points).col(1); // u
167  else if(m_mp.parDim() + 1 == m_mp.targetDim()) // Surface
168  {
169  gsMatrix<T> N, ev;
170  m_auxPatches[np].getPatchRotated().jacobian_into(points, ev);
171  N.setZero(3,1);
172  N(0,0) = ev(1,0)*ev(2,1)-ev(2,0)*ev(1,1);
173  N(1,0) = ev(2,0)*ev(0,1)-ev(0,0)*ev(2,1);
174  N(2,0) = ev(0,0)*ev(1,1)-ev(1,0)*ev(0,1);
175  matrix_det.col(1) = ev.col(1);
176  }
177  }
178  if(mp_vertex.isInterface(patchSide(np,1)) && mp_vertex.isInterface(patchSide(np,3)))
179  {
180  if (!m_optionList.getSwitch("second"))
181  dim_mat += 4;
182  else
183  dim_mat += 1;
184  patchID_iFace.push_back(np);
185  dim_u_iFace.push_back(basisVertexResult[np].basis(0).component(0).size());
186  }
187 
188  }
189  if (patchID.size() != 2)
190  gsInfo << "Something went wrong \n";
191 
192  index_t dofsCorner = 1;
193  if (m_optionList.getSwitch("second"))
194  dofsCorner = 3; // No Neumann
195 
196  if (matrix_det.determinant()*matrix_det.determinant() > 1e-15) // There is (numerically) a kink
197  {
198  if (!m_optionList.getSwitch("second"))
199  dofsCorner = 0; // With Neumann
200  else
201  dofsCorner = 1;
202  }
203 
204 
205  if (m_optionList.getSwitch("info"))
206  gsInfo << "Det: " << matrix_det.determinant() << "\n";
207 
208  gsMatrix<T> coefs_corner(dim_mat, 6);
209  coefs_corner.setZero();
210 
211  index_t shift_row = 0;
212  for (size_t bdy_index = 0; bdy_index < patchID.size(); ++bdy_index)
213  {
214  if (side[bdy_index] < 3) // v
215  {
216  index_t shift_row_neumann = dim_v[bdy_index];
217  for (index_t i = 0; i < dim_v[bdy_index]; ++i)
218  {
219  for (index_t j = 0; j < 6; ++j)
220  {
221  T coef_temp = basisVertexResult[patchID[bdy_index]].patch(j).coef(i*dim_u[bdy_index], 0); // v = 0
222  if (coef_temp * coef_temp > 1e-25)
223  coefs_corner(shift_row+i, j) = coef_temp;
224  if (!m_optionList.getSwitch("second"))
225  {
226  T coef_temp = basisVertexResult[patchID[bdy_index]].patch(j).coef(i*dim_u[bdy_index] +1, 0); // v = 0
227  if (coef_temp * coef_temp > 1e-25)
228  coefs_corner(shift_row_neumann + shift_row+i, j) = coef_temp;
229  }
230  }
231  }
232  shift_row += dim_v[bdy_index];
233  if (!m_optionList.getSwitch("second"))
234  shift_row += dim_v[bdy_index];
235  }
236  else // u
237  {
238  index_t shift_row_neumann = dim_u[bdy_index];
239  for (index_t i = 0; i < dim_u[bdy_index]; ++i)
240  {
241  for (index_t j = 0; j < 6; ++j)
242  {
243  T coef_temp = basisVertexResult[patchID[bdy_index]].patch(j).coef(i, 0); // v = 0
244  if (coef_temp * coef_temp > 1e-25)
245  coefs_corner(shift_row+i, j) = coef_temp;
246  if (!m_optionList.getSwitch("second"))
247  {
248  T coef_temp = basisVertexResult[patchID[bdy_index]].patch(j).coef(i+dim_u[bdy_index], 0); // v = 0
249  if (coef_temp * coef_temp > 1e-25)
250  coefs_corner(shift_row_neumann + shift_row+i, j) = coef_temp;
251  }
252  }
253  }
254  shift_row += dim_u[bdy_index];
255  if (!m_optionList.getSwitch("second"))
256  shift_row += dim_u[bdy_index];
257  }
258  }
259  for (size_t iFace_index = 0; iFace_index < patchID_iFace.size(); ++iFace_index)
260  {
261  if (!m_optionList.getSwitch("second")) {
262  for (index_t i = 0; i < 2; ++i) // Only the first two
263  {
264  for (index_t j = 0; j < 6; ++j) {
265  T coef_temp = basisVertexResult[patchID_iFace[iFace_index]].patch(j).coef(i, 0); // v = 0
266  if (coef_temp * coef_temp > 1e-25)
267  coefs_corner(shift_row + i, j) = coef_temp;
268  coef_temp = basisVertexResult[patchID_iFace[iFace_index]].patch(j).coef(i + dim_u_iFace[iFace_index],
269  0); // v = 0
270  if (coef_temp * coef_temp > 1e-25)
271  coefs_corner(shift_row + 2 + i, j) = coef_temp; // +2 bcs of the previous adding
272  }
273  }
274  shift_row += 4;
275  }
276  else
277  {
278  for (index_t j = 0; j < 6; ++j)
279  {
280  T coef_temp = basisVertexResult[patchID[iFace_index]].patch(j).coef(0, 0); // v = 0
281  if (coef_temp * coef_temp > 1e-25)
282  coefs_corner(shift_row, j) = coef_temp;
283  }
284  shift_row += 1;
285  }
286  }
287 
288  gsMatrix<T> kernel;
289  if (dofsCorner > 0)
290  {
291  T threshold = 1e-10;
292  gsEigen::FullPivLU<gsMatrix<T>> KernelCorner(coefs_corner);
293  KernelCorner.setThreshold(threshold);
294  //gsInfo << "Coefs: " << coefs_corner << "\n";
295  while (KernelCorner.dimensionOfKernel() < dofsCorner) {
296  threshold += 1e-8;
297  KernelCorner.setThreshold(threshold);
298  }
299  if (m_optionList.getSwitch("info"))
300  gsInfo << "Dimension of Kernel: " << KernelCorner.dimensionOfKernel() << " With " << threshold << "\n";
301 
302  gsMatrix<T> vertBas;
303  vertBas.setIdentity(6, 6);
304 
305  kernel = KernelCorner.kernel();
306 
307  size_t count = 0;
308  while (kernel.cols() < 6) {
309  kernel.conservativeResize(kernel.rows(), kernel.cols() + 1);
310  kernel.col(kernel.cols() - 1) = vertBas.col(count);
311 
312  gsEigen::FullPivLU<gsMatrix<T>> ker_temp(kernel);
313  ker_temp.setThreshold(1e-6);
314  if (ker_temp.dimensionOfKernel() != 0) {
315  kernel = kernel.block(0, 0, kernel.rows(), kernel.cols() - 1);
316  }
317  count++;
318  }
319  }
320  else
321  kernel.setIdentity(6, 6);
322 
323  if (m_optionList.getSwitch("info"))
324  gsInfo << "NumDofs: " << dofsCorner << " with Kernel: \n" << kernel << "\n";
325 
326  for(size_t np = 0; np < m_patchesAroundVertex.size(); np++)
327  {
328  gsMultiPatch<T> temp_result_0 = basisVertexResult[np];
329 
330  for (size_t j = 0; j < 6; ++j)
331  {
332  index_t dim_uv = temp_result_0.basis(j).size();
333  gsMatrix<T> coef_bf;
334  coef_bf.setZero(dim_uv, 1);
335  for (index_t i = 0; i < 6; ++i)
336  if (kernel(i, j) * kernel(i, j) > 1e-25)
337  coef_bf += temp_result_0.patch(i).coefs() * kernel(i, j);
338 
339  basisVertexResult[np].patch(j).setCoefs(coef_bf);
340  }
341  }
342 
343 }
344 
345 } // 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::gsTensorBasis::degree
short_t degree(short_t i) const
Returns the degree of the basis wrt variable i.
Definition: gsTensorBasis.h:465
gsApproxC1Vertex.h
Creates the (approx.) C1 Vertex space.
gismo::patchSide
Struct which represents a certain side of a patch.
Definition: gsBoundary.h:231
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::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
gsInfo
#define gsInfo
Definition: gsDebug.h:43
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition: gsMultiPatch.h:33
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