gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
BackgroundMesh.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include "GmshMessage.h"
7 #include "BackgroundMesh.h"
8 #include "Numeric.h"
9 #include "Context.h"
10 #include "GVertex.h"
11 #include "GEdge.h"
12 #include "GFace.h"
13 #include "GModel.h"
14 #include "OS.h"
15 #include "Field.h"
16 #include "MElement.h"
17 #include "MElementOctree.h"
18 #include "MLine.h"
19 #include "MTriangle.h"
20 #include "MQuadrangle.h"
21 #include "MVertex.h"
22 
23 #if defined(HAVE_SOLVER)
24 #include "dofManager.h"
25 #include "laplaceTerm.h"
26 #include "linearSystemCSR.h"
27 #include "linearSystemFull.h"
28 #include "linearSystemPETSc.h"
29 #endif
30 
31 #if defined(HAVE_ANN)
32 static const int NBANN = 2;
33 #endif
34 
35 static const int MAX_THREADS = 256;
36 
37 std::vector<backgroundMesh *> backgroundMesh::_current =
38  std::vector<backgroundMesh *>(MAX_THREADS, (backgroundMesh *)nullptr);
39 
40 void backgroundMesh::set(GFace *gf)
41 {
42  int t = Msg::GetThreadNum();
43  if(t >= MAX_THREADS) {
44  Msg::Error("Maximum number of threads (%d) exceeded in background mesh",
45  MAX_THREADS);
46  return;
47  }
48  if(_current[t]) delete _current[t];
49  _current[t] = new backgroundMesh(gf);
50 }
51 
52 void backgroundMesh::setCrossFieldsByDistance(GFace *gf)
53 {
54  int t = Msg::GetThreadNum();
55  if(t >= MAX_THREADS) return;
56  if(_current[t]) delete _current[t];
57  _current[t] = new backgroundMesh(gf, true);
58 }
59 
60 void backgroundMesh::unset()
61 {
62  int t = Msg::GetThreadNum();
63  if(t >= MAX_THREADS) return;
64  if(_current[t]) delete _current[t];
65  _current[t] = nullptr;
66 }
67 
68 backgroundMesh *backgroundMesh::current()
69 {
70  int t = Msg::GetThreadNum();
71  if(t >= MAX_THREADS) return nullptr;
72  return _current[t];
73 }
74 
75 backgroundMesh::backgroundMesh(GFace *_gf, bool cfd)
76 #if defined(HAVE_ANN)
77  : _octree(nullptr), uv_kdtree(nullptr), nodes(nullptr), angle_nodes(nullptr),
78  angle_kdtree(nullptr)
79 #endif
80 {
81  if(cfd) {
82  Msg::Debug("Building cross field using closest distance");
83  propagateCrossFieldByDistance(_gf);
84  return;
85  }
86 
87  // create a bunch of triangles on the parametric space
88  // those triangles are local to the backgroundMesh so that
89  // they do not depend on the actual mesh that can be deleted
90 
91  std::set<SPoint2> myBCNodes;
92  for(std::size_t i = 0; i < _gf->triangles.size(); i++) {
93  MTriangle *e = _gf->triangles[i];
94  MVertex *news[3];
95  for(int j = 0; j < 3; j++) {
96  MVertex *v = e->getVertex(j);
97  auto it = _3Dto2D.find(v);
98  MVertex *newv = nullptr;
99  if(it == _3Dto2D.end()) {
100  SPoint2 p;
101  reparamMeshVertexOnFace(v, _gf, p);
102  newv = new MVertex(p.x(), p.y(), 0.0);
103  _vertices.push_back(newv);
104  _3Dto2D[v] = newv;
105  _2Dto3D[newv] = v;
106  if(v->onWhat()->dim() < 2) myBCNodes.insert(p);
107  }
108  else
109  newv = it->second;
110  news[j] = newv;
111  }
112  MTriangle *T2D = new MTriangle(news[0], news[1], news[2]);
113  _triangles.push_back(T2D);
114  }
115 
116 #if defined(HAVE_ANN)
117  index = new ANNidx[2];
118  dist = new ANNdist[2];
119  nodes = annAllocPts(myBCNodes.size(), 3);
120  auto itp = myBCNodes.begin();
121  int ind = 0;
122  while(itp != myBCNodes.end()) {
123  SPoint2 pt = *itp;
124  nodes[ind][0] = pt.x();
125  nodes[ind][1] = pt.y();
126  nodes[ind][2] = 0.0;
127  itp++;
128  ind++;
129  }
130  uv_kdtree = new ANNkd_tree(nodes, myBCNodes.size(), 3);
131 #endif
132 
133  // build a search structure
134  _octree = new MElementOctree(_triangles);
135 
136  // compute the mesh sizes at nodes
137  if(CTX::instance()->mesh.lcFromPoints) { propagate1dMesh(_gf); }
138  else {
139  auto itv2 = _2Dto3D.begin();
140  for(; itv2 != _2Dto3D.end(); ++itv2) {
141  _sizes[itv2->first] = CTX::instance()->mesh.lcMax;
142  }
143  }
144  // ensure that other criteria are fulfilled
145  updateSizes(_gf);
146 
147  // compute optimal mesh orientations
148  propagateCrossField(_gf);
149 
150  _3Dto2D.clear();
151  _2Dto3D.clear();
152 }
153 
154 backgroundMesh::~backgroundMesh()
155 {
156  for(std::size_t i = 0; i < _vertices.size(); i++) delete _vertices[i];
157  for(std::size_t i = 0; i < _triangles.size(); i++) delete _triangles[i];
158  if(_octree) delete _octree;
159 #if defined(HAVE_ANN)
160  if(uv_kdtree) delete uv_kdtree;
161  if(angle_kdtree) delete angle_kdtree;
162  if(nodes) annDeallocPts(nodes);
163  if(angle_nodes) annDeallocPts(angle_nodes);
164  delete[] index;
165  delete[] dist;
166 #endif
167 }
168 
169 static void propagateValuesOnFace(GFace *_gf,
170  std::map<MVertex *, double> &dirichlet,
171  simpleFunction<double> *ONE,
172  bool in_parametric_plane = false)
173 {
174 #if defined(HAVE_SOLVER)
175 #if defined(HAVE_PETSC)
176  linearSystemPETSc<double> *_lsys = new linearSystemPETSc<double>;
177 #elif defined(HAVE_GMM)
178  linearSystemCSRGmm<double> *_lsys = new linearSystemCSRGmm<double>;
179 #else
180  linearSystemFull<double> *_lsys = new linearSystemFull<double>;
181 #endif
182 
183  dofManager<double> myAssembler(_lsys);
184 
185  // fix boundary conditions
186  auto itv = dirichlet.begin();
187  for(; itv != dirichlet.end(); ++itv) {
188  myAssembler.fixVertex(itv->first, 0, 1, itv->second);
189  }
190 
191  // Number vertices
192  std::set<MVertex *> vs;
193  for(std::size_t k = 0; k < _gf->triangles.size(); k++)
194  for(int j = 0; j < 3; j++) vs.insert(_gf->triangles[k]->getVertex(j));
195  for(std::size_t k = 0; k < _gf->quadrangles.size(); k++)
196  for(int j = 0; j < 4; j++) vs.insert(_gf->quadrangles[k]->getVertex(j));
197 
198  std::map<MVertex *, SPoint3> theMap;
199  if(in_parametric_plane) {
200  for(auto it = vs.begin(); it != vs.end(); ++it) {
201  SPoint2 p;
202  reparamMeshVertexOnFace(*it, _gf, p);
203  theMap[*it] = SPoint3((*it)->x(), (*it)->y(), (*it)->z());
204  (*it)->setXYZ(p.x(), p.y(), 0.0);
205  }
206  }
207 
208  for(auto it = vs.begin(); it != vs.end(); ++it)
209  myAssembler.numberVertex(*it, 0, 1);
210 
211  // Assemble
212  laplaceTerm l(nullptr, 1, ONE);
213  for(std::size_t k = 0; k < _gf->triangles.size(); k++) {
214  MTriangle *t = _gf->triangles[k];
215  SElement se(t);
216  l.addToMatrix(myAssembler, &se);
217  }
218 
219  // Solve
220  if(myAssembler.sizeOfR()) { _lsys->systemSolve(); }
221 
222  // save solution
223  for(auto it = vs.begin(); it != vs.end(); ++it) {
224  myAssembler.getDofValue(*it, 0, 1, dirichlet[*it]);
225  }
226 
227  if(in_parametric_plane) {
228  for(auto it = vs.begin(); it != vs.end(); ++it) {
229  SPoint3 p = theMap[(*it)];
230  (*it)->setXYZ(p.x(), p.y(), p.z());
231  }
232  }
233  delete _lsys;
234 #endif
235 }
236 
237 void backgroundMesh::propagate1dMesh(GFace *_gf)
238 {
239  std::vector<GEdge *> const &e = _gf->edges();
240  auto it = e.begin();
241  std::map<MVertex *, double> sizes;
242 
243  for(; it != e.end(); ++it) {
244  if(!(*it)->isSeam(_gf)) {
245  for(std::size_t i = 0; i < (*it)->lines.size(); i++) {
246  MVertex *v1 = (*it)->lines[i]->getVertex(0);
247  MVertex *v2 = (*it)->lines[i]->getVertex(1);
248  if(v1 != v2) {
249  double d = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
250  (v1->y() - v2->y()) * (v1->y() - v2->y()) +
251  (v1->z() - v2->z()) * (v1->z() - v2->z()));
252  for(int k = 0; k < 2; k++) {
253  MVertex *v = (*it)->lines[i]->getVertex(k);
254  auto itv = sizes.find(v);
255  if(itv == sizes.end())
256  sizes[v] = log(d);
257  else
258  itv->second = 0.5 * (itv->second + log(d));
259  }
260  }
261  }
262  }
263  }
264 
265  simpleFunction<double> ONE(1.0);
266  propagateValuesOnFace(_gf, sizes, &ONE);
267 
268  auto itv2 = _2Dto3D.begin();
269  for(; itv2 != _2Dto3D.end(); ++itv2) {
270  MVertex *v_2D = itv2->first;
271  MVertex *v_3D = itv2->second;
272  _sizes[v_2D] = exp(sizes[v_3D]);
273  }
274 }
275 
276 crossField2d::crossField2d(MVertex *v, GEdge *ge)
277 {
278  double p;
279  bool success = reparamMeshVertexOnEdge(v, ge, p);
280  if(!success) {
281  Msg::Warning("cannot reparametrize a point in crossField");
282  _angle = 0;
283  return;
284  }
285  SVector3 t = ge->firstDer(p);
286  t.normalize();
287  _angle = atan2(t.y(), t.x());
288  crossField2d::normalizeAngle(_angle);
289 }
290 
291 void backgroundMesh::propagateCrossFieldByDistance(GFace *_gf)
292 {
293  std::vector<GEdge *> const &e = _gf->edges();
294  auto it = e.begin();
295  std::map<MVertex *, double> _cosines4, _sines4;
296  std::map<MVertex *, SPoint2> _param;
297 
298  for(; it != e.end(); ++it) {
299  if(!(*it)->isSeam(_gf)) {
300  for(std::size_t i = 0; i < (*it)->lines.size(); i++) {
301  MVertex *v[2];
302  v[0] = (*it)->lines[i]->getVertex(0);
303  v[1] = (*it)->lines[i]->getVertex(1);
304  SPoint2 p1, p2;
305  reparamMeshEdgeOnFace(v[0], v[1], _gf, p1, p2);
306  /* a correct way of computing angles */
307  Pair<SVector3, SVector3> der = _gf->firstDer((p1 + p2) * .5);
308  SVector3 t1 = der.first();
309  SVector3 t2(v[1]->x() - v[0]->x(), v[1]->y() - v[0]->y(),
310  v[1]->z() - v[0]->z());
311  t1.normalize();
312  t2.normalize();
313  double _angle = angle(t1, t2);
314  // double angle = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
315  crossField2d::normalizeAngle(_angle);
316  for(int i = 0; i < 2; i++) {
317  auto itc = _cosines4.find(v[i]);
318  auto its = _sines4.find(v[i]);
319  if(itc != _cosines4.end()) {
320  itc->second = 0.5 * (itc->second + cos(4 * _angle));
321  its->second = 0.5 * (its->second + sin(4 * _angle));
322  }
323  else {
324  _param[v[i]] = (i == 0) ? p1 : p2;
325  _cosines4[v[i]] = cos(4 * _angle);
326  _sines4[v[i]] = sin(4 * _angle);
327  }
328  }
329  }
330  }
331  }
332 
333 #if defined(HAVE_ANN)
334  index = new ANNidx[NBANN];
335  dist = new ANNdist[NBANN];
336  angle_nodes = annAllocPts(_cosines4.size(), 3);
337  auto itp = _cosines4.begin();
338  int ind = 0;
339  _sin.clear();
340  _cos.clear();
341  while(itp != _cosines4.end()) {
342  MVertex *v = itp->first;
343  double c = itp->second;
344  SPoint2 pt = _param[v];
345  double s = _sines4[v];
346  angle_nodes[ind][0] = pt.x();
347  angle_nodes[ind][1] = pt.y();
348  angle_nodes[ind][2] = 0.0;
349  _cos.push_back(c);
350  _sin.push_back(s);
351  itp++;
352  ind++;
353  }
354  angle_kdtree = new ANNkd_tree(angle_nodes, _cosines4.size(), 3);
355 #endif
356 }
357 
358 inline double myAngle(const SVector3 &a, const SVector3 &b, const SVector3 &d)
359 {
360  double cosTheta = dot(a, b);
361  double sinTheta = dot(crossprod(a, b), d);
362  return atan2(sinTheta, cosTheta);
363 }
364 
365 // smoothness = h * (|grad (cos 4 a)| + |grad (sin 4 a)|)
366 // smoothness is of order 1 if not smooth
367 // smoothness is of order h/L if smooth
368 // h --> mesh size
369 // L --> domain size
370 double backgroundMesh::getSmoothness(MElement *e)
371 {
372  MVertex *v0 = _3Dto2D[e->getVertex(0)];
373  MVertex *v1 = _3Dto2D[e->getVertex(1)];
374  MVertex *v2 = _3Dto2D[e->getVertex(2)];
375  auto i0 = _angles.find(v0);
376  auto i1 = _angles.find(v1);
377  auto i2 = _angles.find(v2);
378  double a[3] = {cos(4 * i0->second), cos(4 * i1->second), cos(4 * i2->second)};
379  double b[3] = {sin(4 * i0->second), sin(4 * i1->second), sin(4 * i2->second)};
380  double f[3];
381  e->interpolateGrad(a, 0, 0, 0, f);
382  const double gradcos = sqrt(f[0] * f[0] + f[1] * f[1] + f[2] * f[2]);
383  e->interpolateGrad(b, 0, 0, 0, f);
384  // const double gradsin = sqrt (f[0]*f[0]+f[1]*f[1]+f[2]*f[2]);
385  const double h = e->maxEdge();
386  return (gradcos /*+ gradsin*/) * h;
387 }
388 
389 double backgroundMesh::getSmoothness(double u, double v, double w)
390 {
391  if(!_octree) return 0.;
392  MElement *e = _octree->find(u, v, w, 2, true);
393  if(!e) return -1.0;
394  MVertex *v0 = e->getVertex(0);
395  MVertex *v1 = e->getVertex(1);
396  MVertex *v2 = e->getVertex(2);
397  auto i0 = _angles.find(v0);
398  auto i1 = _angles.find(v1);
399  auto i2 = _angles.find(v2);
400  double a[3] = {cos(4 * i0->second), cos(4 * i1->second), cos(4 * i2->second)};
401  double b[3] = {sin(4 * i0->second), sin(4 * i1->second), sin(4 * i2->second)};
402  double f[3];
403  e->interpolateGrad(a, 0, 0, 0, f);
404  const double gradcos = sqrt(f[0] * f[0] + f[1] * f[1] + f[2] * f[2]);
405  e->interpolateGrad(b, 0, 0, 0, f);
406  // const double gradsin = sqrt (f[0]*f[0]+f[1]*f[1]+f[2]*f[2]);
407  const double h = e->maxEdge();
408  return (gradcos /*+ gradsin*/) * h;
409 }
410 
411 void backgroundMesh::propagateCrossField(GFace *_gf)
412 {
413  propagateCrossFieldHJ(_gf);
414  // solve the non liear problem
415  constantPerElement<double> C;
416  int ITER = 0;
417  // int NSMOOTH = _gf->triangles.size();
418  while(0) {
419  // int NSMOOTH_NOW = 0;
420  for(std::size_t i = 0; i < _gf->triangles.size(); i++) {
421  double smoothness = getSmoothness(_gf->triangles[i]);
422  double val = smoothness < .5 ? 1.0 : 1.e-3; // exp(-absf/10);
423  C.set(_gf->triangles[i], val);
424  }
425  // if(NSMOOTH_NOW == NSMOOTH) break;
426  // NSMOOTH = NSMOOTH_NOW;
427  // break;
428  _angles.clear();
429  propagateCrossField(_gf, &C);
430  if(++ITER > 0) break;
431  }
432  // printf("converged in %d iterations\n", ITER);
433 #if 0
434  char name[256];
435  sprintf(name, "cross-%d-%d.pos", _gf->tag(), ITER);
436  print(name, 0, 1);
437  sprintf(name, "smooth-%d-%d.pos", _gf->tag(), ITER);
438  print(name, _gf, 2);
439 #endif
440 }
441 
442 void backgroundMesh::propagateCrossFieldHJ(GFace *_gf)
443 {
444  simpleFunction<double> ONE(1.0);
445  propagateCrossField(_gf, &ONE);
446 }
447 
448 void backgroundMesh::propagateCrossField(GFace *_gf,
449  simpleFunction<double> *ONE)
450 {
451  std::map<MVertex *, double> _cosines4, _sines4;
452  std::vector<GEdge *> const &e = _gf->edges();
453  auto it = e.begin();
454  for(; it != e.end(); ++it) {
455  if(!(*it)->isSeam(_gf)) {
456  for(std::size_t i = 0; i < (*it)->lines.size(); i++) {
457  MVertex *v[2];
458  v[0] = (*it)->lines[i]->getVertex(0);
459  v[1] = (*it)->lines[i]->getVertex(1);
460  SPoint2 p1, p2;
461  reparamMeshEdgeOnFace(v[0], v[1], _gf, p1, p2);
462  Pair<SVector3, SVector3> der = _gf->firstDer((p1 + p2) * .5);
463  SVector3 t1 = der.first();
464  SVector3 t2 = der.second();
465  SVector3 n = crossprod(t1, t2);
466  n.normalize();
467  SVector3 d1(v[1]->x() - v[0]->x(), v[1]->y() - v[0]->y(),
468  v[1]->z() - v[0]->z());
469  t1.normalize();
470  d1.normalize();
471  double _angle = myAngle(t1, d1, n);
472  crossField2d::normalizeAngle(_angle);
473  for(int i = 0; i < 2; i++) {
474  auto itc = _cosines4.find(v[i]);
475  auto its = _sines4.find(v[i]);
476  if(itc != _cosines4.end()) {
477  itc->second = 0.5 * (itc->second + cos(4 * _angle));
478  its->second = 0.5 * (its->second + sin(4 * _angle));
479  }
480  else {
481  _cosines4[v[i]] = cos(4 * _angle);
482  _sines4[v[i]] = sin(4 * _angle);
483  }
484  }
485  }
486  }
487  }
488 
489  propagateValuesOnFace(_gf, _cosines4, ONE, false);
490  propagateValuesOnFace(_gf, _sines4, ONE, false);
491 
492  // print("cos4.pos",0,_cosines4,0);
493  // print("sin4.pos",0,_sines4,0);
494 
495  auto itv2 = _2Dto3D.begin();
496  for(; itv2 != _2Dto3D.end(); ++itv2) {
497  MVertex *v_2D = itv2->first;
498  MVertex *v_3D = itv2->second;
499  double angle = atan2(_sines4[v_3D], _cosines4[v_3D]) / 4.0;
500  crossField2d::normalizeAngle(angle);
501  _angles[v_2D] = angle;
502  }
503 }
504 
505 void backgroundMesh::updateSizes(GFace *_gf)
506 {
507  auto itv = _sizes.begin();
508  for(; itv != _sizes.end(); ++itv) {
509  SPoint2 p;
510  MVertex *v = _2Dto3D[itv->first];
511  double lc;
512  if(v->onWhat()->dim() == 0) {
513  lc = BGM_MeshSize(v->onWhat(), 0, 0, v->x(), v->y(), v->z());
514  }
515  else if(v->onWhat()->dim() == 1) {
516  double u;
517  v->getParameter(0, u);
518  lc = BGM_MeshSize(v->onWhat(), u, 0, v->x(), v->y(), v->z());
519  }
520  else {
521  reparamMeshVertexOnFace(v, _gf, p);
522  lc = BGM_MeshSize(_gf, p.x(), p.y(), v->x(), v->y(), v->z());
523  }
524  itv->second = std::min(lc, itv->second);
525  itv->second = std::max(itv->second, CTX::instance()->mesh.lcMin);
526  itv->second = std::min(itv->second, CTX::instance()->mesh.lcMax);
527  }
528  // do not allow large variations in the size field
529  // (Int. J. Numer. Meth. Engng. 43, 1143-1165 (1998) MESH GRADATION
530  // CONTROL, BOROUCHAKI, HECHT, FREY)
531  std::set<MEdge, MEdgeLessThan> edges;
532  for(std::size_t i = 0; i < _triangles.size(); i++) {
533  for(int j = 0; j < _triangles[i]->getNumEdges(); j++) {
534  edges.insert(_triangles[i]->getEdge(j));
535  }
536  }
537  const double _beta = 1.3;
538  for(int i = 0; i < 3; i++) {
539  auto it = edges.begin();
540  for(; it != edges.end(); ++it) {
541  MVertex *v0 = it->getVertex(0);
542  MVertex *v1 = it->getVertex(1);
543  MVertex *V0 = _2Dto3D[v0];
544  MVertex *V1 = _2Dto3D[v1];
545  auto s0 = _sizes.find(V0);
546  auto s1 = _sizes.find(V1);
547  if(s0 == _sizes.end() || s1 == _sizes.end()) {
548  // Note: this Warning is disabled because it's too verbose in logs
549  // is this behavior normal ? Was doing undefined behavior before
550  // (detected by valgrind)
551  // Msg::Warning("in backgroundMesh::updateSizes(), vertex not found in maps");
552  continue;
553  }
554  if(s0->second < s1->second)
555  s1->second = std::min(s1->second, _beta * s0->second);
556  else
557  s0->second = std::min(s0->second, _beta * s1->second);
558  }
559  }
560 }
561 
562 bool backgroundMesh::inDomain(double u, double v, double w) const
563 {
564  if(!_octree) return false;
565  return _octree->find(u, v, w, 2, true) != nullptr;
566 }
567 
568 double backgroundMesh::operator()(double u, double v, double w) const
569 {
570  if(!_octree) {
571  Msg::Error("No octree in background mesh");
572  return 0.;
573  }
574  double uv[3] = {u, v, w};
575  double uv2[3];
576  MElement *e = _octree->find(u, v, w, 2, true);
577  if(!e) {
578 #if defined(HAVE_ANN)
579  if(uv_kdtree->nPoints() < 2) return -1000.;
580  double pt[3] = {u, v, 0.0};
581 #pragma omp critical // just to avoid crash (still incorrect) - should use nanoflann
582  uv_kdtree->annkSearch(pt, 2, index, dist);
583  SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
584  SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
585  SPoint3 pnew;
586  double d;
587  signedDistancePointLine(p1, p2, SPoint3(u, v, 0.), d, pnew);
588  e = _octree->find(pnew.x(), pnew.y(), 0.0, 2, true);
589 #endif
590  if(!e) {
591  Msg::Error("BGM octree: cannot find UVW=%g %g %g", u, v, w);
592  return -1000.0; // 0.4;
593  }
594  }
595  e->xyz2uvw(uv, uv2);
596  auto itv1 = _sizes.find(e->getVertex(0));
597  auto itv2 = _sizes.find(e->getVertex(1));
598  auto itv3 = _sizes.find(e->getVertex(2));
599  return itv1->second * (1 - uv2[0] - uv2[1]) + itv2->second * uv2[0] +
600  itv3->second * uv2[1];
601 }
602 
603 double backgroundMesh::getAngle(double u, double v, double w) const
604 {
605  // use closest point for computing cross field angles: this allows NOT to
606  // generate a spurious mesh and solve a PDE
607  if(!_octree) {
608 #if defined(HAVE_ANN)
609  double angle = 0.;
610  if(angle_kdtree->nPoints() >= NBANN) {
611  double pt[3] = {u, v, 0.0};
612 #pragma omp critical // just to avoid crash (still incorrect) - should use nanoflann
613  angle_kdtree->annkSearch(pt, NBANN, index, dist);
614  double SINE = 0.0, COSINE = 0.0;
615  for(int i = 0; i < NBANN; i++) {
616  SINE += _sin[index[i]];
617  COSINE += _cos[index[i]];
618  }
619  angle = atan2(SINE, COSINE) / 4.0;
620  }
621  crossField2d::normalizeAngle(angle);
622  return angle;
623 #endif
624  }
625 
626  // HACK FOR LEWIS
627  // h = 1+30(y-x^2)^2 + (1-x)^2
628  // double x = u;
629  // double y = v;
630  // double dhdx = 30 * 2 * (y-x*x) * (-2*x) - 2 * (1-x);
631  // double dhdy = 30 * 2 * (y-x*x);
632  // double angles = atan2(y,x*x);
633  // crossField2d::normalizeAngle (angles);
634  // return angles;
635 
636  double uv[3] = {u, v, w};
637  double uv2[3];
638  MElement *e = _octree->find(u, v, w, 2, true);
639  if(!e) {
640 #if defined(HAVE_ANN)
641  if(uv_kdtree->nPoints() < 2) return -1000.0;
642  double pt[3] = {u, v, 0.0};
643 #pragma omp critical // just to avoid crash (still incorrect) - should use nanoflann
644  uv_kdtree->annkSearch(pt, 2, index, dist);
645  SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
646  SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
647  SPoint3 pnew;
648  double d;
649  signedDistancePointLine(p1, p2, SPoint3(u, v, 0.), d, pnew);
650  e = _octree->find(pnew.x(), pnew.y(), 0., 2, true);
651 #endif
652  if(!e) {
653  Msg::Error("BGM octree angle: cannot find UVW=%g %g %g", u, v, w);
654  return -1000.0;
655  }
656  }
657  e->xyz2uvw(uv, uv2);
658  auto itv1 = _angles.find(e->getVertex(0));
659  auto itv2 = _angles.find(e->getVertex(1));
660  auto itv3 = _angles.find(e->getVertex(2));
661 
662  double cos4 = cos(4 * itv1->second) * (1 - uv2[0] - uv2[1]) +
663  cos(4 * itv2->second) * uv2[0] + cos(4 * itv3->second) * uv2[1];
664  double sin4 = sin(4 * itv1->second) * (1 - uv2[0] - uv2[1]) +
665  sin(4 * itv2->second) * uv2[0] + sin(4 * itv3->second) * uv2[1];
666  double angle = atan2(sin4, cos4) / 4.0;
667  crossField2d::normalizeAngle(angle);
668 
669  return angle;
670 }
671 
672 void backgroundMesh::print(const std::string &filename, GFace *gf,
673  const std::map<MVertex *, double> &_whatToPrint,
674  int smooth)
675 {
676  FILE *f = Fopen(filename.c_str(), "w");
677  if(!f) {
678  Msg::Error("Could not open file '%s'", filename.c_str());
679  return;
680  }
681  fprintf(f, "View \"Background Mesh\"{\n");
682  if(smooth) {
683  for(std::size_t i = 0; i < gf->triangles.size(); i++) {
684  MVertex *v1 = gf->triangles[i]->getVertex(0);
685  MVertex *v2 = gf->triangles[i]->getVertex(1);
686  MVertex *v3 = gf->triangles[i]->getVertex(2);
687  double x = getSmoothness(gf->triangles[i]);
688  fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n", v1->x(),
689  v1->y(), v1->z(), v2->x(), v2->y(), v2->z(), v3->x(), v3->y(),
690  v3->z(), x, x, x);
691  }
692  }
693  else {
694  for(std::size_t i = 0; i < _triangles.size(); i++) {
695  MVertex *v1 = _triangles[i]->getVertex(0);
696  MVertex *v2 = _triangles[i]->getVertex(1);
697  MVertex *v3 = _triangles[i]->getVertex(2);
698  auto itv1 = _whatToPrint.find(v1);
699  auto itv2 = _whatToPrint.find(v2);
700  auto itv3 = _whatToPrint.find(v3);
701  if(!gf) {
702  fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n", v1->x(),
703  v1->y(), v1->z(), v2->x(), v2->y(), v2->z(), v3->x(), v3->y(),
704  v3->z(), itv1->second, itv2->second, itv3->second);
705  }
706  else {
707  GPoint p1 = gf->point(SPoint2(v1->x(), v1->y()));
708  GPoint p2 = gf->point(SPoint2(v2->x(), v2->y()));
709  GPoint p3 = gf->point(SPoint2(v3->x(), v3->y()));
710  fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n", p1.x(),
711  p1.y(), p1.z(), p2.x(), p2.y(), p2.z(), p3.x(), p3.y(), p3.z(),
712  itv1->second, itv2->second, itv3->second);
713  }
714  }
715  }
716  fprintf(f, "};\n");
717  fclose(f);
718 }
719 
720 MElement *backgroundMesh::getMeshElementByCoord(double u, double v, double w,
721  bool strict)
722 {
723  if(!_octree) {
724  Msg::Debug("Rebuilding BackgroundMesh element octree");
725  _octree = new MElementOctree(_triangles);
726  }
727  return _octree->find(u, v, w, 2, strict);
728 }
729 
730 /* Global variable instanciation */
731 std::vector<std::unique_ptr<GlobalBackgroundMesh> > global_bmeshes;
732 
733 bool backgroudMeshExists(const std::string &name)
734 {
735  for(size_t i = 0; i < global_bmeshes.size(); ++i) {
736  if(global_bmeshes[i]->name == name) return true;
737  }
738  return false;
739 }
740 
741 GlobalBackgroundMesh &getBackgroundMesh(const std::string &name)
742 {
743  for(size_t i = 0; i < global_bmeshes.size(); ++i) {
744  if(global_bmeshes[i]->name == name) return *(global_bmeshes[i]);
745  }
746  global_bmeshes.push_back(
747  std::unique_ptr<GlobalBackgroundMesh>(new GlobalBackgroundMesh(name)));
748  return *(global_bmeshes.back());
749 }
750 
751 GlobalBackgroundMesh::~GlobalBackgroundMesh()
752 {
753  Msg::Debug("GlobalBackgroundMesh destructor call");
754  for(MVertex *v : mesh_vertices)
755  if(v) {
756  delete v;
757  v = nullptr;
758  }
759  mesh_vertices.clear();
760  edgeBackgroundMeshes.clear();
761  faceBackgroundMeshes.clear();
762 }
763 
764 int GlobalBackgroundMesh::importGModelMeshes(GModel *_gm,
765  bool overwriteExisting)
766 {
767  Msg::Debug("GlobalBackgroundMesh: import GModel mesh ...");
768  gm = _gm;
769 
770  if(overwriteExisting && mesh_vertices.size() > 0) { /* Clear mesh */
771  Msg::Debug("- delete existing mesh (for overwrite)");
772  for(MVertex *v : mesh_vertices)
773  if(v) {
774  delete v;
775  v = nullptr;
776  }
777  mesh_vertices.clear();
778  edgeBackgroundMeshes.clear();
779  faceBackgroundMeshes.clear();
780  Msg::Debug("- import mesh from GModel");
781  }
782 
783  std::unordered_map<MVertex *, MVertex *> old2new;
784 
785  std::set<GVertex *, GEntityPtrLessThan> vertices = gm->getVertices();
786  std::set<GEdge *, GEntityPtrLessThan> edges = gm->getEdges();
787  for(GFace *gf : gm->getFaces()) {
788  for(GVertex *e : gf->embeddedVertices()) vertices.insert(e);
789  for(GEdge *e : gf->embeddedEdges()) edges.insert(e);
790  }
791 
792  /* MVertex from GVertex */
793  mesh_vertices.reserve(mesh_vertices.size() + vertices.size());
794  for(GVertex *gv : vertices) {
795  for(size_t i = 0; i < gv->mesh_vertices.size(); ++i) {
796  MVertex *v = gv->mesh_vertices[i];
797  auto it = old2new.find(v);
798  if(it == old2new.end()) {
799  mesh_vertices.push_back(new MVertex(v->x(), v->y(), v->z(), gv));
800  old2new[v] = mesh_vertices.back();
801  }
802  }
803  }
804 
805  /* MLine from GEdge */
806  size_t nlines = 0;
807  for(GEdge *ge : edges) {
808  edgeBackgroundMeshes[ge] = BackgroundMeshGEdge();
809  BackgroundMeshGEdge &bm = edgeBackgroundMeshes[ge];
810  bm.ge = ge;
811  for(size_t i = 0; i < ge->mesh_vertices.size(); ++i) {
812  MVertex *v = ge->mesh_vertices[i];
813  double t = 0.;
814  v->setParameter(0, t);
815  auto it = old2new.find(v);
816  if(it == old2new.end()) {
817  mesh_vertices.push_back(new MEdgeVertex(v->x(), v->y(), v->z(), ge, t));
818  old2new[v] = mesh_vertices.back();
819  }
820  }
821  bm.lines.reserve(ge->lines.size());
822  for(size_t i = 0; i < ge->lines.size(); ++i) {
823  MElement *elt = ge->lines[i];
824  MVertex *v0 = old2new[elt->getVertex(0)];
825  MVertex *v1 = old2new[elt->getVertex(1)];
826  if(v0 == NULL || v1 == NULL) {
827  Msg::Error("GlobalBackgroundMesh: failed to import mesh (1)");
828  // Note: some vertex are stored outside the GModel ?
829  return -1;
830  }
831  bm.lines.push_back(MLine(v0, v1));
832  nlines += 1;
833  }
834  }
835 
836  /* MTriangle from GFace */
837  size_t ntris = 0;
838  for(GFace *gf : gm->getFaces()) {
839  faceBackgroundMeshes[gf] = BackgroundMeshGFace();
840  BackgroundMeshGFace &bm = faceBackgroundMeshes[gf];
841  bm.gf = gf;
842  /* Vertices */
843  for(size_t i = 0; i < gf->mesh_vertices.size(); ++i) {
844  MVertex *v = gf->mesh_vertices[i];
845  auto it = old2new.find(v);
846  if(it == old2new.end()) {
847  double uv0 = 0.;
848  double uv1 = 0.;
849  v->getParameter(0, uv0);
850  v->getParameter(1, uv1);
851  mesh_vertices.push_back(
852  new MFaceVertex(v->x(), v->y(), v->z(), gf, uv0, uv1));
853  old2new[v] = mesh_vertices.back();
854  }
855  }
856 
857  bm.triangles.reserve(gf->triangles.size() + 2 * gf->quadrangles.size());
858  /* Triangles */
859  for(size_t i = 0; i < gf->triangles.size(); ++i) {
860  MElement *elt = gf->triangles[i];
861  MVertex *v0 = old2new[elt->getVertex(0)];
862  MVertex *v1 = old2new[elt->getVertex(1)];
863  MVertex *v2 = old2new[elt->getVertex(2)];
864  if(v0 == NULL || v1 == NULL || v2 == NULL) {
865  Msg::Error("GlobalBackgroundMesh: failed to import mesh (2)");
866  // Note: some vertex are stored outside the GModel ?
867  return -1;
868  }
869  bm.triangles.push_back(MTriangle(v0, v1, v2));
870  ntris += 1;
871  }
872 
873  /* Quads */
874  if(gf->quadrangles.size() > 0) {
875  for(size_t i = 0; i < gf->quadrangles.size(); ++i) {
876  MElement *elt = gf->quadrangles[i];
877  MVertex *v0 = old2new[elt->getVertex(0)];
878  MVertex *v1 = old2new[elt->getVertex(1)];
879  MVertex *v2 = old2new[elt->getVertex(2)];
880  MVertex *v3 = old2new[elt->getVertex(3)];
881  if(v0 == NULL || v1 == NULL || v2 == NULL || v3 == NULL) {
882  Msg::Error("GlobalBackgroundMesh: failed to import mesh (3)");
883  // Note: some vertex are stored outside the GModel ?
884  return -1;
885  }
886  bm.triangles.push_back(MTriangle(v0, v1, v2));
887  bm.triangles.push_back(MTriangle(v0, v2, v3));
888  ntris += 2;
889  }
890  }
891  }
892 
893  Msg::Debug("- imported: %li vertices, %li lines, %li triangles",
894  old2new.size(), nlines, ntris);
895 
896  return 0;
897 }
backgroundMesh::getMeshElementByCoord
MElement * getMeshElementByCoord(double u, double v, double w, bool strict=true)
Definition: BackgroundMesh.cpp:720
crossprod
SVector3 crossprod(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:150
backgroundMesh::_sizes
std::map< MVertex *, double > _sizes
Definition: BackgroundMesh.h:50
MElement::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MElement.cpp:1128
MTriangle.h
linearSystemFull
Definition: linearSystemFull.h:17
backgroundMesh::_angles
std::map< MVertex *, double > _angles
Definition: BackgroundMesh.h:54
GFace::edges
virtual std::vector< GEdge * > const & edges() const
Definition: GFace.h:121
backgroundMesh::propagateCrossFieldByDistance
void propagateCrossFieldByDistance(GFace *)
Definition: BackgroundMesh.cpp:291
Field.h
GFace::firstDer
virtual Pair< SVector3, SVector3 > firstDer(const SPoint2 &param) const =0
constantPerElement
Definition: simpleFunction.h:42
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
backgroundMesh::print
void print(const std::string &filename, GFace *gf, const std::map< MVertex *, double > &, int smooth=0)
Definition: BackgroundMesh.cpp:672
GPoint::y
double y() const
Definition: GPoint.h:22
GFace.h
backgroundMesh::operator()
double operator()(double u, double v, double w) const
Definition: BackgroundMesh.cpp:568
linearSystemPETSc
Definition: linearSystemPETSc.h:150
GModel::getFaces
std::set< GFace *, GEntityPtrLessThan > getFaces() const
Definition: GModel.h:376
GFace
Definition: GFace.h:33
GFace::triangles
std::vector< MTriangle * > triangles
Definition: GFace.h:428
backgroundMesh::propagate1dMesh
void propagate1dMesh(GFace *)
Definition: BackgroundMesh.cpp:237
backgroundMesh::getSmoothness
double getSmoothness(double u, double v, double w)
Definition: BackgroundMesh.cpp:389
angle
double angle(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:157
SPoint2
Definition: SPoint2.h:12
BackgroundMeshGEdge::ge
GEdge * ge
Definition: BackgroundMesh.h:130
OS.h
Msg::Debug
static void Debug(const char *fmt,...)
Definition: GmshMessage.cpp:752
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
BackgroundMeshGEdge::lines
std::vector< MLine > lines
Definition: BackgroundMesh.h:131
MVertex
Definition: MVertex.h:24
MElementOctree.h
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
backgroundMesh::_2Dto3D
std::map< MVertex *, MVertex * > _2Dto3D
Definition: BackgroundMesh.h:52
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
BackgroundMesh.h
reparamMeshVertexOnFace
bool reparamMeshVertexOnFace(MVertex const *v, const GFace *gf, SPoint2 &param, bool onSurface, bool failOnSeam, int dir)
Definition: MVertex.cpp:522
GlobalBackgroundMesh::~GlobalBackgroundMesh
~GlobalBackgroundMesh()
Definition: BackgroundMesh.cpp:751
backgroundMesh::_octree
MElementOctree * _octree
Definition: BackgroundMesh.h:47
backgroudMeshExists
bool backgroudMeshExists(const std::string &name)
Definition: BackgroundMesh.cpp:733
linearSystemPETSc.h
SPoint3
Definition: SPoint3.h:14
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
laplaceTerm
Definition: laplaceTerm.h:12
GFace::point
virtual GPoint point(double par1, double par2) const =0
contextMeshOptions::lcMin
double lcMin
Definition: Context.h:25
BackgroundMeshGFace::gf
GFace * gf
Definition: BackgroundMesh.h:135
SVector3
Definition: SVector3.h:16
backgroundMesh::propagateCrossField
void propagateCrossField(GFace *, simpleFunction< double > *)
Definition: BackgroundMesh.cpp:448
GlobalBackgroundMesh::edgeBackgroundMeshes
std::unordered_map< GEdge *, BackgroundMeshGEdge > edgeBackgroundMeshes
Definition: BackgroundMesh.h:148
MTriangle::getVertex
virtual MVertex * getVertex(int num)
Definition: MTriangle.h:62
linearSystemCSR.h
MVertex::onWhat
GEntity * onWhat() const
Definition: MVertex.h:82
laplaceTerm.h
backgroundMesh::inDomain
bool inDomain(double u, double v, double w) const
Definition: BackgroundMesh.cpp:562
GmshMessage.h
MLine.h
edges
static int edges[6][2]
Definition: meshGRegionLocalMeshMod.cpp:23
global_bmeshes
std::vector< std::unique_ptr< GlobalBackgroundMesh > > global_bmeshes
Definition: BackgroundMesh.cpp:731
femTerm::addToMatrix
void addToMatrix(dofManager< dataVec > &dm, groupOfElements &L, groupOfElements &C) const
Definition: femTerm.h:53
GPoint
Definition: GPoint.h:13
dofManager::sizeOfR
virtual int sizeOfR() const
Definition: dofManager.h:606
MAX_THREADS
static const int MAX_THREADS
Definition: BackgroundMesh.cpp:35
reparamMeshVertexOnEdge
bool reparamMeshVertexOnEdge(MVertex *v, const GEdge *ge, double &param)
Definition: MVertex.cpp:592
GFace::quadrangles
std::vector< MQuadrangle * > quadrangles
Definition: GFace.h:429
contextMeshOptions::lcMax
double lcMax
Definition: Context.h:25
Fopen
FILE * Fopen(const char *f, const char *mode)
Definition: OS.cpp:273
backgroundMesh::setCrossFieldsByDistance
static void setCrossFieldsByDistance(GFace *)
Definition: BackgroundMesh.cpp:52
backgroundMesh::backgroundMesh
backgroundMesh(GFace *, bool dist=false)
Definition: BackgroundMesh.cpp:75
MLine
Definition: MLine.h:21
Msg::GetThreadNum
static int GetThreadNum()
Definition: GmshMessage.cpp:1638
GPoint::z
double z() const
Definition: GPoint.h:23
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
Pair::first
L first() const
Definition: Pair.h:22
MElement::getVertex
virtual const MVertex * getVertex(int num) const =0
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
GEdge.h
GEdge::firstDer
virtual SVector3 firstDer(double par) const =0
MVertex::setParameter
virtual bool setParameter(int i, double par)
Definition: MVertex.h:102
GModel::getVertices
std::set< GVertex *, GEntityPtrLessThan > getVertices() const
Definition: GModel.h:378
crossField2d::crossField2d
crossField2d(MVertex *, GEdge *)
Definition: BackgroundMesh.cpp:276
backgroundMesh::propagateCrossFieldHJ
void propagateCrossFieldHJ(GFace *)
Definition: BackgroundMesh.cpp:442
dofManager::numberVertex
void numberVertex(MVertex *v, int iComp, int iField)
Definition: dofManager.h:231
Pair::second
R second() const
Definition: Pair.h:23
MElement::interpolateGrad
void interpolateGrad(double val[], double u, double v, double w, double f[], int stride=1, double invjac[3][3]=nullptr, int order=-1)
Definition: MElement.cpp:1230
simpleFunction< double >
MEdgeVertex
Definition: MVertex.h:137
GVertex
Definition: GVertex.h:23
MVertex.h
linearSystemFull::systemSolve
virtual int systemSolve()
Definition: linearSystemFull.h:77
Numeric.h
GModel
Definition: GModel.h:44
reparamMeshEdgeOnFace
bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf, SPoint2 &param1, SPoint2 &param2)
Definition: MVertex.cpp:474
BGM_MeshSize
double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double Z)
Definition: BackgroundMeshTools.cpp:255
MElementOctree
Definition: MElementOctree.h:15
SElement
Definition: SElement.h:18
signedDistancePointLine
void signedDistancePointLine(const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p, double &d, SPoint3 &closePt)
Definition: Numeric.cpp:908
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
CTX::mesh
contextMeshOptions mesh
Definition: Context.h:313
GlobalBackgroundMesh
Definition: BackgroundMesh.h:144
dofManager< double >
MElement
Definition: MElement.h:30
crossField2d::normalizeAngle
static void normalizeAngle(double &angle)
Definition: BackgroundMesh.h:34
dofManager.h
GEntity::tag
int tag() const
Definition: GEntity.h:280
MFaceVertex
Definition: MVertex.h:168
linearSystemFull.h
SVector3::x
double x(void) const
Definition: SVector3.h:30
linearSystemCSRGmm
Definition: linearSystemCSR.h:176
GModel::getEdges
std::set< GEdge *, GEntityPtrLessThan > getEdges() const
Definition: GModel.h:377
MElementOctree::find
MElement * find(double x, double y, double z, int dim=-1, bool strict=false) const
Definition: MElementOctree.cpp:205
BackgroundMeshGFace
Definition: BackgroundMesh.h:134
GlobalBackgroundMesh::gm
GModel * gm
Definition: BackgroundMesh.h:147
BackgroundMeshGEdge
Definition: BackgroundMesh.h:129
MTriangle
Definition: MTriangle.h:26
dofManager::getDofValue
virtual void getDofValue(std::vector< Dof > &keys, std::vector< dataVec > &Vals)
Definition: dofManager.h:235
GlobalBackgroundMesh::importGModelMeshes
int importGModelMeshes(GModel *gm, bool overwriteExisting=true)
Fill the entityMesh map by copying the meshes in the GModel. New MVertex, MLine and MTriangle instanc...
Definition: BackgroundMesh.cpp:764
MVertex::getParameter
virtual bool getParameter(int i, double &par) const
Definition: MVertex.h:97
SPoint2::x
double x(void) const
Definition: SPoint2.h:86
Context.h
backgroundMesh::unset
static void unset()
Definition: BackgroundMesh.cpp:60
backgroundMesh::_vertices
std::vector< MVertex * > _vertices
Definition: BackgroundMesh.h:48
SVector3::y
double y(void) const
Definition: SVector3.h:31
GVertex.h
backgroundMesh::~backgroundMesh
~backgroundMesh()
Definition: BackgroundMesh.cpp:154
z
const double z
Definition: GaussQuadratureQuad.cpp:56
MQuadrangle.h
crossField2d::_angle
double _angle
Definition: BackgroundMesh.h:33
MElement.h
Pair
Definition: Pair.h:10
backgroundMesh::updateSizes
void updateSizes(GFace *)
Definition: BackgroundMesh.cpp:505
GEdge
Definition: GEdge.h:26
BackgroundMeshGFace::triangles
std::vector< MTriangle > triangles
Definition: BackgroundMesh.h:136
MElement::maxEdge
virtual double maxEdge()
Definition: MElement.cpp:246
backgroundMesh::current
static backgroundMesh * current()
Definition: BackgroundMesh.cpp:68
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
GModel.h
constantPerElement::set
void set(MElement *e, scalar v)
Definition: simpleFunction.h:48
MVertex::y
double y() const
Definition: MVertex.h:61
backgroundMesh::_3Dto2D
std::map< MVertex *, MVertex * > _3Dto2D
Definition: BackgroundMesh.h:51
backgroundMesh::getAngle
double getAngle(double u, double v, double w) const
Definition: BackgroundMesh.cpp:603
GlobalBackgroundMesh::faceBackgroundMeshes
std::unordered_map< GFace *, BackgroundMeshGFace > faceBackgroundMeshes
Definition: BackgroundMesh.h:149
GEntity::dim
virtual int dim() const
Definition: GEntity.h:196
GPoint::x
double x() const
Definition: GPoint.h:21
dofManager::fixVertex
void fixVertex(MVertex const *v, int iComp, int iField, const dataVec &value)
Definition: dofManager.h:170
GlobalBackgroundMesh::mesh_vertices
std::vector< MVertex * > mesh_vertices
Definition: BackgroundMesh.h:150
backgroundMesh::_current
static std::vector< backgroundMesh * > _current
Definition: BackgroundMesh.h:55
backgroundMesh
Definition: BackgroundMesh.h:46
myAngle
double myAngle(const SVector3 &a, const SVector3 &b, const SVector3 &d)
Definition: BackgroundMesh.cpp:358
backgroundMesh::_triangles
std::vector< MElement * > _triangles
Definition: BackgroundMesh.h:49
MVertex::x
double x() const
Definition: MVertex.h:60
getBackgroundMesh
GlobalBackgroundMesh & getBackgroundMesh(const std::string &name)
Definition: BackgroundMesh.cpp:741
SVector3::normalize
double normalize()
Definition: SVector3.h:38
propagateValuesOnFace
static void propagateValuesOnFace(GFace *_gf, std::map< MVertex *, double > &dirichlet, simpleFunction< double > *ONE, bool in_parametric_plane=false)
Definition: BackgroundMesh.cpp:169
SPoint2::y
double y(void) const
Definition: SPoint2.h:88
backgroundMesh::set
static void set(GFace *)
Definition: BackgroundMesh.cpp:40