G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsParametrization.hpp
Go to the documentation of this file.
1 
15 #include <gsIO/gsOptionList.h>
16 #include <gsIO/gsWriteParaview.h>
19 
20 namespace gismo
21 {
22 
23 template<class T>
24 bool gsParametrization<T>::rangeCheck(const std::vector<index_t> &corners, const size_t minimum, const size_t maximum)
25 {
26  for (std::vector<index_t>::const_iterator it = corners.begin(); it != corners.end(); it++)
27  {
28  if ((size_t)*it < minimum || (size_t)*it > maximum)
29  { return false; }
30  }
31  return true;
32 }
33 
34 template<class T>
36 {
37  gsOptionList opt;
38  opt.addInt("boundaryMethod", "boundary methodes: {1:chords, 2:corners, 3:smallest, 4:restrict, 5:opposite, 6:distributed}", 4);
39  opt.addInt("parametrizationMethod", "parametrization methods: {1:shape, 2:uniform, 3:distance}", 1);
40  std::vector<index_t> corners;
41  opt.addMultiInt("corners", "vector for corners", corners);
42  opt.addReal("range", "in case of restrict or opposite", 0.1);
43  opt.addInt("number", "number of corners, in case of corners", 4);
44  opt.addReal("precision", "precision to calculate", 1E-8);
45  return opt;
46 }
47 
48 template<class T>
49 gsParametrization<T>::gsParametrization(const gsMesh<T> &mesh, const gsOptionList & list) : m_mesh(mesh)
50 {
51  m_options.update(list, gsOptionList::addIfUnknown);
52 }
53 
54 template<class T>
55 void gsParametrization<T>::calculate(const size_t boundaryMethod,
56  const size_t paraMethod,
57  const std::vector<index_t> &cornersInput,
58  const T rangeInput,
59  const size_t numberInput)
60 {
61  GISMO_ASSERT(boundaryMethod >= 1 && boundaryMethod <= 6,
62  "The boundary method " << boundaryMethod << " is not valid.");
63  GISMO_ASSERT(paraMethod >= 1 && paraMethod <= 3, "The parametrization method " << paraMethod << " is not valid.");
64  size_t n = m_mesh.getNumberOfInnerVertices();
65  size_t N = m_mesh.getNumberOfVertices();
66  size_t B = m_mesh.getNumberOfBoundaryVertices();
67  Neighbourhood neighbourhood(m_mesh, paraMethod);
68 
69  T w = 0;
70  std::vector<T> halfedgeLengths = m_mesh.getBoundaryChordLengths();
71  std::vector<index_t> corners;
72  std::vector<T> lengths;
73 
74  switch (boundaryMethod)
75  {
76  case 1:
77  m_parameterPoints.reserve(n + B);
78  for (size_t i = 1; i <= n + 1; i++)
79  {
80  m_parameterPoints.push_back(Point2D(0, 0, i));
81  }
82  for (size_t i = 0; i < B - 1; i++)
83  {
84  w += halfedgeLengths[i] * (T(1) / m_mesh.getBoundaryLength()) * (T)(4);
85  m_parameterPoints.push_back(Neighbourhood::findPointOnBoundary(w, n + i + 2));
86  }
87  break;
88  case 2:
89  corners = cornersInput;
90  case 3:
91  case 4:
92  case 5:
93  case 6: // N
94  if (boundaryMethod != 2)
95  corners = neighbourhood.getBoundaryCorners(boundaryMethod, rangeInput, numberInput);
96 
97  m_parameterPoints.reserve(N);
98  for (size_t i = 1; i <= N; i++)
99  {
100  m_parameterPoints.push_back(Point2D(0, 0, i));
101  }
102 
103  lengths = m_mesh.getCornerLengths(corners);
104  m_parameterPoints[n + corners[0] - 1] = Point2D(0, 0, n + corners[0]);
105 
106  for (size_t i = corners[0] + 1; i < corners[0] + B; i++)
107  {
108  w += halfedgeLengths[(i - 2) % B]
109  / findLengthOfPositionPart(i > B ? i - B : i, B, corners, lengths);
110  m_parameterPoints[(n + i - 1) > N - 1 ? n + i - 1 - B : n + i - 1] =
111  Neighbourhood::findPointOnBoundary(w, n + i > N ? n + i - B : n + i);
112  }
113  break;
114  default:
115  GISMO_ERROR("boundaryMethod not valid: " << boundaryMethod);
116  }
117 
118  constructAndSolveEquationSystem(neighbourhood, n, N);
119 }
120 
121 template<class T>
123  const size_t n,
124  const size_t N)
125 {
126  gsMatrix<T> A;
127  A.resize(n, n);
128  std::vector<T> lambdas;
129  gsVector<T> b1(n), b2(n);
130  b1.setZero(); b2.setZero();
131 
132  for (size_t i = 0; i < n; i++)
133  {
134  lambdas = neighbourhood.getLambdas(i);
135  for (size_t j = 0; j < n; j++)
136  {
137  A(i, j) = ( i==j ? (T)(1) : -lambdas[j] );
138  }
139 
140  for (size_t j = n; j < N; j++)
141  {
142  b1(i) += (lambdas[j]) * (m_parameterPoints[j][0]);
143  b2(i) += (lambdas[j]) * (m_parameterPoints[j][1]);
144  }
145  }
146 
147  gsVector<T> u(n), v(n);
148  gsEigen::PartialPivLU<typename gsMatrix<T>::Base> LU = A.partialPivLu();
149  u = LU.solve(b1);
150  v = LU.solve(b2);
151 
152  for (size_t i = 0; i < n; i++)
153  m_parameterPoints[i] << u(i), v(i);
154 }
155 
156 template<class T>
158 {
159  return m_parameterPoints[vertexIndex - 1];
160 }
161 
162 template<class T>
164 {
165  gsMatrix<T> m(2, m_mesh.getNumberOfVertices());
166  for (size_t i = 1; i <= m_mesh.getNumberOfVertices(); i++)
167  {
168  m.col(i - 1) << this->getParameterPoint(i)[0], this->getParameterPoint(i)[1];
169  }
170  return m;
171 }
172 
173 template<class T>
175 {
176  gsMatrix<T> m(3, m_mesh.getNumberOfVertices());
177  for (size_t i = 1; i <= m_mesh.getNumberOfVertices(); i++)
178  {
179  m.col(i - 1) << m_mesh.getVertex(i)->x(), m_mesh.getVertex(i)->y(), m_mesh.getVertex(i)->z();
180  }
181  return m;
182 }
183 
184 template<class T>
186 {
187  gsMesh<T> mesh;
188  mesh.reserve(3 * m_mesh.getNumberOfTriangles(), m_mesh.getNumberOfTriangles(), 0);
189  for (size_t i = 0; i < m_mesh.getNumberOfTriangles(); i++)
190  {
191  typename gsMesh<T>::VertexHandle v[3];
192  for (size_t j = 1; j <= 3; ++j)
193  {
194  v[j - 1] = mesh.addVertex(getParameterPoint(m_mesh.getGlobalVertexIndex(j, i))[0],
195  getParameterPoint(m_mesh.getGlobalVertexIndex(j, i))[1]);
196  }
197  mesh.addFace(v[0], v[1], v[2]);
198  }
199  return mesh.cleanMesh();
200 }
201 
202 
203 template <class T>
204 void gsParametrization<T>::writeTexturedMesh(std::string filename) const
205 {
206  gsMatrix<T> params(2, m_mesh.numVertices());
207 
208  for(size_t i=0; i<m_mesh.numVertices(); i++)
209  {
210  size_t index = m_mesh.unsorted(i);
211  params.col(i) = getParameterPoint(index).transpose();
212  }
213  gsWriteParaview(m_mesh, filename, params);
214 }
215 
216 template<class T>
218 {
219  m_options.update(list, gsOptionList::addIfUnknown);
220  return *this;
221 }
222 
223 template<class T>
225 {
226  calculate(m_options.getInt("boundaryMethod"),
227  m_options.getInt("parametrizationMethod"),
228  m_options.getMultiInt("corners"),
229  m_options.getReal("range"),
230  m_options.getInt("number"));
231 }
232 
233 template<class T>
234 T gsParametrization<T>::findLengthOfPositionPart(const size_t position,
235  const size_t numberOfPositions,
236  const std::vector<index_t> &bounds,
237  const std::vector<T> &lengths)
238 {
239  GISMO_UNUSED(numberOfPositions);
240  GISMO_ASSERT(1 <= position && position <= numberOfPositions, "The position " << position
241  << " is not a valid input. There are only " << numberOfPositions << " possible positions.");
242  GISMO_ASSERT(rangeCheck(bounds, 1, numberOfPositions), "The bounds are not a valid input. They have to be out of the possible positions, which only are "
243  << numberOfPositions << ". ");
244  size_t numberOfBounds = bounds.size();
245  size_t s = lengths.size();
246  if (position > (size_t)bounds[numberOfBounds - 1] || position <= (size_t)bounds[0])
247  return lengths[s - 1];
248  for (size_t i = 0; i < numberOfBounds; i++)
249  {
250  if (position - (size_t)bounds[0] + 1 > (size_t)bounds[i] - (size_t)bounds[0] + 1
251  && position - (size_t)bounds[0] + 1 <= (size_t)bounds[(i + 1) % numberOfBounds] - (size_t)bounds[0] + 1)
252  return lengths[i];
253  }
254  return 0;
255 }
256 
257 
258 //******************************************************************************************
259 //******************************* nested class Neighbourhood *******************************
260 //******************************************************************************************
261 
262 template<class T>
263 gsParametrization<T>::Neighbourhood::Neighbourhood(const gsHalfEdgeMesh<T> & meshInfo, const size_t parametrizationMethod) : m_basicInfos(meshInfo)
264 {
265  m_localParametrizations.reserve(meshInfo.getNumberOfInnerVertices());
266  for(size_t i=1; i <= meshInfo.getNumberOfInnerVertices(); i++)
267  {
268  m_localParametrizations.push_back(LocalParametrization(meshInfo, LocalNeighbourhood(meshInfo, i), parametrizationMethod));
269  }
270 
271  m_localBoundaryNeighbourhoods.reserve(meshInfo.getNumberOfVertices() - meshInfo.getNumberOfInnerVertices());
272  for(size_t i=meshInfo.getNumberOfInnerVertices()+1; i<= meshInfo.getNumberOfVertices(); i++)
273  {
274  m_localBoundaryNeighbourhoods.push_back(LocalNeighbourhood(meshInfo, i, 0));
275  }
276 }
277 
278 template<class T>
279 const std::vector<T>& gsParametrization<T>::Neighbourhood::getLambdas(const size_t i) const
280 {
281  return m_localParametrizations[i].getLambdas();
282 }
283 
284 template<class T>
285 const std::vector<index_t> gsParametrization<T>::Neighbourhood::getBoundaryCorners(const size_t method, const T range, const size_t number) const
286 {
287  std::vector<std::pair<T , size_t> > angles;
288  std::vector<index_t> corners;
289  angles.reserve(m_localBoundaryNeighbourhoods.size());
290  for(typename std::vector<LocalNeighbourhood>::const_iterator it=m_localBoundaryNeighbourhoods.begin(); it!=m_localBoundaryNeighbourhoods.end(); it++)
291  {
292  angles.push_back(std::pair<T , size_t>(it->getInnerAngle(), it->getVertexIndex() - m_basicInfos.getNumberOfInnerVertices()));
293  }
294  std::sort(angles.begin(), angles.end());
295  if(method == 3)
296  {
297  this->takeCornersWithSmallestAngles(4, angles, corners);
298  std::sort(corners.begin(), corners.end());
299  gsDebug << "According to the method 'smallest inner angles' the following corners were chosen:\n";
300  for(std::vector<index_t>::iterator it=corners.begin(); it!=corners.end(); it++)
301  {
302  gsDebug << (*it) << "\n";
303  }
304  }
305  else if(method == 5)
306  {
307  searchAreas(range, angles, corners);
308  gsDebug << "According to the method 'nearly opposite corners' the following corners were chosen:\n";
309  for(size_t i=0; i<corners.size(); i++)
310  {
311  gsDebug << corners[i] << "\n";
312  }
313  }
314  else if(method == 4)
315  {
316  bool flag = true;
317  corners.reserve(4);
318  corners.push_back(angles.front().second);
319  angles.erase(angles.begin());
320  while(corners.size() < 4)
321  {
322  flag = true;
323  for(std::vector<index_t>::iterator it=corners.begin(); it!=corners.end(); it++)
324  {
325  if(m_basicInfos.getShortestBoundaryDistanceBetween(angles.front().second, *it) < range*m_basicInfos.getBoundaryLength())
326  flag = false;
327  }
328  if(flag)
329  corners.push_back(angles.front().second);
330  angles.erase(angles.begin());
331  }
332  std::sort(corners.begin(), corners.end());
333  for(size_t i=0; i<corners.size(); i++)
334  {
335  gsDebug << corners[i] << "\n";
336  }
337  }
338  else if(method == 6)
339  {
340  T oldDifference = 0;
341  T newDifference = 0;
342  std::vector<index_t> newCorners;
343  std::vector<T> lengths;
344  angles.erase(angles.begin()+number, angles.end());
345  gsDebug << "Angles:\n";
346  for(size_t i=0; i<angles.size(); i++)
347  {
348  gsDebug << angles[i].first << ", " << angles[i].second << "\n";
349  }
350  newCorners.reserve((angles.size()*(angles.size()-1)*(angles.size()-2)*(angles.size()-3))/6);
351  corners.reserve((angles.size()*(angles.size()-1)*(angles.size()-2)*(angles.size()-3))/6);
352  for(size_t i=0; i<angles.size(); i++) // n
353  {
354  for(size_t j=i+1; j<angles.size(); j++) // * (n-1)/2
355  {
356  for(size_t k=j+1; k<angles.size(); k++) // * (n-2)/3
357  {
358  for(size_t l=k+1; l<angles.size(); l++) // * (n-3)/4
359  {
360  newCorners.push_back(angles[i].second);
361  newCorners.push_back(angles[j].second);
362  newCorners.push_back(angles[k].second);
363  newCorners.push_back(angles[l].second);
364  std::sort(newCorners.begin(), newCorners.end());
365  lengths = m_basicInfos.getCornerLengths(newCorners);
366  std::sort(lengths.begin(), lengths.end());
367  newDifference = math::abs(lengths[0] - lengths[3]);
368  if(oldDifference == 0 || newDifference < oldDifference)
369  {
370  corners.erase(corners.begin(), corners.end());
371  corners.push_back(angles[i].second);
372  corners.push_back(angles[j].second);
373  corners.push_back(angles[k].second);
374  corners.push_back(angles[l].second);
375  std::sort(corners.begin(), corners.end());
376  }
377  newCorners.erase(newCorners.begin(), newCorners.end());
378  oldDifference = newDifference;
379  }
380  }
381  }
382  }
383  gsDebug << "According to the method 'evenly distributed corners' the following corners were chosen:\n";
384  for(size_t i=0; i<corners.size(); i++)
385  {
386  gsDebug << corners[i] << "\n";
387  }
388  }
389  return corners;
390 }
391 
392 template<class T>
394 {
395  GISMO_ASSERT(0 <= w && w <= 4, "Wrong value for w.");
396  if(0 <= w && w <=1)
397  return Point2D(w,0, vertexIndex);
398  else if(1<w && w<=2)
399  return Point2D(1,w-(T)(1), vertexIndex);
400  else if(2<w && w<=3)
401  return Point2D(T(1)-w+(T)(2),1, vertexIndex);
402  else if(3<w && w<=4)
403  return Point2D(0,T(1)-w+(T)(3), vertexIndex);
404  return Point2D();
405 }
406 
407 //*****************************************************************************************************
408 //*****************************************************************************************************
409 //*******************THE******INTERN******FUNCTIONS******ARE******NOW******FOLLOWING*******************
410 //*****************************************************************************************************
411 //*****************************************************************************************************
412 
413 template<class T>
414 void gsParametrization<T>::Neighbourhood::takeCornersWithSmallestAngles(size_t number, std::vector<std::pair<T , size_t> >& sortedAngles, std::vector<index_t>& corners) const
415 {
416  sortedAngles.erase(sortedAngles.begin()+number, sortedAngles.end());
417 
418  corners.clear();
419  corners.reserve(sortedAngles.size());
420  for(typename std::vector<std::pair<T, size_t> >::iterator it=sortedAngles.begin(); it!=sortedAngles.end(); it++)
421  corners.push_back(it->second);
422 }
423 
424 template<class T>
425 std::vector<T> gsParametrization<T>::Neighbourhood::midpoints(const size_t numberOfCorners, const T length) const
426 {
427  std::vector<T> midpoints;
428  midpoints.reserve(numberOfCorners-1);
429  T n = 1./numberOfCorners;
430  for(size_t i=1; i<numberOfCorners; i++)
431  {
432  midpoints.push_back(T(i)*length*n);
433  }
434  return midpoints;
435 }
436 
437 template<class T>
438 void gsParametrization<T>::Neighbourhood::searchAreas(const T range, std::vector<std::pair<T, size_t> >& sortedAngles, std::vector<index_t>& corners) const
439 {
440  T l = m_basicInfos.getBoundaryLength();
441  std::vector<T> h = m_basicInfos.getBoundaryChordLengths();
442  this->takeCornersWithSmallestAngles(1,sortedAngles, corners);
443  std::vector<std::vector<std::pair<T , size_t> > > areas;
444  areas.reserve(3);
445  for(size_t i=0; i<3; i++)
446  {
447  areas.push_back(std::vector<std::pair<T , size_t> >());
448  }
449  std::vector<T> midpoints = this->midpoints(4, l);
450 
451  T walkAlong = 0;
452  for(size_t i=0; i<h.size(); i++)
453  {
454  walkAlong += h[(corners[0]+i-1)%h.size()];
455  for(int j = 2; j>=0; j--)
456  {
457  if(math::abs(walkAlong-midpoints[j]) <= l*range)
458  {
459  areas[j].push_back(std::pair<T , size_t>(m_localBoundaryNeighbourhoods[(corners[0]+i)%(h.size())].getInnerAngle(), (corners[0]+i)%h.size() + 1));
460  break;
461  }
462  }
463  }
464  std::sort(areas[0].begin(), areas[0].end());
465  std::sort(areas[1].begin(), areas[1].end());
466  std::sort(areas[2].begin(), areas[2].end());
467  bool smaller = false;
468  //corners.reserve(3);
469  for(size_t i=0; i<areas[0].size(); i++)
470  {
471  if(areas[0][i].second > (size_t)corners[0] || areas[0][i].second < (size_t)corners[0])
472  {
473  corners.push_back(areas[0][i].second);
474  if(areas[0][i].second < (size_t)corners[0])
475  {
476  smaller = true;
477  }
478  break;
479  }
480  }
481  for(size_t i=0; i<areas[1].size(); i++)
482  {
483  if(smaller)
484  {
485  if(areas[1][i].second > (size_t)corners[1] && areas[1][i].second < (size_t)corners[0])
486  {
487  corners.push_back(areas[1][i].second);
488  break;
489  }
490  }
491  else if(areas[1][i].second > (size_t)corners[1] || areas[1][i].second < (size_t)corners[0])
492  {
493  corners.push_back(areas[1][i].second);
494  if(areas[1][i].second < (size_t)corners[0])
495  {
496  smaller = true;
497  }
498  break;
499  }
500  }
501  for(size_t i=0; i<areas[2].size(); i++)
502  {
503  if(smaller)
504  {
505  if(areas[2][i].second > (size_t)corners[2] && areas[2][i].second < (size_t)corners[0])
506  {
507  corners.push_back(areas[2][i].second);
508  break;
509  }
510  }
511  else if(areas[2][i].second > (size_t)corners[2] || areas[2][i].second < (size_t)corners[0])
512  {
513  corners.push_back(areas[2][i].second);
514  break;
515  }
516  }
517 }
518 
519 //*******************************************************************************************
520 //**************************** nested class LocalParametrization ****************************
521 //*******************************************************************************************
522 
523 template<class T>
524 gsParametrization<T>::LocalParametrization::LocalParametrization(const gsHalfEdgeMesh<T>& meshInfo, const LocalNeighbourhood& localNeighbourhood, const size_t parametrizationMethod)
525 {
526  m_vertexIndex = localNeighbourhood.getVertexIndex();
527  std::list<size_t> indices = localNeighbourhood.getVertexIndicesOfNeighbours();
528  size_t d = localNeighbourhood.getNumberOfNeighbours();
529  switch (parametrizationMethod)
530  {
531  case 1:
532  {
533  std::list<T> angles = localNeighbourhood.getAngles();
534  VectorType points;
535  T theta = 0;
536  T nextAngle = 0;
537  for(typename std::list<T>::iterator it = angles.begin(); it!=angles.end(); ++it)
538  {
539  theta += *it;
540  }
541  Point2D p(0, 0, 0);
542  T length = (*meshInfo.getVertex(indices.front()) - *meshInfo.getVertex(m_vertexIndex)).norm();
543  Point2D nextPoint(length, 0, indices.front());
544  points.reserve(indices.size());
545  points.push_back(nextPoint);
546  gsVector<T> actualVector = nextPoint - p;
547  indices.pop_front();
548  T thetaInv = 1./theta;
549  gsVector<T> nextVector;
550  while(!indices.empty())
551  {
552  length = (*meshInfo.getVertex(indices.front()) - *meshInfo.getVertex(m_vertexIndex)).norm();
553  //length = (meshInfo.getVertex(indices.front()) - meshInfo.getVertex(m_vertexIndex) ).norm();
554  nextAngle = angles.front()*thetaInv * (T)(2) * (T)(EIGEN_PI);
555  nextVector = (gsEigen::Rotation2D<T>(nextAngle).operator*(actualVector).normalized()*length) + p;
556  nextPoint = Point2D(nextVector[0], nextVector[1], indices.front());
557  points.push_back(nextPoint);
558  actualVector = nextPoint - p;
559  angles.pop_front();
560  indices.pop_front();
561  }
562  calculateLambdas(meshInfo.getNumberOfVertices(), points);
563  }
564  break;
565  case 2:
566  m_lambdas.reserve(meshInfo.getNumberOfVertices());
567  for(size_t j=1; j <= meshInfo.getNumberOfVertices(); j++)
568  {
569  m_lambdas.push_back(0); // Lambda(m_vertexIndex, j, 0)
570  }
571  while(!indices.empty())
572  {
573  m_lambdas[indices.front()-1] += (1./d);
574  indices.pop_front();
575  }
576  break;
577  case 3:
578  {
579  std::list<T> neighbourDistances = localNeighbourhood.getNeighbourDistances();
580  T sumOfDistances = 0;
581  for(typename std::list<T>::iterator it = neighbourDistances.begin(); it != neighbourDistances.end(); it++)
582  {
583  sumOfDistances += *it;
584  }
585  T sumOfDistancesInv = 1./sumOfDistances;
586  m_lambdas.reserve(meshInfo.getNumberOfVertices());
587  for(size_t j=1; j <= meshInfo.getNumberOfVertices(); j++)
588  {
589  m_lambdas.push_back(0); //Lambda(m_vertexIndex, j, 0)
590  }
591  for(typename std::list<T>::iterator it = neighbourDistances.begin(); it != neighbourDistances.end(); it++)
592  {
593  m_lambdas[indices.front()-1] += ((*it)*sumOfDistancesInv);
594  indices.pop_front();
595  }
596  }
597  break;
598  default:
599  GISMO_ERROR("parametrizationMethod not valid: " << parametrizationMethod);
600  }
601 }
602 
603 template<class T>
605 {
606  return m_lambdas;
607 }
608 
609 //*****************************************************************************************************
610 //*****************************************************************************************************
611 //*******************THE******INTERN******FUNCTIONS******ARE******NOW******FOLLOWING*******************
612 //*****************************************************************************************************
613 //*****************************************************************************************************
614 
615 template<class T>
616 void gsParametrization<T>::LocalParametrization::calculateLambdas(const size_t N, VectorType& points)
617 {
618  m_lambdas.reserve(N);
619  for(size_t j=1; j <= N; j++)
620  {
621  m_lambdas.push_back(0); //Lambda(m_vertexIndex, j, 0)
622  }
623  Point2D p(0, 0, 0);
624  size_t d = points.size();
625  std::vector<T> my(d, 0);
626  size_t l=1;
627  size_t steps = 0;
628  //size_t checkOption = 0;
629  for(typename VectorType::const_iterator it=points.begin(); it != points.end(); it++)
630  {
631  gsLineSegment<2,T> actualLine(p, *it);
632  for(size_t i=1; i < d-1; i++)
633  {
634  if(l+i == d)
635  steps = d - 1;
636  else
637  steps = (l+i)%d - 1;
638  //checkoption is set to another number, in case mu is negativ
639  if(actualLine.intersectSegment(*(points.begin()+steps), *(points.begin()+(steps+1)%d)/*, checkOption*/))
640  {
641  //BarycentricCoordinates b(p, *it, *(points.begin()+steps), *(points.begin()+(steps+1)%d));
642  // calculating Barycentric Coordinates
643  gsMatrix<T, 3, 3> matrix;
644  matrix.topRows(2).col(0) = *it;
645  matrix.topRows(2).col(1) = *(points.begin()+steps);
646  matrix.topRows(2).col(2) = *(points.begin()+(steps+1)%d);
647  matrix.row(2).setOnes();
648 
649  gsVector3d<T> vector3d;
650  vector3d << p, 1;
651  gsVector3d<T> delta = matrix.partialPivLu().solve(vector3d);
652  my[l-1] = delta(0);
653  my[steps] = delta(1);
654  my[(steps + 1)%d] = delta(2);
655  break;
656  }
657  }
658  for(size_t k = 1; k <= d; k++)
659  {
660  m_lambdas[points[k-1].getVertexIndex()-1] += (my[k-1]);
661  }
662  std::fill(my.begin(), my.end(), 0);
663  l++;
664  }
665  for(typename std::vector<T>::iterator it=m_lambdas.begin(); it != m_lambdas.end(); it++)
666  {
667  *it /= d;
668  }
669  for(typename std::vector<T>::iterator it=m_lambdas.begin(); it != m_lambdas.end(); it++)
670  {
671  if(*it < 0)
672  gsInfo << *it << "\n";
673  }
674 }
675 
676 //*******************************************************************************************
677 //***************************** nested class LocalNeighbourhood *****************************
678 //*******************************************************************************************
679 
680 template<class T>
681 gsParametrization<T>::LocalNeighbourhood::LocalNeighbourhood(const gsHalfEdgeMesh<T>& meshInfo, const size_t vertexIndex, const bool innerVertex)
682 {
683  GISMO_ASSERT(!((innerVertex && vertexIndex > meshInfo.getNumberOfInnerVertices()) || vertexIndex < 1),
684  "Vertex with index " << vertexIndex << " does either not exist (< 1) or is not an inner vertex (> "
685  << meshInfo.getNumberOfInnerVertices() << ").");
686 
687  gsVertex<T> tmp;
688  m_vertexIndex = vertexIndex;
689  std::queue<typename gsHalfEdgeMesh<T>::Halfedge>
690  allHalfedges = meshInfo.getOppositeHalfedges(m_vertexIndex, innerVertex);
691  std::queue<typename gsHalfEdgeMesh<T>::Halfedge> nonFittingHalfedges;
692  m_neighbours.appendNextHalfedge(allHalfedges.front());
693  tmp = *meshInfo.getVertex(allHalfedges.front().getOrigin()) - *meshInfo.getVertex(m_vertexIndex);
694  m_angles.push_back(tmp.angle((*meshInfo.getVertex(allHalfedges.front().getEnd())
695  - *meshInfo.getVertex(vertexIndex))));
696  m_neighbourDistances.push_back(allHalfedges.front().getLength());
697  allHalfedges.pop();
698  while (!allHalfedges.empty())
699  {
700  if (m_neighbours.isAppendableAsNext(allHalfedges.front()))
701  {
702  m_neighbours.appendNextHalfedge(allHalfedges.front());
703  tmp = *meshInfo.getVertex(allHalfedges.front().getOrigin()) - *meshInfo.getVertex(m_vertexIndex);
704  m_angles
705  .push_back(tmp.angle((*meshInfo.getVertex(allHalfedges.front().getEnd())
706  - *meshInfo.getVertex(m_vertexIndex))));
707  m_neighbourDistances.push_back(allHalfedges.front().getLength());
708  allHalfedges.pop();
709  while (!nonFittingHalfedges.empty())
710  {
711  allHalfedges.push(nonFittingHalfedges.front());
712  nonFittingHalfedges.pop();
713  }
714  }
715  else if (m_neighbours.isAppendableAsPrev(allHalfedges.front()))
716  {
717  m_neighbours.appendPrevHalfedge(allHalfedges.front());
718  tmp = *meshInfo.getVertex(allHalfedges.front().getOrigin()) - *meshInfo.getVertex(m_vertexIndex);
719  m_angles
720  .push_front(tmp.angle((*meshInfo.getVertex(allHalfedges.front().getEnd())
721  - *meshInfo.getVertex(m_vertexIndex))));
722  m_neighbourDistances.push_back(allHalfedges.front().getLength());
723  allHalfedges.pop();
724  while (!nonFittingHalfedges.empty())
725  {
726  allHalfedges.push(nonFittingHalfedges.front());
727  nonFittingHalfedges.pop();
728  }
729  }
730  else
731  {
732  nonFittingHalfedges.push(allHalfedges.front());
733  allHalfedges.pop();
734  }
735  }
736 }
737 
738 template<class T>
740 {
741  return m_vertexIndex;
742 }
743 
744 template<class T>
746 {
747  return m_neighbours.getNumberOfVertices();
748 }
749 
750 template<class T>
752 {
753  return m_neighbours.getVertexIndices();
754 }
755 
756 template<class T>
758 {
759  return m_angles;
760 }
761 
762 template<class T>
764 {
765  T angle = 0;
766  for(typename std::list<T>::const_iterator it=m_angles.begin(); it!=m_angles.end(); it++)
767  {
768  angle += (*it);
769  }
770  return angle;
771 }
772 
773 template<class T>
775 {
776  return m_neighbourDistances;
777 }
778 
779 } // namespace gismo
Class that maintains parametrization This class Parametrization stores the mesh information and the t...
Definition: gsParametrization.h:46
const std::list< T > & getAngles() const
Get angles.
Definition: gsParametrization.hpp:757
A fixed-size, statically allocated 3D vector.
Definition: gsVector.h:218
const std::queue< Halfedge > getOppositeHalfedges(const size_t vertexIndex, const bool innerVertex=1) const
Returns queue of all opposite halfedges of vertex The opposite halfedge of a point in a triangle is m...
Definition: gsHalfEdgeMesh.hpp:189
Class that maintains neighbourhood information of triangle mesh. Represents the neighbourhood propert...
Definition: gsParametrization.h:285
size_t getVertexIndex() const
Get vertex index.
Definition: gsParametrization.hpp:739
#define gsDebug
Definition: gsDebug.h:61
void addMultiInt(const std::string &label, const std::string &desc, const std::vector< index_t > &values)
Adds an option-group gn containing values of a std::vector.
Definition: gsOptionList.cpp:221
const Point2D & getParameterPoint(size_t vertexIndex) const
Get parameter point Returns the parameter point with given vertex index.
Definition: gsParametrization.hpp:157
const std::vector< T > & getLambdas() const
Get lambdas The lambdas are returned.
Definition: gsParametrization.hpp:604
void writeTexturedMesh(std::string filename) const
Definition: gsParametrization.hpp:204
T getInnerAngle() const
Get inner angle.
Definition: gsParametrization.hpp:763
Class that maintains the local neighbourhood properties.
Definition: gsParametrization.h:127
bool intersectSegment(const gsPoint< d, T > &origin, const gsPoint< d, T > &end)
Tells wheter line intersects segment.
Definition: gsLineSegment.h:40
gsMatrix< T > createUVmatrix()
Definition: gsParametrization.hpp:163
Provides declaration of functions writing Paraview files.
std::list< T > getNeighbourDistances() const
Get neighbour distances.
Definition: gsParametrization.hpp:774
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
size_t getNumberOfInnerVertices() const
Get number of inner vertices The number of inner vertices of the triangle mesh is returned...
Definition: gsHalfEdgeMesh.hpp:69
LocalParametrization(const gsHalfEdgeMesh< T > &meshInfo, const LocalNeighbourhood &localNeighbourhood, const size_t parametrizationMethod=2)
Constructor Using this constructor one needs to input mesh information, a local neighbourhood and a p...
Definition: gsParametrization.hpp:524
void calculateLambdas(const size_t N, VectorType &points)
Calculate lambdas The lambdas according to Floater&#39;s algorithm are calculated.
Definition: gsParametrization.hpp:616
virtual gsMesh< T > createFlatMesh() const
Definition: gsParametrization.hpp:185
const std::list< size_t > getVertexIndicesOfNeighbours() const
Get vertex indices of neighbours.
Definition: gsParametrization.hpp:751
void addInt(const std::string &label, const std::string &desc, const index_t &value)
Adds a option named label, with description desc and value value.
Definition: gsOptionList.cpp:201
Class Representing a triangle mesh with 3D vertices.
Definition: gsMesh.h:31
Provides a list of labeled parameters/options that can be set and accessed easily.
size_t getNumberOfVertices() const
Get number of vertices The number of vertices of the triangle mesh is returned.
Definition: gsHalfEdgeMesh.hpp:57
gsMesh & cleanMesh()
reorders the vertices of all faces of an .stl mesh, such that only 1 vertex is used instead of #(adja...
Definition: gsMesh.hpp:250
Class maintains local parametrization This class represents a local parametrization for a point in th...
Definition: gsParametrization.h:221
const std::vector< T > & getLambdas(const size_t i) const
Get vector of lambdas.
Definition: gsParametrization.hpp:279
#define gsInfo
Definition: gsDebug.h:43
gsMatrix< T > createXYZmatrix()
Definition: gsParametrization.hpp:174
Class that maintains parametrization.
Neighbourhood(const gsHalfEdgeMesh< T > &meshInfo, const size_t parametrizationMethod=2)
Default constructor.
Definition: gsParametrization.hpp:263
void update(const gsOptionList &other, updateType type=ignoreIfUnknown)
Updates the object using the data from other.
Definition: gsOptionList.cpp:253
const gsMesh< T >::gsVertexHandle & getVertex(const size_t vertexIndex) const
Get vertex The vertex with index &#39;vertexIndex&#39; is returned.
Definition: gsHalfEdgeMesh.hpp:81
static gsOptionList defaultOptions()
Returns the list of default options for gsParametrization.
Definition: gsParametrization.hpp:35
A Point in T^d, with an index number.
Definition: gsPoint.h:26
Provides declaration of Line class.
void addReal(const std::string &label, const std::string &desc, const Real &value)
Adds a option named label, with description desc and value value.
Definition: gsOptionList.cpp:211
#define GISMO_UNUSED(x)
Definition: gsDebug.h:112
LocalNeighbourhood(const gsHalfEdgeMesh< T > &meshInfo, const size_t vertexIndex, const bool innerVertex=1)
Constructor.
Definition: gsParametrization.hpp:681
void constructAndSolveEquationSystem(const Neighbourhood &neighbourhood, const size_t n, const size_t N)
Constructs linear equation system and solves it.
Definition: gsParametrization.hpp:122
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
const std::vector< index_t > getBoundaryCorners(const size_t method, const T range=0.1, const size_t number=4) const
Get boundary corners depending on the method.
Definition: gsParametrization.hpp:285
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4488
Represents a line segment in d dimensions.
Definition: gsLineSegment.h:21
Class which holds a list of parameters/options, and provides easy access to them. ...
Definition: gsOptionList.h:32
virtual void compute()
Main function which performs the computation.
Definition: gsParametrization.hpp:224
gsHalfEdgeMesh is a gsMesh implementation that handles Halfedges
Definition: gsHalfEdgeMesh.h:46
size_t getNumberOfNeighbours() const
Get number of neighbours.
Definition: gsParametrization.hpp:745
gsParametrization(const gsMesh< T > &mesh, const gsOptionList &list=defaultOptions())
Constructor using the input mesh and (possibly) options.
Definition: gsParametrization.hpp:49
gsVertex class that represents a 3D vertex for a gsMesh.
Definition: gsVertex.h:26