gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
meshGFace.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 "meshGFace.h"
7 #include "BDS.h"
8 #include "BackgroundMesh.h"
9 #include "Context.h"
10 #include "DivideAndConquer.h"
11 #include "GEdge.h"
12 #include "GFace.h"
13 #include "GModel.h"
14 #include "GPoint.h"
15 #include "GVertex.h"
16 #include "GmshMessage.h"
17 #include "HighOrder.h"
18 #include "MElementOctree.h"
19 #include "MLine.h"
20 #include "MQuadrangle.h"
21 #include "MTriangle.h"
22 #include "MVertex.h"
23 #include "Numeric.h"
24 #include "OS.h"
25 #include "boundaryLayersData.h"
26 #include "discreteEdge.h"
27 #include "discreteFace.h"
28 #include "filterElements.h"
29 #include "meshGEdge.h"
30 #include "meshGFaceBDS.h"
31 #include "meshGFaceBamg.h"
32 #include "meshGFaceBipartiteLabelling.h"
33 #include "meshGFaceDelaunayInsertion.h"
34 #include "meshGFaceOptimize.h"
35 #include "meshTriangulation.h"
36 #include "qualityMeasures.h"
37 #include "robustPredicates.h"
38 #include <limits>
39 #include <map>
40 #include <sstream>
41 #include <stdlib.h>
42 
43 bool pointInsideParametricDomain(std::vector<SPoint2> &bnd, SPoint2 &p, SPoint2 &out, int &N)
44 {
45  int count = 0;
46  for (size_t i = 0; i < bnd.size(); i += 2)
47  {
48  SPoint2 p1 = bnd[i];
49  SPoint2 p2 = bnd[i + 1];
50  double a = robustPredicates::orient2d(p1, p2, p);
51  double b = robustPredicates::orient2d(p1, p2, out);
52  if (a * b < 0)
53  {
54  a = robustPredicates::orient2d(p, out, p1);
55  b = robustPredicates::orient2d(p, out, p2);
56  if (a * b < 0)
57  count++;
58  }
59  }
60  N = count;
61  if (count % 2 == 0)
62  return false;
63  return true;
64 }
65 
66 static void trueBoundary(GFace *gf, std::vector<SPoint2> &bnd, int debug)
67 {
68  FILE *view_t = nullptr;
69  if (debug)
70  {
71  char name[245];
72  sprintf(name, "trueBoundary%d.pos", gf->tag());
73  view_t = Fopen(name, "w");
74  if (view_t)
75  fprintf(view_t, "View \"True Boundary\"{\n");
76  }
77  std::vector<GEdge *> edg = gf->edges();
78  std::set<GEdge *> edges(edg.begin(), edg.end());
79 
80  for (auto it = edges.begin(); it != edges.end(); ++it)
81  {
82  GEdge *ge = *it;
83  Range<double> r = ge->parBoundsOnFace(gf);
84  SPoint2 p[300];
85  int NITER = ge->isSeam(gf) ? 2 : 1;
86  for (int i = 0; i < NITER; i++)
87  {
88  int count = NITER == 2 ? 300 : 300;
89  for (int k = 0; k < count; k++)
90  {
91  double t = (double)k / (count - 1);
92  double xi = r.low() + (r.high() - r.low()) * t;
93  p[k] = ge->reparamOnFace(gf, xi, i);
94  if (k > 0)
95  {
96  if (view_t)
97  {
98  fprintf(view_t, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n", p[k - 1].x(), p[k - 1].y(), 0.0, p[k].x(),
99  p[k].y(), 0.0, ge->tag(), ge->tag());
100  }
101  bnd.push_back(p[k - 1]);
102  bnd.push_back(p[k]);
103  }
104  }
105  }
106  }
107  if (view_t)
108  {
109  fprintf(view_t, "};\n");
110  fclose(view_t);
111  }
112 }
113 
114 static void computeElementShapes(GFace *gf, double &worst, double &avg, double &best, int &nT, int &greaterThan)
115 {
116  worst = 1.e22;
117  avg = 0.0;
118  best = 0.0;
119  nT = 0;
120  greaterThan = 0;
121  for (std::size_t i = 0; i < gf->triangles.size(); i++)
122  {
123  double q = qmTriangle::gamma(gf->triangles[i]);
124  if (q > .9)
125  greaterThan++;
126  avg += q;
127  worst = std::min(worst, q);
128  best = std::max(best, q);
129  nT++;
130  }
131  avg /= nT;
132 }
133 
134 class quadMeshRemoveHalfOfOneDMesh
135 {
136  private:
137  GFace *_gf;
138  std::map<GEdge *, std::vector<MLine *>> _backup;
139  std::map<MEdge, MVertex *, MEdgeLessThan> _middle;
140  void _subdivide()
141  {
142  std::vector<MQuadrangle *> qnew;
143  std::map<MEdge, MVertex *, MEdgeLessThan> eds;
144  for (std::size_t i = 0; i < _gf->triangles.size(); i++)
145  {
146  MVertex *v[3];
147  SPoint2 m[3];
148  for (int j = 0; j < 3; j++)
149  {
150  MEdge E = _gf->triangles[i]->getEdge(j);
151  SPoint2 p1, p2;
152  reparamMeshEdgeOnFace(E.getVertex(0), E.getVertex(1), _gf, p1, p2);
153  auto it = _middle.find(E);
154  auto it2 = eds.find(E);
155  m[j] = p1;
156  if (it == _middle.end() && it2 == eds.end())
157  {
158  GPoint gp = _gf->point((p1 + p2) * 0.5);
159  double XX = 0.5 * (E.getVertex(0)->x() + E.getVertex(1)->x());
160  double YY = 0.5 * (E.getVertex(0)->y() + E.getVertex(1)->y());
161  double ZZ = 0.5 * (E.getVertex(0)->z() + E.getVertex(1)->z());
162  v[j] = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
163  _gf->mesh_vertices.push_back(v[j]);
164  eds[E] = v[j];
165  }
166  else if (it == _middle.end())
167  {
168  v[j] = it2->second;
169  }
170  else
171  {
172  v[j] = it->second;
173  v[j]->onWhat()->mesh_vertices.push_back(v[j]);
174  if (!CTX::instance()->mesh.secondOrderLinear)
175  {
176  // re-push middle vertex on the curve (this can of course lead to an
177  // invalid mesh)
178  double u = 0.;
179  if (v[j]->getParameter(0, u) && v[j]->onWhat()->dim() == 1)
180  {
181  GEdge *ge = static_cast<GEdge *>(v[j]->onWhat());
182  GPoint p = ge->point(u);
183  v[j]->x() = p.x();
184  v[j]->y() = p.y();
185  v[j]->z() = p.z();
186  }
187  }
188  }
189  }
190  GPoint gp = _gf->point((m[0] + m[1] + m[2]) * (1. / 3.));
191  double XX = (v[0]->x() + v[1]->x() + v[2]->x()) * (1. / 3.);
192  double YY = (v[0]->y() + v[1]->y() + v[2]->y()) * (1. / 3.);
193  double ZZ = (v[0]->z() + v[1]->z() + v[2]->z()) * (1. / 3.);
194  MFaceVertex *vmid = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
195  _gf->mesh_vertices.push_back(vmid);
196  qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(0), v[0], vmid, v[2]));
197  qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(1), v[1], vmid, v[0]));
198  qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(2), v[2], vmid, v[1]));
199  delete _gf->triangles[i];
200  }
201  _gf->triangles.clear();
202  for (std::size_t i = 0; i < _gf->quadrangles.size(); i++)
203  {
204  MVertex *v[4];
205  SPoint2 m[4];
206  for (int j = 0; j < 4; j++)
207  {
208  MEdge E = _gf->quadrangles[i]->getEdge(j);
209  SPoint2 p1, p2;
210  reparamMeshEdgeOnFace(E.getVertex(0), E.getVertex(1), _gf, p1, p2);
211  auto it = _middle.find(E);
212  auto it2 = eds.find(E);
213  m[j] = p1;
214  if (it == _middle.end() && it2 == eds.end())
215  {
216  GPoint gp = _gf->point((p1 + p2) * 0.5);
217  double XX = 0.5 * (E.getVertex(0)->x() + E.getVertex(1)->x());
218  double YY = 0.5 * (E.getVertex(0)->y() + E.getVertex(1)->y());
219  double ZZ = 0.5 * (E.getVertex(0)->z() + E.getVertex(1)->z());
220  v[j] = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
221  _gf->mesh_vertices.push_back(v[j]);
222  eds[E] = v[j];
223  }
224  else if (it == _middle.end())
225  {
226  v[j] = it2->second;
227  }
228  else
229  {
230  v[j] = it->second;
231  v[j]->onWhat()->mesh_vertices.push_back(v[j]);
232  if (!CTX::instance()->mesh.secondOrderLinear)
233  {
234  // re-push middle vertex on the curve (this can of course lead to an
235  // invalid mesh)
236  double u = 0.;
237  if (v[j]->getParameter(0, u) && v[j]->onWhat()->dim() == 1)
238  {
239  GEdge *ge = static_cast<GEdge *>(v[j]->onWhat());
240  GPoint p = ge->point(u);
241  v[j]->x() = p.x();
242  v[j]->y() = p.y();
243  v[j]->z() = p.z();
244  }
245  }
246  }
247  }
248  GPoint gp = _gf->point((m[0] + m[1] + m[2] + m[3]) * 0.25);
249  // FIXME: not exactly correct, but that's the place where we want the
250  // point to reside
251  double XX = 0.25 * (v[0]->x() + v[1]->x() + v[2]->x() + v[3]->x());
252  double YY = 0.25 * (v[0]->y() + v[1]->y() + v[2]->y() + v[3]->y());
253  double ZZ = 0.25 * (v[0]->z() + v[1]->z() + v[2]->z() + v[3]->z());
254  MVertex *vmid = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
255  _gf->mesh_vertices.push_back(vmid);
256  qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(0), v[0], vmid, v[3]));
257  qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(1), v[1], vmid, v[0]));
258  qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(2), v[2], vmid, v[1]));
259  qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(3), v[3], vmid, v[2]));
260  delete _gf->quadrangles[i];
261  }
262  _gf->quadrangles = qnew;
263  }
264  void _restore()
265  {
266  std::vector<GEdge *> const &edges = _gf->edges();
267  auto ite = edges.begin();
268  while (ite != edges.end())
269  {
270  for (std::size_t i = 0; i < (*ite)->lines.size(); i++)
271  {
272  delete (*ite)->lines[i];
273  }
274  (*ite)->lines = _backup[*ite];
275  ++ite;
276  }
277  }
278 
279  public:
280  // remove one point every two and remember middle points
281  quadMeshRemoveHalfOfOneDMesh(GFace *gf, bool periodic) : _gf(gf)
282  {
283  // only do it if a full recombination has to (and can) be done
284  if (!CTX::instance()->mesh.recombineAll && !gf->meshAttributes.recombine)
285  return;
286  if (CTX::instance()->mesh.algoRecombine < 2 && CTX::instance()->mesh.algoRecombine != 4)
287  return;
288  if (gf->compound.size())
289  return;
290  if (periodic)
291  {
292  Msg::Error("Full-quad recombination not ready yet for periodic surfaces");
293  return;
294  }
295  std::vector<GEdge *> const &edges = gf->edges();
296  auto ite = edges.begin();
297  while (ite != edges.end())
298  {
299  if ((*ite)->meshAttributes.method == MESH_TRANSFINITE)
300  {
301  Msg::Warning("Full-quad recombination only compatible with "
302  "transfinite meshes if those are performed first");
303  }
304  if (!(*ite)->isMeshDegenerated())
305  {
306  std::vector<MLine *> temp;
307  (*ite)->mesh_vertices.clear();
308  for (std::size_t i = 0; i < (*ite)->lines.size(); i += 2)
309  {
310  if (i + 1 >= (*ite)->lines.size())
311  {
312  Msg::Error("1D mesh cannot be divided by 2");
313  break;
314  }
315  MVertex *v1 = (*ite)->lines[i]->getVertex(0);
316  MVertex *v2 = (*ite)->lines[i]->getVertex(1);
317  MVertex *v3 = (*ite)->lines[i + 1]->getVertex(1);
318  v2->x() = 0.5 * (v1->x() + v3->x());
319  v2->y() = 0.5 * (v1->y() + v3->y());
320  v2->z() = 0.5 * (v1->z() + v3->z());
321  temp.push_back(new MLine(v1, v3));
322  if (v1->onWhat() == *ite)
323  (*ite)->mesh_vertices.push_back(v1);
324  _middle[MEdge(v1, v3)] = v2;
325  }
326  _backup[*ite] = (*ite)->lines;
327  (*ite)->lines = temp;
328  }
329  ++ite;
330  }
331  CTX::instance()->mesh.lcFactor *= 2.0;
332  }
333  void finish()
334  {
335  if (_backup.empty())
336  return;
337  // recombine the elements on the half mesh
338  CTX::instance()->mesh.lcFactor /= 2.0;
339  bool blossom = (CTX::instance()->mesh.algoRecombine == 3);
340  int topo = CTX::instance()->mesh.recombineOptimizeTopology;
341  recombineIntoQuads(_gf, blossom, topo, true, 0.1);
342  _subdivide();
343  _restore();
344  recombineIntoQuads(_gf, blossom, topo, true, 1.e-3);
345  computeElementShapes(_gf, _gf->meshStatistics.worst_element_shape, _gf->meshStatistics.average_element_shape,
346  _gf->meshStatistics.best_element_shape, _gf->meshStatistics.nbTriangle,
347  _gf->meshStatistics.nbGoodQuality);
348  }
349 };
350 
351 static void copyMesh(GFace *source, GFace *target)
352 {
353  std::map<MVertex *, MVertex *> vs2vt;
354 
355  // add principal GVertex pairs
356 
357  std::vector<GVertex *> s_vtcs = source->vertices();
358  s_vtcs.insert(s_vtcs.end(), source->embeddedVertices().begin(), source->embeddedVertices().end());
359  for (auto it = source->embeddedEdges().begin(); it != source->embeddedEdges().end(); it++)
360  {
361  if ((*it)->getBeginVertex())
362  s_vtcs.push_back((*it)->getBeginVertex());
363  if ((*it)->getEndVertex())
364  s_vtcs.push_back((*it)->getEndVertex());
365  }
366  std::vector<GVertex *> t_vtcs = target->vertices();
367  t_vtcs.insert(t_vtcs.end(), target->embeddedVertices().begin(), target->embeddedVertices().end());
368  for (auto it = target->embeddedEdges().begin(); it != target->embeddedEdges().end(); it++)
369  {
370  if ((*it)->getBeginVertex())
371  t_vtcs.push_back((*it)->getBeginVertex());
372  if ((*it)->getEndVertex())
373  t_vtcs.push_back((*it)->getEndVertex());
374  }
375 
376  if (s_vtcs.size() != t_vtcs.size())
377  {
378  Msg::Info("Periodicity imposed on topologically incompatible surfaces"
379  "(%d vs %d points)",
380  s_vtcs.size(), t_vtcs.size());
381  }
382 
383  std::set<GVertex *> checkVtcs(s_vtcs.begin(), s_vtcs.end());
384 
385  for (auto tvIter = t_vtcs.begin(); tvIter != t_vtcs.end(); ++tvIter)
386  {
387  GVertex *gvt = *tvIter;
388  auto gvsIter = target->vertexCounterparts.find(gvt);
389 
390  if (gvsIter == target->vertexCounterparts.end())
391  {
392  Msg::Error("Periodic meshing of surface %d with surface %d: "
393  "point %d has no periodic counterpart",
394  target->tag(), source->tag(), gvt->tag());
395  }
396  else
397  {
398  GVertex *gvs = gvsIter->second;
399  if (checkVtcs.find(gvs) == checkVtcs.end())
400  {
401  if (gvs)
402  Msg::Error("Periodic meshing of surface %d with surface %d: "
403  "point %d has periodic counterpart %d outside of source surface",
404  target->tag(), source->tag(), gvt->tag(), gvs->tag());
405 
406  else
407  Msg::Error("Periodic meshing of surface %d with surface %d: "
408  "point %d has no periodic counterpart",
409  target->tag(), source->tag(), gvt->tag());
410  }
411  if (gvs)
412  {
413  MVertex *vs = gvs->mesh_vertices[0];
414  MVertex *vt = gvt->mesh_vertices[0];
415  vs2vt[vs] = vt;
416  target->correspondingVertices[vt] = vs;
417  }
418  }
419  }
420 
421  // add corresponding curve nodes assuming curves were correctly meshed already
422 
423  std::vector<GEdge *> s_edges = source->edges();
424  s_edges.insert(s_edges.end(), source->embeddedEdges().begin(), source->embeddedEdges().end());
425  std::vector<GEdge *> t_edges = target->edges();
426  t_edges.insert(t_edges.end(), target->embeddedEdges().begin(), target->embeddedEdges().end());
427 
428  std::set<GEdge *> checkEdges;
429  checkEdges.insert(s_edges.begin(), s_edges.end());
430 
431  for (auto te_iter = t_edges.begin(); te_iter != t_edges.end(); ++te_iter)
432  {
433  GEdge *get = *te_iter;
434 
435  auto gesIter = target->edgeCounterparts.find(get);
436  if (gesIter == target->edgeCounterparts.end())
437  {
438  Msg::Error("Periodic meshing of surface %d with surface %d: "
439  "curve %d has no periodic counterpart",
440  target->tag(), source->tag(), get->tag());
441  }
442  else
443  {
444  GEdge *ges = gesIter->second.first;
445  if (checkEdges.find(ges) == checkEdges.end())
446  {
447  Msg::Error("Periodic meshing of surface %d with surface %d: "
448  "curve %d has periodic counterpart %d",
449  target->tag(), source->tag(), get->tag(), ges->tag());
450  }
451  if (get->mesh_vertices.size() != ges->mesh_vertices.size())
452  {
453  Msg::Error("Periodic meshing of surface %d with surface %d: "
454  "curve %d has %d vertices, whereas correspondant %d has %d",
455  target->tag(), source->tag(), get->tag(), get->mesh_vertices.size(), ges->tag(),
456  ges->mesh_vertices.size());
457  }
458  int orientation = gesIter->second.second;
459  int is = orientation == 1 ? 0 : get->mesh_vertices.size() - 1;
460  for (unsigned it = 0; it < get->mesh_vertices.size(); it++, is += orientation)
461  {
462  MVertex *vs = ges->mesh_vertices[is];
463  MVertex *vt = get->mesh_vertices[it];
464  vs2vt[vs] = vt;
465  target->correspondingVertices[vt] = vs;
466  }
467  }
468  }
469 
470  // transform interior nodes
471  std::vector<double> &tfo = target->affineTransform;
472 
473  for (std::size_t i = 0; i < source->mesh_vertices.size(); i++)
474  {
475  MVertex *vs = source->mesh_vertices[i];
476  SPoint2 XXX;
477 
478  double ps[4] = {vs->x(), vs->y(), vs->z(), 1.};
479  double res[4] = {0., 0., 0., 0.};
480  int idx = 0;
481  for (int i = 0; i < 4; i++)
482  for (int j = 0; j < 4; j++)
483  res[i] += tfo[idx++] * ps[j];
484 
485  SPoint3 tp(res[0], res[1], res[2]);
486  XXX = target->parFromPoint(tp);
487 
488  GPoint gp = target->point(XXX);
489  MVertex *vt = new MFaceVertex(gp.x(), gp.y(), gp.z(), target, gp.u(), gp.v());
490  target->mesh_vertices.push_back(vt);
491  target->correspondingVertices[vt] = vs;
492  vs2vt[vs] = vt;
493  }
494 
495  // create new elements
496  for (unsigned i = 0; i < source->triangles.size(); i++)
497  {
498  MVertex *v1 = vs2vt[source->triangles[i]->getVertex(0)];
499  MVertex *v2 = vs2vt[source->triangles[i]->getVertex(1)];
500  MVertex *v3 = vs2vt[source->triangles[i]->getVertex(2)];
501  if (v1 && v2 && v3)
502  {
503  target->triangles.push_back(new MTriangle(v1, v2, v3));
504  }
505  else
506  {
507  Msg::Error("Could not find periodic counterpart of triangle nodes "
508  "%lu %lu %lu",
509  source->triangles[i]->getVertex(0)->getNum(), source->triangles[i]->getVertex(1)->getNum(),
510  source->triangles[i]->getVertex(2)->getNum());
511  }
512  }
513 
514  for (unsigned i = 0; i < source->quadrangles.size(); i++)
515  {
516  MVertex *v1 = vs2vt[source->quadrangles[i]->getVertex(0)];
517  MVertex *v2 = vs2vt[source->quadrangles[i]->getVertex(1)];
518  MVertex *v3 = vs2vt[source->quadrangles[i]->getVertex(2)];
519  MVertex *v4 = vs2vt[source->quadrangles[i]->getVertex(3)];
520  if (v1 && v2 && v3 && v4)
521  {
522  target->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
523  }
524  else
525  {
526  Msg::Error("Could not find periodic counterpart of quadrangle nodes "
527  "%lu %lu %lu %lu",
528  source->quadrangles[i]->getVertex(0)->getNum(), source->quadrangles[i]->getVertex(1)->getNum(),
529  source->quadrangles[i]->getVertex(2)->getNum(), source->quadrangles[i]->getVertex(3)->getNum());
530  }
531  }
532 }
533 
534 static void remeshUnrecoveredEdges(std::multimap<MVertex *, BDS_Point *> &recoverMultiMapInv,
535  std::set<EdgeToRecover> &edgesNotRecovered, bool all)
536 {
537  deMeshGFace dem;
538 
539  auto itr = edgesNotRecovered.begin();
540  for (; itr != edgesNotRecovered.end(); ++itr)
541  {
542  std::vector<GFace *> l_faces = itr->ge->faces();
543  // un-mesh model faces adjacent to the model edge
544  for (auto it = l_faces.begin(); it != l_faces.end(); ++it)
545  {
546  if ((*it)->triangles.size() || (*it)->quadrangles.size())
547  {
548  (*it)->meshStatistics.status = GFace::PENDING;
549  dem(*it);
550  }
551  }
552 
553  // add a new point in the middle of the intersecting segment
554  int p1 = itr->p1;
555  int p2 = itr->p2;
556  int N = itr->ge->lines.size();
557  GVertex *g1 = itr->ge->getBeginVertex();
558  GVertex *g2 = itr->ge->getEndVertex();
559  Range<double> bb = itr->ge->parBounds(0);
560 
561  std::vector<MLine *> newLines;
562 
563  for (int i = 0; i < N; i++)
564  {
565  MVertex *v1 = itr->ge->lines[i]->getVertex(0);
566  MVertex *v2 = itr->ge->lines[i]->getVertex(1);
567 
568  auto itp1 = recoverMultiMapInv.lower_bound(v1);
569  auto itp2 = recoverMultiMapInv.lower_bound(v2);
570 
571  if (itp1 != recoverMultiMapInv.end() && itp2 != recoverMultiMapInv.end())
572  {
573  BDS_Point *pp1 = itp1->second;
574  BDS_Point *pp2 = itp2->second;
575  BDS_Point *pp1b = itp1->second;
576  BDS_Point *pp2b = itp2->second;
577  if (recoverMultiMapInv.count(v1) == 2)
578  {
579  itp1++;
580  pp1b = itp1->second;
581  }
582  if (recoverMultiMapInv.count(v2) == 2)
583  {
584  itp2++;
585  pp2b = itp2->second;
586  }
587 
588  if ((pp1->iD == p1 && pp2->iD == p2) || (pp1->iD == p2 && pp2->iD == p1) ||
589  (pp1b->iD == p1 && pp2b->iD == p2) || (pp1b->iD == p2 && pp2b->iD == p1) || all)
590  {
591  double t1;
592  double lc1 = -1;
593  if (v1->onWhat() == g1)
594  t1 = bb.low();
595  else if (v1->onWhat() == g2)
596  t1 = bb.high();
597  else
598  {
599  MEdgeVertex *ev1 = (MEdgeVertex *)v1;
600  lc1 = ev1->getLc();
601  v1->getParameter(0, t1);
602  }
603  double t2;
604  double lc2 = -1;
605  if (v2->onWhat() == g1)
606  t2 = bb.low();
607  else if (v2->onWhat() == g2)
608  t2 = bb.high();
609  else
610  {
611  MEdgeVertex *ev2 = (MEdgeVertex *)v2;
612  lc2 = ev2->getLc();
613  v2->getParameter(0, t2);
614  }
615 
616  // periodic
617  if (v1->onWhat() == g1 && v1->onWhat() == g2)
618  t1 = fabs(t2 - bb.low()) < fabs(t2 - bb.high()) ? bb.low() : bb.high();
619  if (v2->onWhat() == g1 && v2->onWhat() == g2)
620  t2 = fabs(t1 - bb.low()) < fabs(t1 - bb.high()) ? bb.low() : bb.high();
621 
622  if (lc1 == -1)
623  lc1 = BGM_MeshSize(v1->onWhat(), 0, 0, v1->x(), v1->y(), v1->z());
624  if (lc2 == -1)
625  lc2 = BGM_MeshSize(v2->onWhat(), 0, 0, v2->x(), v2->y(), v2->z());
626  // should be better, i.e. equidistant
627  double t = 0.5 * (t2 + t1);
628  double lc = 0.5 * (lc1 + lc2);
629  GPoint V = itr->ge->point(t);
630  MEdgeVertex *newv = new MEdgeVertex(V.x(), V.y(), V.z(), itr->ge, t, 0, lc);
631  newLines.push_back(new MLine(v1, newv));
632  newLines.push_back(new MLine(newv, v2));
633  delete itr->ge->lines[i];
634  }
635  else
636  {
637  newLines.push_back(itr->ge->lines[i]);
638  }
639  }
640  else
641  {
642  newLines.push_back(itr->ge->lines[i]);
643  }
644  }
645 
646  itr->ge->lines = newLines;
647  itr->ge->mesh_vertices.clear();
648  N = itr->ge->lines.size();
649  for (int i = 1; i < N; i++)
650  {
651  itr->ge->mesh_vertices.push_back(itr->ge->lines[i]->getVertex(0));
652  }
653  }
654 }
655 
656 static void remeshUnrecoveredEdges(std::map<MVertex *, BDS_Point *> &recoverMapInv,
657  std::set<EdgeToRecover> &edgesNotRecovered)
658 {
659  deMeshGFace dem;
660 
661  auto itr = edgesNotRecovered.begin();
662  for (; itr != edgesNotRecovered.end(); ++itr)
663  {
664  std::vector<GFace *> l_faces = itr->ge->faces();
665  // un-mesh model faces adjacent to the model edge
666  for (auto it = l_faces.begin(); it != l_faces.end(); ++it)
667  {
668  if ((*it)->triangles.size() || (*it)->quadrangles.size())
669  {
670  (*it)->meshStatistics.status = GFace::PENDING;
671  dem(*it);
672  }
673  }
674 
675  // add a new point in the middle of the intersecting segment
676  int p1 = itr->p1;
677  int p2 = itr->p2;
678  int N = itr->ge->lines.size();
679  GVertex *g1 = itr->ge->getBeginVertex();
680  GVertex *g2 = itr->ge->getEndVertex();
681  Range<double> bb = itr->ge->parBounds(0);
682 
683  std::vector<MLine *> newLines;
684 
685  for (int i = 0; i < N; i++)
686  {
687  MVertex *v1 = itr->ge->lines[i]->getVertex(0);
688  MVertex *v2 = itr->ge->lines[i]->getVertex(1);
689  auto itp1 = recoverMapInv.find(v1);
690  auto itp2 = recoverMapInv.find(v2);
691  if (itp1 != recoverMapInv.end() && itp2 != recoverMapInv.end())
692  {
693  BDS_Point *pp1 = itp1->second;
694  BDS_Point *pp2 = itp2->second;
695  if ((pp1->iD == p1 && pp2->iD == p2) || (pp1->iD == p2 && pp2->iD == p1))
696  {
697  double t1;
698  double lc1 = -1;
699  if (v1->onWhat() == g1)
700  t1 = bb.low();
701  else if (v1->onWhat() == g2)
702  t1 = bb.high();
703  else
704  {
705  MEdgeVertex *ev1 = (MEdgeVertex *)v1;
706  lc1 = ev1->getLc();
707  v1->getParameter(0, t1);
708  }
709  double t2;
710  double lc2 = -1;
711  if (v2->onWhat() == g1)
712  t2 = bb.low();
713  else if (v2->onWhat() == g2)
714  t2 = bb.high();
715  else
716  {
717  MEdgeVertex *ev2 = (MEdgeVertex *)v2;
718  lc2 = ev2->getLc();
719  v2->getParameter(0, t2);
720  }
721 
722  // periodic
723  if (v1->onWhat() == g1 && v1->onWhat() == g2)
724  t1 = fabs(t2 - bb.low()) < fabs(t2 - bb.high()) ? bb.low() : bb.high();
725  if (v2->onWhat() == g1 && v2->onWhat() == g2)
726  t2 = fabs(t1 - bb.low()) < fabs(t1 - bb.high()) ? bb.low() : bb.high();
727 
728  if (lc1 == -1)
729  lc1 = BGM_MeshSize(v1->onWhat(), 0, 0, v1->x(), v1->y(), v1->z());
730  if (lc2 == -1)
731  lc2 = BGM_MeshSize(v2->onWhat(), 0, 0, v2->x(), v2->y(), v2->z());
732  // should be better, i.e. equidistant
733  double t = 0.5 * (t2 + t1);
734  double lc = 0.5 * (lc1 + lc2);
735  GPoint V = itr->ge->point(t);
736  MEdgeVertex *newv = new MEdgeVertex(V.x(), V.y(), V.z(), itr->ge, t, 0, lc);
737  newLines.push_back(new MLine(v1, newv));
738  newLines.push_back(new MLine(newv, v2));
739  delete itr->ge->lines[i];
740  }
741  else
742  {
743  newLines.push_back(itr->ge->lines[i]);
744  }
745  }
746  else
747  {
748  newLines.push_back(itr->ge->lines[i]);
749  }
750  }
751 
752  itr->ge->lines = newLines;
753  itr->ge->mesh_vertices.clear();
754  N = itr->ge->lines.size();
755  for (int i = 1; i < N; i++)
756  {
757  itr->ge->mesh_vertices.push_back(itr->ge->lines[i]->getVertex(0));
758  }
759  }
760 }
761 
762 static bool algoDelaunay2D(GFace *gf)
763 {
764  if (gf->getMeshingAlgo() == ALGO_2D_DELAUNAY || gf->getMeshingAlgo() == ALGO_2D_BAMG ||
765  gf->getMeshingAlgo() == ALGO_2D_FRONTAL || gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
766  gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS || gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR ||
767  gf->getMeshingAlgo() == ALGO_2D_BAMG)
768  return true;
769 
770  if (gf->getMeshingAlgo() == ALGO_2D_AUTO && gf->geomType() == GEntity::Plane)
771  return true;
772 
773  if (gf->getMeshingAlgo() == ALGO_2D_INITIAL_ONLY)
774  return true;
775 
776  return false;
777 }
778 
779 static bool recoverEdge(BDS_Mesh *m, GFace *gf, GEdge *ge, std::map<MVertex *, BDS_Point *> &recoverMapInv,
780  std::set<EdgeToRecover> *e2r, std::set<EdgeToRecover> *notRecovered, int pass)
781 {
782  BDS_GeomEntity *g = nullptr;
783  if (pass == 2)
784  {
785  m->add_geom(ge->tag(), 1);
786  g = m->get_geom(ge->tag(), 1);
787  }
788 
789  bool _fatallyFailed;
790 
791  for (std::size_t i = 0; i < ge->lines.size(); i++)
792  {
793  MVertex *vstart = ge->lines[i]->getVertex(0);
794  MVertex *vend = ge->lines[i]->getVertex(1);
795  auto itpstart = recoverMapInv.find(vstart);
796  auto itpend = recoverMapInv.find(vend);
797  if (itpstart != recoverMapInv.end() && itpend != recoverMapInv.end())
798  {
799  BDS_Point *pstart = itpstart->second;
800  BDS_Point *pend = itpend->second;
801  if (pass == 1)
802  e2r->insert(EdgeToRecover(pstart->iD, pend->iD, ge));
803  else
804  {
805  BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, _fatallyFailed, e2r, notRecovered);
806  if (e)
807  e->g = g;
808  else
809  {
810  if (_fatallyFailed)
811  {
812  Msg::Error("Unable to recover the edge %d (%d/%d) on curve %d (on "
813  "surface %d)",
814  ge->lines[i]->getNum(), i + 1, ge->lines.size(), ge->tag(), gf->tag());
815  if (Msg::GetVerbosity() == 99)
816  {
817  outputScalarField(m->triangles, "wrongmesh.pos", 0);
818  outputScalarField(m->triangles, "wrongparam.pos", 1);
819  }
820  }
821  return !_fatallyFailed;
822  }
823  }
824  }
825  }
826 
827  if (pass == 2 && ge->getBeginVertex())
828  {
829  MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
830  MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
831  auto itpstart = recoverMapInv.find(vstart);
832  auto itpend = recoverMapInv.find(vend);
833  if (itpstart != recoverMapInv.end() && itpend != recoverMapInv.end())
834  {
835  BDS_Point *pstart = itpstart->second;
836  BDS_Point *pend = itpend->second;
837  if (!pstart->g)
838  {
839  m->add_geom(pstart->iD, 0);
840  BDS_GeomEntity *g0 = m->get_geom(pstart->iD, 0);
841  pstart->g = g0;
842  }
843  if (!pend->g)
844  {
845  m->add_geom(pend->iD, 0);
846  BDS_GeomEntity *g0 = m->get_geom(pend->iD, 0);
847  pend->g = g0;
848  }
849  }
850  }
851 
852  return true;
853 }
854 
855 static void addOrRemove(MVertex *v1, MVertex *v2, std::set<MEdge, MEdgeLessThan> &bedges,
856  std::set<MEdge, MEdgeLessThan> &removed)
857 {
858  MEdge e(v1, v2);
859  if (removed.find(e) != removed.end())
860  return;
861  auto it = bedges.find(e);
862  if (it == bedges.end())
863  bedges.insert(e);
864  else
865  {
866  bedges.erase(it);
867  removed.insert(e);
868  }
869 }
870 
871 static bool meshGenerator(GFace *, int, bool, int, bool, std::vector<GEdge *> *);
872 
873 static void modifyInitialMeshForBoundaryLayers(GFace *gf, std::vector<MQuadrangle *> &blQuads,
874  std::vector<MTriangle *> &blTris, std::set<MVertex *> &verts, bool debug)
875 {
876  if (!buildAdditionalPoints2D(gf))
877  return;
878  BoundaryLayerColumns *_columns = gf->getColumns();
879 
880  std::set<MEdge, MEdgeLessThan> bedges;
881  std::set<MEdge, MEdgeLessThan> removed;
882 
883  std::vector<GEdge *> edges = gf->edges();
884  std::vector<GEdge *> embedded_edges = gf->getEmbeddedEdges();
885  edges.insert(edges.begin(), embedded_edges.begin(), embedded_edges.end());
886  auto ite = edges.begin();
887 
888  FILE *ff2 = nullptr;
889  if (debug)
890  ff2 = Fopen("tato.pos", "w");
891  if (ff2)
892  fprintf(ff2, "View \" \"{\n");
893 
894  std::vector<MLine *> _lines;
895 
896  while (ite != edges.end())
897  {
898  for (std::size_t i = 0; i < (*ite)->lines.size(); i++)
899  {
900  _lines.push_back((*ite)->lines[i]);
901  MVertex *v1 = (*ite)->lines[i]->getVertex(0);
902  MVertex *v2 = (*ite)->lines[i]->getVertex(1);
903  MEdge dv(v1, v2);
904  addOrRemove(v1, v2, bedges, removed);
905  for (std::size_t SIDE = 0; SIDE < _columns->_normals.count(dv); SIDE++)
906  {
907  std::vector<MElement *> myCol;
908  edgeColumn ec = _columns->getColumns(v1, v2, SIDE);
909  const BoundaryLayerData &c1 = ec._c1;
910  const BoundaryLayerData &c2 = ec._c2;
911  int N = std::min(c1._column.size(), c2._column.size());
912  if (!N)
913  continue;
914  for (int l = 0; l < N; ++l)
915  {
916  MVertex *v11, *v12, *v21, *v22;
917  v21 = c1._column[l];
918  v22 = c2._column[l];
919  if (l == 0)
920  {
921  v11 = v1;
922  v12 = v2;
923  }
924  else
925  {
926  v11 = c1._column[l - 1];
927  v12 = c2._column[l - 1];
928  }
929  MQuadrangle *qq = new MQuadrangle(v11, v21, v22, v12);
930  qq->setPartition(l + 1);
931  myCol.push_back(qq);
932  blQuads.push_back(qq);
933  if (ff2)
934  fprintf(ff2, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n", v11->x(), v11->y(),
935  v11->z(), v12->x(), v12->y(), v12->z(), v22->x(), v22->y(), v22->z(), v21->x(),
936  v21->y(), v21->z(), l + 1, l + 1, l + 1, l + 1);
937  }
938  if (c1._column.size() != c2._column.size())
939  {
940  MVertex *v11, *v12, *v;
941  v11 = c1._column[N - 1];
942  v12 = c2._column[N - 1];
943  v = c1._column.size() > c2._column.size() ? c1._column[N] : c2._column[N];
944  MTriangle *qq = new MTriangle(v11, v12, v);
945  qq->setPartition(N + 1);
946  myCol.push_back(qq);
947  blTris.push_back(qq);
948  if (ff2)
949  fprintf(ff2, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n", v11->x(), v11->y(), v11->z(),
950  v12->x(), v12->y(), v12->z(), v->x(), v->y(), v->z(), N + 1, N + 1, N + 1);
951  }
952  // int M = std::max(c1._column.size(),c2._column.size());
953  for (std::size_t l = 0; l < myCol.size(); l++)
954  _columns->_toFirst[myCol[l]] = myCol[0];
955  _columns->_elemColumns[myCol[0]] = myCol;
956  }
957  }
958  ++ite;
959  }
960 
961  for (auto itf = _columns->beginf(); itf != _columns->endf(); ++itf)
962  {
963  MVertex *v = itf->first;
964  int nbCol = _columns->getNbColumns(v);
965 
966  for (int i = 0; i < nbCol - 1; i++)
967  {
968  const BoundaryLayerData &c1 = _columns->getColumn(v, i);
969  const BoundaryLayerData &c2 = _columns->getColumn(v, i + 1);
970  int N = std::min(c1._column.size(), c2._column.size());
971  std::vector<MElement *> myCol;
972  for (int l = 0; l < N; ++l)
973  {
974  MVertex *v11, *v12, *v21, *v22;
975  v21 = c1._column[l];
976  v22 = c2._column[l];
977  if (l == 0)
978  {
979  v11 = v;
980  v12 = v;
981  }
982  else
983  {
984  v11 = c1._column[l - 1];
985  v12 = c2._column[l - 1];
986  }
987  if (v11 != v12)
988  {
989  MQuadrangle *qq = new MQuadrangle(v11, v12, v22, v21);
990  qq->setPartition(l + 1);
991  myCol.push_back(qq);
992  blQuads.push_back(qq);
993  if (ff2)
994  fprintf(ff2, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n", v11->x(), v11->y(),
995  v11->z(), v12->x(), v12->y(), v12->z(), v22->x(), v22->y(), v22->z(), v21->x(),
996  v21->y(), v21->z(), l + 1, l + 1, l + 1, l + 1);
997  }
998  else
999  {
1000  MTriangle *qq = new MTriangle(v, v22, v21);
1001  qq->setPartition(l + 1);
1002  myCol.push_back(qq);
1003  blTris.push_back(qq);
1004  if (ff2)
1005  fprintf(ff2, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n", v->x(), v->y(), v->z(), v22->x(),
1006  v22->y(), v22->z(), v21->x(), v21->y(), v21->z(), l + 1, l + 1, l + 1);
1007  }
1008  }
1009  if (myCol.size())
1010  {
1011  for (std::size_t l = 0; l < myCol.size(); l++)
1012  _columns->_toFirst[myCol[l]] = myCol[0];
1013  _columns->_elemColumns[myCol[0]] = myCol;
1014  }
1015  }
1016  }
1017 
1018  if (ff2)
1019  {
1020  fprintf(ff2, "};\n");
1021  fclose(ff2);
1022  }
1023 
1024  filterOverlappingElements(_lines, blTris, blQuads, _columns->_elemColumns, _columns->_toFirst);
1025  for (std::size_t i = 0; i < blQuads.size(); i++)
1026  blQuads[i]->setPartition(0);
1027  for (std::size_t i = 0; i < blTris.size(); i++)
1028  blTris[i]->setPartition(0);
1029 
1030  for (std::size_t i = 0; i < blQuads.size(); i++)
1031  {
1032  addOrRemove(blQuads[i]->getVertex(0), blQuads[i]->getVertex(1), bedges, removed);
1033  addOrRemove(blQuads[i]->getVertex(1), blQuads[i]->getVertex(2), bedges, removed);
1034  addOrRemove(blQuads[i]->getVertex(2), blQuads[i]->getVertex(3), bedges, removed);
1035  addOrRemove(blQuads[i]->getVertex(3), blQuads[i]->getVertex(0), bedges, removed);
1036  for (int j = 0; j < 4; j++)
1037  if (blQuads[i]->getVertex(j)->onWhat() == gf)
1038  verts.insert(blQuads[i]->getVertex(j));
1039  }
1040  for (std::size_t i = 0; i < blTris.size(); i++)
1041  {
1042  addOrRemove(blTris[i]->getVertex(0), blTris[i]->getVertex(1), bedges, removed);
1043  addOrRemove(blTris[i]->getVertex(1), blTris[i]->getVertex(2), bedges, removed);
1044  addOrRemove(blTris[i]->getVertex(2), blTris[i]->getVertex(0), bedges, removed);
1045  for (int j = 0; j < 3; j++)
1046  if (blTris[i]->getVertex(j)->onWhat() == gf)
1047  verts.insert(blTris[i]->getVertex(j));
1048  }
1049 
1050  discreteEdge ne(gf->model(), 444444, nullptr, (*edges.begin())->getEndVertex());
1051  std::vector<GEdge *> hop;
1052  auto it = bedges.begin();
1053 
1054  FILE *ff = nullptr;
1055  if (debug)
1056  ff = Fopen("toto.pos", "w");
1057  if (ff)
1058  fprintf(ff, "View \" \"{\n");
1059  for (; it != bedges.end(); ++it)
1060  {
1061  ne.lines.push_back(new MLine(it->getVertex(0), it->getVertex(1)));
1062  if (ff)
1063  fprintf(ff, "SL(%g,%g,%g,%g,%g,%g){1,1};\n", it->getVertex(0)->x(), it->getVertex(0)->y(),
1064  it->getVertex(0)->z(), it->getVertex(1)->x(), it->getVertex(1)->y(), it->getVertex(1)->z());
1065  }
1066  if (ff)
1067  {
1068  fprintf(ff, "};\n");
1069  fclose(ff);
1070  }
1071  hop.push_back(&ne);
1072 
1073  deMeshGFace kil_;
1074  kil_(gf);
1075  meshGenerator(gf, 0, false, 99, false, &hop);
1076 }
1077 
1078 static bool improved_translate(GFace *gf, MVertex *vertex, SVector3 &v1, SVector3 &v2)
1079 {
1080  double x, y;
1081  double angle;
1082  SPoint2 point;
1083  SVector3 s1, s2;
1084  SVector3 normal;
1085  SVector3 basis_u, basis_v;
1086  Pair<SVector3, SVector3> derivatives;
1087 
1088  reparamMeshVertexOnFace(vertex, gf, point);
1089  x = point.x();
1090  y = point.y();
1091 
1092  angle = backgroundMesh::current()->getAngle(x, y, 0.0);
1093  derivatives = gf->firstDer(point);
1094 
1095  s1 = derivatives.first();
1096  s2 = derivatives.second();
1097  normal = crossprod(s1, s2);
1098 
1099  basis_u = s1;
1100  basis_u.normalize();
1101  basis_v = crossprod(normal, basis_u);
1102  basis_v.normalize();
1103 
1104  v1 = basis_u * cos(angle) + basis_v * sin(angle);
1105  v2 = crossprod(v1, normal);
1106  v2.normalize();
1107 
1108  return 1;
1109 }
1110 
1111 static void directions_storage(GFace *gf)
1112 {
1113  std::set<MVertex *> vertices;
1114  for (std::size_t i = 0; i < gf->getNumMeshElements(); i++)
1115  {
1116  MElement *element = gf->getMeshElement(i);
1117  for (std::size_t j = 0; j < element->getNumVertices(); j++)
1118  {
1119  MVertex *vertex = element->getVertex(j);
1120  vertices.insert(vertex);
1121  }
1122  }
1123 
1124  backgroundMesh::set(gf);
1125 
1126  gf->storage1.clear();
1127  gf->storage2.clear();
1128  gf->storage3.clear();
1129  gf->storage4.clear();
1130 
1131  for (auto it = vertices.begin(); it != vertices.end(); it++)
1132  {
1133  SPoint2 point;
1134  SVector3 v1;
1135  SVector3 v2;
1136  bool ok = improved_translate(gf, *it, v1, v2);
1137  if (ok)
1138  {
1139  gf->storage1.push_back(SPoint3((*it)->x(), (*it)->y(), (*it)->z()));
1140  gf->storage2.push_back(v1);
1141  gf->storage3.push_back(v2);
1142  reparamMeshVertexOnFace(*it, gf, point);
1143  gf->storage4.push_back(backgroundMesh::current()->operator()(point.x(), point.y(), 0.0));
1144  }
1145  }
1146 
1147  backgroundMesh::unset();
1148 }
1149 
1150 static void BDS2GMSH(BDS_Mesh *m, GFace *gf, std::map<BDS_Point *, MVertex *, PointLessThan> &recoverMap)
1151 {
1152  auto itt = m->triangles.begin();
1153  while (itt != m->triangles.end())
1154  {
1155  BDS_Face *t = *itt;
1156  if (!t->deleted)
1157  {
1158  BDS_Point *n[4];
1159  t->getNodes(n);
1160  MVertex *v[4] = {nullptr, nullptr, nullptr, nullptr};
1161  for (int i = 0; i < 4; i++)
1162  {
1163  if (!n[i])
1164  continue;
1165  if (recoverMap.find(n[i]) == recoverMap.end())
1166  {
1167  v[i] = new MFaceVertex(n[i]->X, n[i]->Y, n[i]->Z, gf, n[i]->u, n[i]->v);
1168  recoverMap[n[i]] = v[i];
1169  gf->mesh_vertices.push_back(v[i]);
1170  }
1171  else
1172  v[i] = recoverMap[n[i]];
1173  }
1174  if (!v[3])
1175  {
1176  // when a singular point is present, degenerated triangles may be
1177  // created, for example on a sphere that contains one pole
1178  if (v[0] != v[1] && v[0] != v[2] && v[1] != v[2])
1179  gf->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
1180  }
1181  else
1182  {
1183  gf->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3]));
1184  }
1185  }
1186  ++itt;
1187  }
1188 }
1189 
1190 static void deleteUnusedVertices(GFace *gf)
1191 {
1192  std::set<MVertex *, MVertexPtrLessThan> allverts;
1193  for (std::size_t i = 0; i < gf->triangles.size(); i++)
1194  {
1195  for (int j = 0; j < 3; j++)
1196  {
1197  if (gf->triangles[i]->getVertex(j)->onWhat() == gf)
1198  allverts.insert(gf->triangles[i]->getVertex(j));
1199  }
1200  }
1201  for (std::size_t i = 0; i < gf->quadrangles.size(); i++)
1202  {
1203  for (int j = 0; j < 4; j++)
1204  {
1205  if (gf->quadrangles[i]->getVertex(j)->onWhat() == gf)
1206  allverts.insert(gf->quadrangles[i]->getVertex(j));
1207  }
1208  }
1209  for (std::size_t i = 0; i < gf->mesh_vertices.size(); i++)
1210  {
1211  if (allverts.find(gf->mesh_vertices[i]) == allverts.end())
1212  delete gf->mesh_vertices[i];
1213  }
1214  gf->mesh_vertices.clear();
1215  gf->mesh_vertices.insert(gf->mesh_vertices.end(), allverts.begin(), allverts.end());
1216 }
1217 
1218 // Builds An initial triangular mesh that respects the boundaries of
1219 // the domain, including embedded points and surfaces
1220 static bool meshGenerator(GFace *gf, int RECUR_ITER, bool repairSelfIntersecting1dMesh, int onlyInitialMesh, bool debug,
1221  std::vector<GEdge *> *replacementEdges)
1222 {
1223  if (CTX::instance()->debugSurface > 0 && gf->tag() != CTX::instance()->debugSurface)
1224  {
1225  gf->meshStatistics.status = GFace::DONE;
1226  return true;
1227  }
1228  if (CTX::instance()->debugSurface > 0)
1229  debug = true;
1230 
1231  BDS_GeomEntity CLASS_F(1, 2);
1232  BDS_GeomEntity CLASS_EXTERIOR(1, 3);
1233  std::map<BDS_Point *, MVertex *, PointLessThan> recoverMap;
1234  std::map<MVertex *, BDS_Point *> recoverMapInv;
1235  std::vector<GEdge *> edges = replacementEdges ? *replacementEdges : gf->edges();
1236 
1237  FILE *fdeb = nullptr;
1238  if (debug && RECUR_ITER == 0)
1239  {
1240  char name[245];
1241  sprintf(name, "surface%d-boundary-real.pos", gf->tag());
1242  fdeb = fopen(name, "w");
1243  fprintf(fdeb, "View \"\"{\n");
1244  }
1245 
1246  // build a set with all points of the boundaries
1247  std::set<MVertex *, MVertexPtrLessThan> all_vertices, boundary;
1248  auto ite = edges.begin();
1249  while (ite != edges.end())
1250  {
1251  if ((*ite)->isSeam(gf))
1252  {
1253  if (fdeb != nullptr)
1254  fclose(fdeb);
1255  return false;
1256  }
1257  if (!(*ite)->isMeshDegenerated())
1258  {
1259  for (std::size_t i = 0; i < (*ite)->lines.size(); i++)
1260  {
1261  MVertex *v1 = (*ite)->lines[i]->getVertex(0);
1262  MVertex *v2 = (*ite)->lines[i]->getVertex(1);
1263 
1264  if (fdeb)
1265  {
1266  fprintf(fdeb, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(),
1267  v2->z(), (*ite)->tag(), (*ite)->tag());
1268  }
1269 
1270  all_vertices.insert(v1);
1271  all_vertices.insert(v2);
1272  if (boundary.find(v1) == boundary.end())
1273  boundary.insert(v1);
1274  else
1275  boundary.erase(v1);
1276  if (boundary.find(v2) == boundary.end())
1277  boundary.insert(v2);
1278  else
1279  boundary.erase(v2);
1280  }
1281  }
1282  else
1283  Msg::Debug("Degenerated mesh on edge %d", (*ite)->tag());
1284  ++ite;
1285  }
1286 
1287  if (fdeb)
1288  {
1289  fprintf(fdeb, "};\n");
1290  fclose(fdeb);
1291  }
1292 
1293  if (boundary.size())
1294  {
1295  Msg::Error("The 1D mesh seems not to be forming a closed loop (%d boundary "
1296  "nodes are considered once)",
1297  boundary.size());
1298  for (auto it = boundary.begin(); it != boundary.end(); it++)
1299  {
1300  Msg::Debug("Remaining node %lu", (*it)->getNum());
1301  }
1302  gf->meshStatistics.status = GFace::FAILED;
1303  return false;
1304  }
1305 
1306  std::vector<GEdge *> emb_edges = gf->getEmbeddedEdges();
1307  ite = emb_edges.begin();
1308  while (ite != emb_edges.end())
1309  {
1310  if (!(*ite)->isMeshDegenerated())
1311  {
1312  all_vertices.insert((*ite)->mesh_vertices.begin(), (*ite)->mesh_vertices.end());
1313  if ((*ite)->getBeginVertex())
1314  all_vertices.insert((*ite)->getBeginVertex()->mesh_vertices.begin(),
1315  (*ite)->getBeginVertex()->mesh_vertices.end());
1316  if ((*ite)->getEndVertex())
1317  all_vertices.insert((*ite)->getEndVertex()->mesh_vertices.begin(),
1318  (*ite)->getEndVertex()->mesh_vertices.end());
1319  }
1320  ++ite;
1321  }
1322 
1323  // add embedded vertices
1324  std::vector<GVertex *> emb_vertx = gf->getEmbeddedVertices();
1325  auto itvx = emb_vertx.begin();
1326  while (itvx != emb_vertx.end())
1327  {
1328  all_vertices.insert((*itvx)->mesh_vertices.begin(), (*itvx)->mesh_vertices.end());
1329  ++itvx;
1330  }
1331 
1332  if (all_vertices.size() < 3)
1333  {
1334  Msg::Warning("Mesh generation of surface %d skipped: only %d nodes on "
1335  "the boundary",
1336  gf->tag(), all_vertices.size());
1337  gf->meshStatistics.status = GFace::DONE;
1338  return true;
1339  }
1340  else if (all_vertices.size() == 3)
1341  {
1342  MVertex *vv[3] = {nullptr, nullptr, nullptr};
1343  int i = 0;
1344  for (auto it = all_vertices.begin(); it != all_vertices.end(); it++)
1345  {
1346  vv[i++] = *it;
1347  }
1348  gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
1349  gf->meshStatistics.status = GFace::DONE;
1350  return true;
1351  }
1352 
1353  // Buid a BDS_Mesh structure that is convenient for doing the actual
1354  // meshing procedure
1355  BDS_Mesh *m = new BDS_Mesh;
1356 
1357  std::vector<BDS_Point *> points(all_vertices.size());
1358  SBoundingBox3d bbox;
1359  int count = 0;
1360  for (auto it = all_vertices.begin(); it != all_vertices.end(); it++)
1361  {
1362  MVertex *here = *it;
1363  GEntity *ge = here->onWhat();
1364  SPoint2 param;
1365  reparamMeshVertexOnFace(here, gf, param);
1366  BDS_Point *pp = m->add_point(count, param[0], param[1], gf);
1367  m->add_geom(ge->tag(), ge->dim());
1368  BDS_GeomEntity *g = m->get_geom(ge->tag(), ge->dim());
1369  pp->g = g;
1370  recoverMap[pp] = here;
1371  recoverMapInv[here] = pp;
1372  points[count] = pp;
1373  bbox += SPoint3(param[0], param[1], 0);
1374  count++;
1375  }
1376 
1377  bbox.makeCube();
1378 
1379  // use a divide & conquer type algorithm to create a triangulation.
1380  // We add to the triangulation a box with 4 points that encloses the
1381  // domain.
1382  if (CTX::instance()->mesh.oldInitialDelaunay2D)
1383  {
1384  // compute the bounding box in parametric space
1385  SVector3 dd(bbox.max(), bbox.min());
1386  double LC2D = norm(dd);
1387  DocRecord doc(points.size() + 4);
1388  for (std::size_t i = 0; i < points.size(); i++)
1389  {
1390  double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
1391  double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
1392  doc.points[i].where.h = points[i]->u + XX;
1393  doc.points[i].where.v = points[i]->v + YY;
1394  doc.points[i].data = points[i];
1395  doc.points[i].adjacent = nullptr;
1396  }
1397 
1398  // increase the size of the bounding box
1399  bbox *= 2.5;
1400 
1401  // add 4 points than encloses the domain (use negative number to
1402  // distinguish those fake vertices)
1403  double bb[4][2] = {{bbox.min().x(), bbox.min().y()},
1404  {bbox.min().x(), bbox.max().y()},
1405  {bbox.max().x(), bbox.min().y()},
1406  {bbox.max().x(), bbox.max().y()}};
1407  for (int ip = 0; ip < 4; ip++)
1408  {
1409  BDS_Point *pp = m->add_point(-ip - 1, bb[ip][0], bb[ip][1], gf);
1410  m->add_geom(gf->tag(), 2);
1411  BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
1412  pp->g = g;
1413  doc.points[points.size() + ip].where.h = bb[ip][0];
1414  doc.points[points.size() + ip].where.v = bb[ip][1];
1415  doc.points[points.size() + ip].adjacent = nullptr;
1416  doc.points[points.size() + ip].data = pp;
1417  }
1418 
1419  // Use "fast" inhouse recursive algo to generate the triangulation.
1420  // At this stage the triangulation is not what we need
1421  // -) It does not necessary recover the boundaries
1422  // -) It contains triangles outside the domain (the first edge
1423  // loop is the outer one)
1424  Msg::Debug("Meshing of the convex hull (%d points)", points.size());
1425  try
1426  {
1427  doc.MakeMeshWithPoints();
1428  }
1429  catch (std::runtime_error &e)
1430  {
1431  Msg::Error("%s", e.what());
1432  }
1433  Msg::Debug("Meshing of the convex hull (%d points) done", points.size());
1434 
1435  for (int i = 0; i < doc.numTriangles; i++)
1436  {
1437  int a = doc.triangles[i].a;
1438  int b = doc.triangles[i].b;
1439  int c = doc.triangles[i].c;
1440  int n = doc.numPoints;
1441  if (a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n)
1442  {
1443  Msg::Warning("Skipping bad triangle %d", i);
1444  continue;
1445  }
1446  BDS_Point *p1 = (BDS_Point *)doc.points[doc.triangles[i].a].data;
1447  BDS_Point *p2 = (BDS_Point *)doc.points[doc.triangles[i].b].data;
1448  BDS_Point *p3 = (BDS_Point *)doc.points[doc.triangles[i].c].data;
1449  m->add_triangle(p1->iD, p2->iD, p3->iD);
1450  }
1451  }
1452  else
1453  {
1454  // New initial 2D mesher --> it actually does everything
1455  // -- triangulate points
1456  // -- recover edges
1457  // -- color triangles
1458 
1459  std::vector<GEdge *> temp;
1460  if (replacementEdges)
1461  {
1462  temp = gf->edges();
1463  gf->set(*replacementEdges);
1464  }
1465  // recover and color so most of the code below can go away. Works also for
1466  // periodic faces
1467  PolyMesh *pm = GFaceInitialMesh(gf->tag(), 1);
1468 
1469  if (replacementEdges)
1470  {
1471  gf->set(temp);
1472  }
1473 
1474  std::map<int, BDS_Point *> aaa;
1475  for (auto it = all_vertices.begin(); it != all_vertices.end(); it++)
1476  {
1477  MVertex *here = *it;
1478  aaa[here->getNum()] = recoverMapInv[here];
1479  }
1480 
1481  for (int ip = 0; ip < 4; ip++)
1482  {
1483  PolyMesh::Vertex *v = pm->vertices[ip];
1484  v->data = -ip - 1;
1485  BDS_Point *pp = m->add_point(v->data, v->position.x(), v->position.y(), gf);
1486  m->add_geom(gf->tag(), 2);
1487  BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
1488  pp->g = g;
1489  aaa[v->data] = pp;
1490  }
1491 
1492  for (size_t i = 0; i < pm->faces.size(); i++)
1493  {
1494  PolyMesh::HalfEdge *he = pm->faces[i]->he;
1495  int a = he->v->data;
1496  int b = he->next->v->data;
1497  int c = he->next->next->v->data;
1498  BDS_Point *p1 = aaa[a];
1499  BDS_Point *p2 = aaa[b];
1500  BDS_Point *p3 = aaa[c];
1501  if (p1 && p2 && p3)
1502  m->add_triangle(p1->iD, p2->iD, p3->iD);
1503  }
1504  delete pm;
1505  }
1506 
1507  if (debug && RECUR_ITER == 0)
1508  {
1509  char name[245];
1510  sprintf(name, "surface%d-initial-real.pos", gf->tag());
1511  outputScalarField(m->triangles, name, 0, gf);
1512  sprintf(name, "surface%d-initial-param.pos", gf->tag());
1513  outputScalarField(m->triangles, name, 1, gf);
1514  }
1515 
1516  // Recover the boundary edges and compute characteristic lenghts using mesh
1517  // edge spacing. If two of these edges intersect, then the 1D mesh have to be
1518  // densified
1519  Msg::Debug("Recovering %d model edges", edges.size());
1520  std::set<EdgeToRecover> edgesToRecover;
1521  std::set<EdgeToRecover> edgesNotRecovered;
1522  ite = edges.begin();
1523  while (ite != edges.end())
1524  {
1525  if (!(*ite)->isMeshDegenerated())
1526  recoverEdge(m, gf, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
1527  ++ite;
1528  }
1529  ite = emb_edges.begin();
1530  while (ite != emb_edges.end())
1531  {
1532  if (!(*ite)->isMeshDegenerated())
1533  recoverEdge(m, gf, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
1534  ++ite;
1535  }
1536 
1537  // effectively recover the medge
1538  ite = edges.begin();
1539  while (ite != edges.end())
1540  {
1541  if (!(*ite)->isMeshDegenerated())
1542  {
1543  if (!recoverEdge(m, gf, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2))
1544  {
1545  delete m;
1546  gf->meshStatistics.status = GFace::FAILED;
1547  return false;
1548  }
1549  }
1550  ++ite;
1551  }
1552 
1553  Msg::Debug("Recovering %d mesh edges (%d not recovered)", edgesToRecover.size(), edgesNotRecovered.size());
1554 
1555  if (edgesNotRecovered.size() || gf->meshStatistics.refineAllEdges)
1556  {
1557  std::ostringstream sstream;
1558  for (auto itr = edgesNotRecovered.begin(); itr != edgesNotRecovered.end(); ++itr)
1559  sstream << " " << itr->ge->tag();
1560  if (gf->meshStatistics.refineAllEdges)
1561  {
1562  Msg::Info("8-| Splitting all edges and trying again");
1563  }
1564  else
1565  {
1566  Msg::Info(":-( There are %d intersections in the 1D mesh (curves%s)", edgesNotRecovered.size(),
1567  sstream.str().c_str());
1568  if (repairSelfIntersecting1dMesh)
1569  Msg::Info("8-| Splitting those edges and trying again");
1570  }
1571  if (debug)
1572  {
1573  char name[245];
1574  sprintf(name, "surface%d-not_yet_recovered-real-%d.msh", gf->tag(), RECUR_ITER);
1575  gf->model()->writeMSH(name);
1576  }
1577 
1578  if (repairSelfIntersecting1dMesh)
1579  {
1580  remeshUnrecoveredEdges(recoverMapInv, edgesNotRecovered);
1581  gf->meshStatistics.refineAllEdges = false;
1582  }
1583  else
1584  {
1585  auto itr = edgesNotRecovered.begin();
1586  int I = 0;
1587  for (; itr != edgesNotRecovered.end(); ++itr)
1588  {
1589  int p1 = itr->p1;
1590  int p2 = itr->p2;
1591  int tag = itr->ge->tag();
1592  Msg::Error("Edge not recovered: %d %d %d", p1, p2, tag);
1593  I++;
1594  }
1595  }
1596 
1597  // delete the mesh
1598  delete m;
1599  if (RECUR_ITER < 10)
1600  {
1601  return meshGenerator(gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh, debug,
1602  replacementEdges);
1603  }
1604  return false;
1605  }
1606 
1607  if (RECUR_ITER > 0)
1608  Msg::Info(":-) All edges recovered after %d iteration%s", RECUR_ITER, (RECUR_ITER > 1) ? "s" : "");
1609 
1610  Msg::Debug("Boundary edges recovered for surface %d", gf->tag());
1611 
1612  // look for a triangle that has a negative node and recursively tag all
1613  // exterior triangles
1614  {
1615  auto itt = m->triangles.begin();
1616  while (itt != m->triangles.end())
1617  {
1618  (*itt)->g = nullptr;
1619  ++itt;
1620  }
1621  itt = m->triangles.begin();
1622  while (itt != m->triangles.end())
1623  {
1624  BDS_Face *t = *itt;
1625  BDS_Point *n[4];
1626  if (t->getNodes(n))
1627  {
1628  if (n[0]->iD < 0 || n[1]->iD < 0 || n[2]->iD < 0)
1629  {
1630  recur_tag(t, &CLASS_EXTERIOR);
1631  break;
1632  }
1633  }
1634  ++itt;
1635  }
1636  }
1637 
1638  // now find an edge that has belongs to one of the exterior triangles
1639  {
1640  auto ite = m->edges.begin();
1641  while (ite != m->edges.end())
1642  {
1643  BDS_Edge *e = *ite;
1644  if (e->g && e->numfaces() == 2)
1645  {
1646  if (e->faces(0)->g == &CLASS_EXTERIOR)
1647  {
1648  recur_tag(e->faces(1), &CLASS_F);
1649  break;
1650  }
1651  else if (e->faces(1)->g == &CLASS_EXTERIOR)
1652  {
1653  recur_tag(e->faces(0), &CLASS_F);
1654  break;
1655  }
1656  }
1657  ++ite;
1658  }
1659  auto itt = m->triangles.begin();
1660  while (itt != m->triangles.end())
1661  {
1662  if ((*itt)->g == &CLASS_EXTERIOR)
1663  (*itt)->g = nullptr;
1664  ++itt;
1665  }
1666  }
1667 
1668  {
1669  auto ite = m->edges.begin();
1670  while (ite != m->edges.end())
1671  {
1672  BDS_Edge *e = *ite;
1673  if (e->g && e->numfaces() == 2)
1674  {
1675  BDS_Point *oface[2];
1676  e->oppositeof(oface);
1677  if (oface[0]->iD < 0)
1678  {
1679  recur_tag(e->faces(1), &CLASS_F);
1680  break;
1681  }
1682  else if (oface[1]->iD < 0)
1683  {
1684  recur_tag(e->faces(0), &CLASS_F);
1685  break;
1686  }
1687  }
1688  ++ite;
1689  }
1690  }
1691 
1692  ite = emb_edges.begin();
1693  while (ite != emb_edges.end())
1694  {
1695  if (!(*ite)->isMeshDegenerated())
1696  recoverEdge(m, gf, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
1697  ++ite;
1698  }
1699 
1700  // compute characteristic lengths at vertices
1701  if (CTX::instance()->mesh.algo2d != ALGO_2D_BAMG && !onlyInitialMesh)
1702  {
1703  Msg::Debug("Computing mesh size field at mesh nodes %d", edgesToRecover.size());
1704  auto it = m->points.begin();
1705  for (; it != m->points.end(); ++it)
1706  {
1707  BDS_Point *pp = *it;
1708  auto itv = recoverMap.find(pp);
1709  if (itv != recoverMap.end())
1710  {
1711  MVertex *here = itv->second;
1712  GEntity *ge = here->onWhat();
1713  if (ge->dim() == 0)
1714  {
1715  pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
1716  }
1717  else if (ge->dim() == 1)
1718  {
1719  double u;
1720  here->getParameter(0, u);
1721  pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
1722  }
1723  else
1724  pp->lcBGM() = MAX_LC;
1725  pp->lc() = pp->lcBGM();
1726  }
1727  }
1728  }
1729 
1730  // delete useless stuff
1731  auto itt = m->triangles.begin();
1732  while (itt != m->triangles.end())
1733  {
1734  BDS_Face *t = *itt;
1735  if (!t->g)
1736  m->del_face(t);
1737  ++itt;
1738  }
1739  m->cleanup();
1740 
1741  {
1742  auto ite = m->edges.begin();
1743  while (ite != m->edges.end())
1744  {
1745  BDS_Edge *e = *ite;
1746  if (e->numfaces() == 0)
1747  m->del_edge(e);
1748  else
1749  {
1750  if (!e->g)
1751  e->g = &CLASS_F;
1752  if (!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree)
1753  e->p1->g = e->g;
1754  if (!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree)
1755  e->p2->g = e->g;
1756  }
1757  ++ite;
1758  }
1759  }
1760  m->cleanup();
1761  m->del_point(m->find_point(-1));
1762  m->del_point(m->find_point(-2));
1763  m->del_point(m->find_point(-3));
1764  m->del_point(m->find_point(-4));
1765 
1766  if (debug)
1767  {
1768  char name[245];
1769  sprintf(name, "surface%d-recovered-real.pos", gf->tag());
1770  outputScalarField(m->triangles, name, 0, gf);
1771  sprintf(name, "surface%d-recovered-param.pos", gf->tag());
1772  outputScalarField(m->triangles, name, 1, gf);
1773  }
1774 
1775  if (1)
1776  {
1777  auto itt = m->triangles.begin();
1778  while (itt != m->triangles.end())
1779  {
1780  BDS_Face *t = *itt;
1781  if (!t->deleted)
1782  {
1783  BDS_Point *n[4];
1784  if (t->getNodes(n))
1785  {
1786  MVertex *v1 = recoverMap[n[0]];
1787  MVertex *v2 = recoverMap[n[1]];
1788  MVertex *v3 = recoverMap[n[2]];
1789  if (!n[3])
1790  {
1791  if (v1 != v2 && v1 != v3 && v2 != v3)
1792  gf->triangles.push_back(new MTriangle(v1, v2, v3));
1793  }
1794  else
1795  {
1796  MVertex *v4 = recoverMap[n[3]];
1797  gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
1798  }
1799  }
1800  }
1801  ++itt;
1802  }
1803  }
1804 
1805  {
1806  int nb_swap;
1807  Msg::Debug("Delaunizing the initial mesh");
1808  delaunayizeBDS(gf, *m, nb_swap);
1809  }
1810 
1811  // only delete the mesh data stored in the base GFace class
1812  gf->GFace::deleteMesh();
1813 
1814  Msg::Debug("Starting to add internal nodes");
1815  // start mesh generation
1816  if (!algoDelaunay2D(gf) && !onlyInitialMesh)
1817  {
1818  refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true, &recoverMapInv, nullptr);
1819  refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, false, &recoverMapInv, nullptr);
1820  }
1821 
1822  gf->meshStatistics.status = GFace::DONE;
1823 
1824  // fill the small gmsh structures
1825  BDS2GMSH(m, gf, recoverMap);
1826 
1827  std::vector<MQuadrangle *> blQuads;
1828  std::vector<MTriangle *> blTris;
1829  std::set<MVertex *> verts;
1830 
1831  bool infty = false;
1832  if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD)
1833  {
1834  infty = true;
1835  if (!onlyInitialMesh)
1836  buildBackgroundMesh(gf, CTX::instance()->mesh.crossFieldClosestPoint);
1837  }
1838  else if (gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS || gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
1839  {
1840  infty = true;
1841  /* New version of PACK / QUADQS use a different background mesh */
1842  // if(!onlyInitialMesh) buildBackgroundMesh(gf, false);
1843  }
1844 
1845  if (!onlyInitialMesh)
1846  modifyInitialMeshForBoundaryLayers(gf, blQuads, blTris, verts, debug);
1847 
1848  // the delaunay algo is based directly on internal gmsh structures BDS mesh is
1849  // passed in order not to recompute local coordinates of vertices
1850  if (algoDelaunay2D(gf) && !onlyInitialMesh)
1851  {
1852  if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL)
1853  {
1854  bowyerWatsonFrontal(gf);
1855  }
1856  else if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD)
1857  {
1858  bowyerWatsonFrontalLayers(gf, true);
1859  }
1860  else if (gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
1861  {
1862  bowyerWatsonParallelograms(gf);
1863  }
1864  else if (gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
1865  {
1866  Msg::Error("ALGO_2D_PACK_PRLGRMS_CSTR deprecated");
1867  // bowyerWatsonParallelogramsConstrained(gf, gf->constr_vertices);
1868  }
1869  else if (gf->getMeshingAlgo() == ALGO_2D_DELAUNAY || gf->getMeshingAlgo() == ALGO_2D_AUTO)
1870  {
1871  bowyerWatson(gf);
1872  }
1873  else
1874  {
1875  bowyerWatson(gf, 15000);
1876  meshGFaceBamg(gf);
1877  }
1878 
1879  if (!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
1880  laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, infty);
1881  }
1882 
1883  if (debug)
1884  {
1885  char name[256];
1886  // sprintf(name, "trueBoundary%d.pos", gf->tag());
1887  // std::vector<SPoint2> bnd;
1888  // trueBoundary(name, gf,bnd);
1889  sprintf(name, "real%d.pos", gf->tag());
1890  outputScalarField(m->triangles, name, 0, gf);
1891  sprintf(name, "param%d.pos", gf->tag());
1892  outputScalarField(m->triangles, name, 1, gf);
1893  }
1894 
1895  delete m;
1896 
1897  gf->quadrangles.insert(gf->quadrangles.begin(), blQuads.begin(), blQuads.end());
1898  gf->triangles.insert(gf->triangles.begin(), blTris.begin(), blTris.end());
1899  gf->mesh_vertices.insert(gf->mesh_vertices.begin(), verts.begin(), verts.end());
1900 
1901  splitElementsInBoundaryLayerIfNeeded(gf);
1902 
1903  if ((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) && (onlyInitialMesh != 99) &&
1904  (CTX::instance()->mesh.algoRecombine <= 1 || CTX::instance()->mesh.algoRecombine == 4))
1905  {
1906 
1907  if (CTX::instance()->mesh.algoRecombine == 4)
1908  {
1909  meshGFaceQuadrangulateBipartiteLabelling(gf->tag());
1910  }
1911  else
1912  {
1913  bool blossom = (CTX::instance()->mesh.algoRecombine == 1);
1914  int topo = CTX::instance()->mesh.recombineOptimizeTopology;
1915  int repos = CTX::instance()->mesh.recombineNodeRepositioning;
1916  double minqual = CTX::instance()->mesh.recombineMinimumQuality;
1917  recombineIntoQuads(gf, blossom, topo, repos, minqual);
1918  }
1919  }
1920 
1921  computeElementShapes(gf, gf->meshStatistics.worst_element_shape, gf->meshStatistics.average_element_shape,
1922  gf->meshStatistics.best_element_shape, gf->meshStatistics.nbTriangle,
1923  gf->meshStatistics.nbGoodQuality);
1924 
1925  if (CTX::instance()->mesh.algo3d == ALGO_3D_RTREE)
1926  {
1927  directions_storage(gf);
1928  }
1929 
1930  // remove unused vertices, generated e.g. during background mesh
1931  deleteUnusedVertices(gf);
1932 
1933  return true;
1934 }
1935 
1936 // this function buils a list of BDS points that are consecutive in one given
1937 // edge loop, taking care of periodic surfaces with seams; it also fills the
1938 // recoverMap to link BDS points with mesh nodes
1939 
1940 static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, std::vector<BDS_Point *> &result,
1941  SBoundingBox3d &bbox, BDS_Mesh *m,
1942  std::map<BDS_Point *, MVertex *, PointLessThan> &recoverMap, int &count,
1943  int countTot, double tol, bool seam_the_first = false)
1944 {
1945  result.clear();
1946  count = 0;
1947 
1948  // building the list of nodes in the parametric space for periodic surfaces is
1949  // tricky:
1950  // 1) both reparametrizations of the seams must be tested, as the topological
1951  // representation does not know which one to use
1952  // 2) both orientations of curves must be tested, as the topological
1953  // representation of periodic curves cannot distinguish them
1954  // 3) reparametrization of curves on surfaces can lead to slightly different
1955  // parametric coordinates, especially with OpenCASCADE - so a tolerance is
1956  // needed
1957  std::vector<GEdgeSigned> signedEdges(gel.begin(), gel.end());
1958  std::vector<SPoint2> coords;
1959  std::vector<MVertex *> verts;
1960 
1961 #if 0 // for debugging - don't remove
1962  printf("curve loop for surface %d\n", gf->tag());
1963  for(std::size_t i = 0; i < signedEdges.size(); i++) {
1964  GEdge *ge = signedEdges[i].getEdge();
1965  bool seam = ge->isSeam(gf);
1966  Range<double> range = ge->parBoundsOnFace(gf);
1967  printf(" curve %d: ", ge->tag());
1968  SPoint2 p;
1969  p = ge->reparamOnFace(gf, range.low(), 1);
1970  printf("beg (%g,%g) ", p.x(), p.y());
1971  if(seam) {
1972  p = ge->reparamOnFace(gf, range.low(), -1);
1973  printf("beg_alt (%g,%g) ", p.x(), p.y());
1974  }
1975  p = ge->reparamOnFace(gf, range.high(), 1);
1976  printf("end (%g,%g) ", p.x(), p.y());
1977  if(seam) {
1978  p = ge->reparamOnFace(gf, range.high(), -1);
1979  printf("end_alt (%g,%g) ", p.x(), p.y());
1980  }
1981  printf("\n");
1982  }
1983 #endif
1984 
1985  for (int initial_dir = 0; initial_dir < 2; initial_dir++)
1986  {
1987 
1988  if (coords.size())
1989  break; // we succeeded with initial_dir == 0
1990 
1991  for (std::size_t i = 0; i < signedEdges.size(); i++)
1992  {
1993  std::vector<SPoint2> p, p_alt, p_rev, p_alt_rev;
1994  std::vector<MVertex *> v, v_rev;
1995  GEdge *ge = signedEdges[i].getEdge();
1996  bool seam = ge->isSeam(gf);
1997  Range<double> range = ge->parBoundsOnFace(gf);
1998 
1999  // get parametric coordinates of curve nodes on the surface, with both
2000  // possible values for seams
2001  if (ge->getBeginVertex())
2002  {
2003  p.push_back(ge->reparamOnFace(gf, range.low(), 1));
2004  if (seam)
2005  p_alt.push_back(ge->reparamOnFace(gf, range.low(), -1));
2006  v.push_back(ge->getBeginVertex()->mesh_vertices[0]);
2007  }
2008  for (std::size_t j = 0; j < ge->mesh_vertices.size(); j++)
2009  {
2010  double u;
2011  ge->mesh_vertices[j]->getParameter(0, u);
2012  p.push_back(ge->reparamOnFace(gf, u, 1));
2013  if (seam)
2014  p_alt.push_back(ge->reparamOnFace(gf, u, -1));
2015  v.push_back(ge->mesh_vertices[j]);
2016  }
2017  if (ge->getEndVertex())
2018  {
2019  p.push_back(ge->reparamOnFace(gf, range.high(), 1));
2020  if (seam)
2021  p_alt.push_back(ge->reparamOnFace(gf, range.high(), -1));
2022  v.push_back(ge->getEndVertex()->mesh_vertices[0]);
2023  }
2024 
2025  // get the reverse mesh
2026  p_rev = p;
2027  std::reverse(p_rev.begin(), p_rev.end());
2028  if (seam)
2029  {
2030  p_alt_rev = p_alt;
2031  std::reverse(p_alt_rev.begin(), p_alt_rev.end());
2032  }
2033  v_rev = v;
2034  std::reverse(v_rev.begin(), v_rev.end());
2035 
2036  if (i == 0)
2037  {
2038  // choose which mesh to consider for the first curve in the loop
2039  if (initial_dir == 0)
2040  {
2041  if (seam && seam_the_first)
2042  coords = p_alt;
2043  else
2044  coords = p;
2045  verts = v;
2046  }
2047  else
2048  {
2049  if (seam && seam_the_first)
2050  coords = p_alt_rev;
2051  else
2052  coords = p_rev;
2053  verts = v_rev;
2054  }
2055  }
2056  else
2057  {
2058  // detect which mesh variant to use for the next curve by selecting the
2059  // mesh that starts with the node at the smallest distance, within the
2060  // prescribed tolerance
2061  double dist1 = coords.back().distance(p.front());
2062  double dist2 = coords.back().distance(p_rev.front());
2063  if (!seam)
2064  {
2065  if (dist1 < dist2 && dist1 < tol)
2066  {
2067  coords.pop_back();
2068  coords.insert(coords.end(), p.begin(), p.end());
2069  verts.pop_back();
2070  verts.insert(verts.end(), v.begin(), v.end());
2071  }
2072  else if (dist2 < dist1 && dist2 < tol)
2073  {
2074  coords.pop_back();
2075  coords.insert(coords.end(), p_rev.begin(), p_rev.end());
2076  verts.pop_back();
2077  verts.insert(verts.end(), v_rev.begin(), v_rev.end());
2078  }
2079  else
2080  {
2081  Msg::Debug("Distances (%g, %g) in parametric space larger than "
2082  "tolerance (%g) between end of curve %d and "
2083  "begining of curve %d...",
2084  dist1, dist2, tol, signedEdges[i - 1].getEdge()->tag(), ge->tag());
2085  if (initial_dir == 0)
2086  {
2087  Msg::Debug("... will try with alternate initial orientation");
2088  coords.clear();
2089  verts.clear();
2090  break;
2091  }
2092  else
2093  {
2094  Msg::Debug("... will try with larger tolerance");
2095  return false;
2096  }
2097  }
2098  }
2099  else
2100  {
2101  double dist3 = coords.back().distance(p_alt.front());
2102  double dist4 = coords.back().distance(p_alt_rev.front());
2103  if (dist1 < dist2 && dist1 < dist3 && dist1 < dist4 && dist1 < tol)
2104  {
2105  coords.pop_back();
2106  coords.insert(coords.end(), p.begin(), p.end());
2107  verts.pop_back();
2108  verts.insert(verts.end(), v.begin(), v.end());
2109  }
2110  else if (dist2 < dist1 && dist2 < dist3 && dist2 < dist4 && dist2 < tol)
2111  {
2112  coords.pop_back();
2113  coords.insert(coords.end(), p_rev.begin(), p_rev.end());
2114  verts.pop_back();
2115  verts.insert(verts.end(), v_rev.begin(), v_rev.end());
2116  }
2117  else if (dist3 < dist1 && dist3 < dist2 && dist3 < dist4 && dist3 < tol)
2118  {
2119  coords.pop_back();
2120  coords.insert(coords.end(), p_alt.begin(), p_alt.end());
2121  verts.pop_back();
2122  verts.insert(verts.end(), v.begin(), v.end());
2123  }
2124  else if (dist4 < dist1 && dist4 < dist2 && dist4 < dist3 && dist4 < tol)
2125  {
2126  coords.pop_back();
2127  coords.insert(coords.end(), p_alt_rev.begin(), p_alt_rev.end());
2128  verts.pop_back();
2129  verts.insert(verts.end(), v_rev.begin(), v_rev.end());
2130  }
2131  else
2132  {
2133  Msg::Debug("Distances (%g, %g, %g, %g) in parametric space larger "
2134  "than tolerance (%g) between end of curve %d and "
2135  "begining of seam curve %d...",
2136  dist1, dist2, dist3, dist4, tol, signedEdges[i - 1].getEdge()->tag(), ge->tag());
2137  if (initial_dir == 0)
2138  {
2139  Msg::Debug("... will try with alternate initial orientation");
2140  coords.clear();
2141  verts.clear();
2142  break;
2143  }
2144  else
2145  {
2146  Msg::Debug("... will try with larger tolerance");
2147  return false;
2148  }
2149  }
2150  }
2151  }
2152  }
2153  }
2154 
2155  if (verts.size() != coords.size())
2156  {
2157  Msg::Error("Wrong number of parametric coordinates for boundary nodes "
2158  "on surface %d",
2159  gf->tag());
2160  return false;
2161  }
2162  if (verts.empty())
2163  {
2164  Msg::Debug("No nodes in 1D mesh on surface %d", gf->tag());
2165  return true;
2166  }
2167  double dist = coords.back().distance(coords.front());
2168  if (dist < tol)
2169  {
2170  coords.pop_back();
2171  verts.pop_back();
2172  }
2173  else
2174  {
2175  Msg::Debug("Distance %g between first and last node in 1D mesh of surface "
2176  "%d exceeds tolerance %g",
2177  dist, gf->tag(), tol);
2178  return false;
2179  }
2180 
2181  std::map<BDS_Point *, MVertex *, PointLessThan> recoverMapLocal;
2182 
2183  for (std::size_t i = 0; i < verts.size(); i++)
2184  {
2185  MVertex *here = verts[i];
2186  GEntity *ge = here->onWhat();
2187  BDS_Point *pp = nullptr;
2188  if (ge->dim() == 0)
2189  {
2190  // point might already be part of another loop in the same surface
2191  for (auto it = recoverMap.begin(); it != recoverMap.end(); ++it)
2192  {
2193  if (it->second == here)
2194  {
2195  SPoint2 param(it->first->u, it->first->v);
2196  double dist = coords[i].distance(param);
2197  if (dist < tol)
2198  {
2199  pp = it->first;
2200  break;
2201  }
2202  }
2203  }
2204  }
2205  if (pp == nullptr)
2206  {
2207  double U, V;
2208  SPoint2 param = coords[i];
2209  U = param.x();
2210  V = param.y();
2211  pp = m->add_point(count + countTot, U, V, gf);
2212  if (ge->dim() == 0)
2213  {
2214  pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
2215  }
2216  else if (ge->dim() == 1)
2217  {
2218  double u;
2219  here->getParameter(0, u);
2220  pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
2221  }
2222  else
2223  {
2224  pp->lcBGM() = MAX_LC;
2225  }
2226  pp->lc() = pp->lcBGM();
2227  m->add_geom(ge->tag(), ge->dim());
2228  BDS_GeomEntity *g = m->get_geom(ge->tag(), ge->dim());
2229  pp->g = g;
2230  bbox += SPoint3(U, V, 0);
2231  }
2232  // printf("node %d coord %g %g\n", here->getNum(), pp->u, pp->v);
2233  result.push_back(pp);
2234  recoverMapLocal[pp] = here;
2235  count++;
2236  }
2237 
2238  Msg::Debug("Succeeded finding consecutive list of nodes on surface "
2239  "%d, with tolerance %g",
2240  gf->tag(), tol);
2241 
2242  // we're all set!
2243  recoverMap.insert(recoverMapLocal.begin(), recoverMapLocal.end());
2244 
2245  return true;
2246 }
2247 
2248 static GEdge *getGEdge(GFace *gf, MVertex *v1, MVertex *v2)
2249 {
2250  std::vector<GEdge *> e = gf->edges();
2251  for (size_t i = 0; i < e.size(); i++)
2252  {
2253  GEdge *ge = e[i];
2254  for (size_t j = 0; j < ge->lines.size(); j++)
2255  {
2256  MLine *l = ge->lines[j];
2257  if (l->getVertex(0) == v1 && l->getVertex(1) == v2)
2258  return ge;
2259  if (l->getVertex(1) == v1 && l->getVertex(0) == v2)
2260  return ge;
2261  }
2262  }
2263  return nullptr;
2264 }
2265 
2266 static bool meshGeneratorPeriodic(GFace *gf, int RECUR_ITER, bool repairSelfIntersecting1dMesh, bool debug = true)
2267 {
2268  if (CTX::instance()->debugSurface > 0 && gf->tag() != CTX::instance()->debugSurface)
2269  {
2270  gf->meshStatistics.status = GFace::DONE;
2271  return true;
2272  }
2273  if (CTX::instance()->debugSurface > 0)
2274  debug = true;
2275 
2276  // TEST !!!
2277  // PolyMesh * pm = GFaceInitialMesh(gf->tag(), 1);
2278 
2279  std::map<BDS_Point *, MVertex *, PointLessThan> recoverMap;
2280  std::multimap<MVertex *, BDS_Point *> recoverMultiMapInv;
2281 
2282  std::vector<SPoint2> true_boundary;
2283  trueBoundary(gf, true_boundary, debug);
2284 
2285  Range<double> rangeU = gf->parBounds(0);
2286  Range<double> rangeV = gf->parBounds(1);
2287 
2288  double du = rangeU.high() - rangeU.low();
2289  double dv = rangeV.high() - rangeV.low();
2290 
2291  const double LC2D = std::sqrt(du * du + dv * dv);
2292 
2293  // Buid a BDS_Mesh structure that is convenient for doing the actual meshing
2294  // procedure
2295  BDS_Mesh *m = new BDS_Mesh;
2296 
2297  std::vector<std::vector<BDS_Point *>> edgeLoops_BDS;
2298  SBoundingBox3d bbox;
2299  int nbPointsTotal = 0;
2300  {
2301  for (auto it = gf->edgeLoops.begin(); it != gf->edgeLoops.end(); it++)
2302  {
2303  std::vector<BDS_Point *> edgeLoop_BDS;
2304  int nbPointsLocal;
2305  const double fact[5] = {1.e-12, 1.e-9, 1.e-6, 1.e-3, 1.e-1};
2306  bool ok = false;
2307  for (int i = 0; i < 5; i++)
2308  {
2309  if (buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m, recoverMap, nbPointsLocal,
2310  nbPointsTotal, fact[i] * LC2D))
2311  {
2312  ok = true;
2313  break;
2314  }
2315  if (buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m, recoverMap, nbPointsLocal,
2316  nbPointsTotal, fact[i] * LC2D, true))
2317  {
2318  ok = true;
2319  break;
2320  }
2321  }
2322  if (!ok)
2323  {
2324  gf->meshStatistics.status = GFace::FAILED;
2325  Msg::Error("The 1D mesh seems not to be forming a closed loop");
2326  delete m;
2327  return false;
2328  }
2329  nbPointsTotal += nbPointsLocal;
2330  edgeLoops_BDS.push_back(edgeLoop_BDS);
2331  }
2332  }
2333 
2334  {
2335  auto it = recoverMap.begin();
2336  std::map<MVertex *, BDS_Point *> INV;
2337  for (; it != recoverMap.end(); ++it)
2338  {
2339  recoverMultiMapInv.insert(std::make_pair(it->second, it->first));
2340 
2341  auto it2 = INV.find(it->second);
2342  if (it2 != INV.end())
2343  {
2344  it->first->_periodicCounterpart = it2->second;
2345  it2->second->_periodicCounterpart = it->first;
2346  }
2347  INV[it->second] = it->first;
2348  }
2349  }
2350 
2351  if (nbPointsTotal < 3)
2352  {
2353  Msg::Warning("Mesh generation of surface %d skipped: only %d nodes on "
2354  "the boundary",
2355  gf->tag(), nbPointsTotal);
2356  gf->meshStatistics.status = GFace::DONE;
2357  delete m;
2358  return true;
2359  }
2360  else if (nbPointsTotal == 3)
2361  {
2362  MVertex *vv[3] = {nullptr, nullptr, nullptr};
2363  int i = 0;
2364  for (auto it = recoverMap.begin(); it != recoverMap.end(); it++)
2365  {
2366  vv[i++] = it->second;
2367  }
2368  if (vv[0] && vv[1] && vv[2])
2369  gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
2370  gf->meshStatistics.status = GFace::DONE;
2371  delete m;
2372  return true;
2373  }
2374 
2375  // Use a divide & conquer type algorithm to create a triangulation. We add to
2376  // the triangulation a box with 4 points that encloses the domain.
2377 
2378  std::vector<int> edgesEmbedded;
2379 
2380  {
2381  int count = 0;
2382 
2383  // Embedded Vertices
2384  // add embedded vertices
2385  std::vector<GVertex *> emb_vertx = gf->getEmbeddedVertices();
2386  auto itvx = emb_vertx.begin();
2387 
2388  std::map<MVertex *, std::set<BDS_Point *>> invertedRecoverMap;
2389  for (auto it = recoverMap.begin(); it != recoverMap.end(); it++)
2390  {
2391  invertedRecoverMap[it->second].insert(it->first);
2392  }
2393 
2394  int pNum = m->MAXPOINTNUMBER;
2395  nbPointsTotal += emb_vertx.size();
2396  {
2397  std::vector<GEdge *> emb_edges = gf->getEmbeddedEdges();
2398  std::set<MVertex *> vs;
2399  auto ite = emb_edges.begin();
2400  while (ite != emb_edges.end())
2401  {
2402  for (std::size_t i = 0; i < (*ite)->lines.size(); i++)
2403  {
2404  for (std::size_t j = 0; j < 2; j++)
2405  {
2406  MVertex *v = (*ite)->lines[i]->getVertex(j);
2407  if (invertedRecoverMap.find(v) == invertedRecoverMap.end() && vs.find(v) == vs.end())
2408  {
2409  vs.insert(v);
2410  }
2411  }
2412  }
2413  ++ite;
2414  }
2415  nbPointsTotal += vs.size();
2416  }
2417  DocRecord doc(nbPointsTotal + 4);
2418 
2419  while (itvx != emb_vertx.end())
2420  {
2421  MVertex *v = (*itvx)->mesh_vertices[0];
2422  double uv[2] = {0, 0};
2423  GPoint gp = gf->closestPoint(SPoint3(v->x(), v->y(), v->z()), uv);
2424  BDS_Point *pp = m->add_point(++pNum, gp.u(), gp.v(), gf);
2425  m->add_geom(-(*itvx)->tag(), 0);
2426  pp->g = m->get_geom(-(*itvx)->tag(), 0);
2427  pp->lcBGM() = BGM_MeshSize(*itvx, 0, 0, v->x(), v->y(), v->z());
2428  pp->lc() = pp->lcBGM();
2429  recoverMap[pp] = v;
2430  double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
2431  double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
2432  doc.points[count].where.h = pp->u + XX;
2433  doc.points[count].where.v = pp->v + YY;
2434  doc.points[count].adjacent = nullptr;
2435  doc.points[count].data = pp;
2436  count++;
2437  ++itvx;
2438  }
2439 
2440  std::vector<GEdge *> emb_edges = gf->getEmbeddedEdges();
2441  std::set<MVertex *> vs;
2442  std::map<MVertex *, BDS_Point *> facile;
2443  auto ite = emb_edges.begin();
2444  while (ite != emb_edges.end())
2445  {
2446  m->add_geom(-(*ite)->tag(), 1);
2447  for (std::size_t i = 0; i < (*ite)->lines.size(); i++)
2448  {
2449  for (std::size_t j = 0; j < 2; j++)
2450  {
2451  MVertex *v = (*ite)->lines[i]->getVertex(j);
2452  BDS_Point *pp = nullptr;
2453  const auto it = invertedRecoverMap.find(v);
2454  if (it != invertedRecoverMap.end())
2455  {
2456  if (it->second.size() > 1)
2457  {
2458  const GEdge *edge = (*ite);
2459  const Range<double> parBounds = edge->parBoundsOnFace(gf);
2460  GPoint firstPoint = edge->point(parBounds.low());
2461  GPoint lastPoint = edge->point(parBounds.high());
2462  double param;
2463  if (v->point().distance(SPoint3(firstPoint.x(), firstPoint.y(), firstPoint.z())) <
2464  v->point().distance(SPoint3(lastPoint.x(), lastPoint.y(), lastPoint.z())))
2465  {
2466  // Vertex lies on first point of edge
2467  param = parBounds.low();
2468  }
2469  else
2470  {
2471  // Vertex lies on last point of edge
2472  param = parBounds.high();
2473  }
2474  SPoint2 pointOnSurface = edge->reparamOnFace(gf, param, 1);
2475 
2476  const std::set<BDS_Point *> &possiblePoints = it->second;
2477  for (auto pntIt = possiblePoints.begin(); pntIt != possiblePoints.end(); ++pntIt)
2478  {
2479  if (pointOnSurface.distance(SPoint2((*pntIt)->u, (*pntIt)->v)) < 1e-10)
2480  {
2481  pp = (*pntIt);
2482  break;
2483  }
2484  }
2485  if (pp == nullptr)
2486  {
2487  Msg::Error("Embedded edge node %d is on the seam edge of "
2488  "surface %d and no appropriate point could be "
2489  "found!",
2490  v->getNum(), gf->tag());
2491  }
2492  }
2493  else
2494  {
2495  pp = *(it->second.begin());
2496  }
2497  facile[v] = pp;
2498  }
2499  if (pp == nullptr && vs.find(v) == vs.end())
2500  {
2501  vs.insert(v);
2502  double uv[2] = {0, 0};
2503  GPoint gp = gf->closestPoint(SPoint3(v->x(), v->y(), v->z()), uv);
2504  BDS_Point *pp = m->add_point(++pNum, gp.u(), gp.v(), gf);
2505  pp->g = m->get_geom(-(*ite)->tag(), 1);
2506  if (v->onWhat()->dim() == 0)
2507  pp->lcBGM() = BGM_MeshSize(v->onWhat(), 0, 0, v->x(), v->y(), v->z());
2508  else
2509  {
2510  double uu;
2511  v->getParameter(0, uu);
2512  pp->lcBGM() = BGM_MeshSize(*ite, uu, 0, v->x(), v->y(), v->z());
2513  }
2514  pp->lc() = pp->lcBGM();
2515  recoverMap[pp] = v;
2516  facile[v] = pp;
2517  double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
2518  double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
2519  doc.points[count].where.h = pp->u + XX;
2520  doc.points[count].where.v = pp->v + YY;
2521  doc.points[count].adjacent = nullptr;
2522  doc.points[count].data = pp;
2523  count++;
2524  }
2525  }
2526  }
2527  for (std::size_t i = 0; i < (*ite)->lines.size(); i++)
2528  {
2529  BDS_Point *p0 = facile[(*ite)->lines[i]->getVertex(0)];
2530  BDS_Point *p1 = facile[(*ite)->lines[i]->getVertex(1)];
2531  edgesEmbedded.push_back(p0->iD);
2532  edgesEmbedded.push_back(p1->iD);
2533  }
2534  ++ite;
2535  }
2536 
2537  for (std::size_t i = 0; i < edgeLoops_BDS.size(); i++)
2538  {
2539  std::vector<BDS_Point *> &edgeLoop_BDS = edgeLoops_BDS[i];
2540  for (std::size_t j = 0; j < edgeLoop_BDS.size(); j++)
2541  {
2542  BDS_Point *pp = edgeLoop_BDS[j];
2543  double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
2544  double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX;
2545  doc.points[count].where.h = pp->u + XX;
2546  doc.points[count].where.v = pp->v + YY;
2547  doc.points[count].adjacent = nullptr;
2548  doc.points[count].data = pp;
2549  count++;
2550  }
2551  }
2552 
2553  // Increase the size of the bounding box, add 4 points that enclose
2554  // the domain, use negative number to distinguish those fake
2555  // vertices
2556 
2557  if (du / dv < 1200 && dv / du < 1200)
2558  {
2559  // Fix a bug here if the size of the box is zero
2560  bbox.makeCube();
2561  }
2562 
2563  bbox *= 3.5;
2564  GPoint bb[4] = {GPoint(bbox.min().x(), bbox.min().y(), 0), GPoint(bbox.min().x(), bbox.max().y(), 0),
2565  GPoint(bbox.max().x(), bbox.min().y(), 0), GPoint(bbox.max().x(), bbox.max().y(), 0)};
2566  for (int ip = 0; ip < 4; ip++)
2567  {
2568  BDS_Point *pp = m->add_point(-ip - 1, bb[ip].x(), bb[ip].y(), gf);
2569  m->add_geom(gf->tag(), 2);
2570  BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
2571  pp->g = g;
2572  doc.points[nbPointsTotal + ip].where.h = bb[ip].x();
2573  doc.points[nbPointsTotal + ip].where.v = bb[ip].y();
2574  doc.points[nbPointsTotal + ip].adjacent = nullptr;
2575  doc.points[nbPointsTotal + ip].data = pp;
2576  }
2577 
2578  // Use "fast" inhouse recursive algo to generate the triangulation
2579  // At this stage the triangulation is not what we need
2580  // -) It does not necessary recover the boundaries
2581  // -) It contains triangles outside the domain (the first edge
2582  // loop is the outer one)
2583  Msg::Debug("Meshing of the convex hull (%d nodes)", nbPointsTotal);
2584 
2585  try
2586  {
2587  doc.MakeMeshWithPoints();
2588  }
2589  catch (std::runtime_error &e)
2590  {
2591  Msg::Error("%s", e.what());
2592  }
2593 
2594  for (int i = 0; i < doc.numTriangles; i++)
2595  {
2596  int a = doc.triangles[i].a;
2597  int b = doc.triangles[i].b;
2598  int c = doc.triangles[i].c;
2599  int n = doc.numPoints;
2600  if (a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n)
2601  {
2602  Msg::Warning("Skipping bad triangle %d", i);
2603  continue;
2604  }
2605  BDS_Point *p1 = (BDS_Point *)doc.points[doc.triangles[i].a].data;
2606  BDS_Point *p2 = (BDS_Point *)doc.points[doc.triangles[i].b].data;
2607  BDS_Point *p3 = (BDS_Point *)doc.points[doc.triangles[i].c].data;
2608  m->add_triangle(p1->iD, p2->iD, p3->iD);
2609  }
2610  }
2611 
2612  // Recover the boundary edges and compute characteristic lenghts using mesh
2613  // edge spacing
2614  BDS_GeomEntity CLASS_F(1, 2);
2615  BDS_GeomEntity CLASS_E(1, 1);
2616  BDS_GeomEntity CLASS_EXTERIOR(3, 2);
2617 
2618  if (debug)
2619  {
2620  char name[245];
2621  sprintf(name, "surface%d-initial-real.pos", gf->tag());
2622  outputScalarField(m->triangles, name, 0, gf);
2623  sprintf(name, "surface%d-initial-param.pos", gf->tag());
2624  outputScalarField(m->triangles, name, 1, gf);
2625  }
2626 
2627  bool _fatallyFailed;
2628 
2629  for (std::size_t i = 0; i < edgesEmbedded.size() / 2; i++)
2630  {
2631  BDS_Edge *e = m->recover_edge(edgesEmbedded[2 * i], edgesEmbedded[2 * i + 1], _fatallyFailed);
2632  if (!e)
2633  {
2634  Msg::Error("Impossible to recover the edge %d %d", edgesEmbedded[2 * i], edgesEmbedded[2 * i + 1]);
2635  gf->meshStatistics.status = GFace::FAILED;
2636  delete m;
2637  return false;
2638  }
2639  else
2640  e->g = &CLASS_E;
2641  }
2642 
2643  std::set<EdgeToRecover> edgesNotRecovered;
2644 
2645  bool doItAgain = gf->meshStatistics.refineAllEdges;
2646  for (std::size_t i = 0; i < edgeLoops_BDS.size(); i++)
2647  {
2648  std::vector<BDS_Point *> &edgeLoop_BDS = edgeLoops_BDS[i];
2649  for (std::size_t j = 0; j < edgeLoop_BDS.size(); j++)
2650  {
2651  int num1 = edgeLoop_BDS[j]->iD;
2652  int num2 = edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD;
2653  BDS_Edge *e = m->find_edge(num1, num2);
2654  if (!e)
2655  {
2656  // printf("recovering %d %d\n",num1,num2);
2657  e = m->recover_edge(num1, num2, _fatallyFailed);
2658  BDS_Point *p1 = m->find_point(num1);
2659  BDS_Point *p2 = m->find_point(num2);
2660  MVertex *v1 = recoverMap[p1];
2661  MVertex *v2 = recoverMap[p2];
2662  GEdge *ge = getGEdge(gf, v1, v2);
2663  // if (!ge){
2664  // Msg::Error("cannot find GEdge with mesh edge %d %d (%d %d)\n",
2665  // num1, num2, v1->getNum(), v2->getNum());
2666  // }
2667  if (ge)
2668  edgesNotRecovered.insert(EdgeToRecover(num1, num2, ge));
2669  if (!e)
2670  {
2671  // what is before does not seem to work properly
2672  // Msg::Warning("ITER %d Impossible to recover the edge %d %d",
2673  // RECUR_ITER, num1, num2);
2674  // Msg::Warning("Will split model edge %d and try again", ge->tag());
2675  doItAgain = true;
2676  // gf->meshStatistics.status = GFace::FAILED;
2677  // delete m;
2678  // return false;
2679  }
2680  else
2681  {
2682  // Msg::Warning("ITER %d edge %d %d RECOVERED",RECUR_ITER, num1,num2);
2683  e->g = &CLASS_E;
2684  }
2685  }
2686  else
2687  e->g = &CLASS_E;
2688  }
2689  }
2690 
2691  if (doItAgain)
2692  {
2693  // this block is not thread safe. 2D mesh will be serialized for surfaces
2694  // that have their 1D mesh self-intersect
2695  if (Msg::GetNumThreads() != 1)
2696  {
2697  gf->meshStatistics.status = GFace::PENDING;
2698  delete m;
2699  Msg::Info("Surface %d has self-intersections in its 1D mesh: serializing "
2700  "this one",
2701  gf->tag());
2702  return true;
2703  }
2704 
2705  if (!gf->meshStatistics.refineAllEdges)
2706  {
2707  std::ostringstream sstream;
2708  for (auto itr = edgesNotRecovered.begin(); itr != edgesNotRecovered.end(); ++itr)
2709  sstream << " " << itr->ge->tag() << "("
2710  << (itr->ge->getBeginVertex() ? itr->ge->getBeginVertex()->tag() : -1) << ","
2711  << (itr->ge->getEndVertex() ? itr->ge->getEndVertex()->tag() : -1) << ")";
2712  Msg::Info(":-( There are %d intersections in the 1D mesh (curves%s)", edgesNotRecovered.size(),
2713  sstream.str().c_str());
2714  }
2715  if (repairSelfIntersecting1dMesh)
2716  {
2717  Msg::Info("8-| Splitting those edges and trying again");
2718  if (gf->meshStatistics.refineAllEdges)
2719  {
2720  std::vector<GEdge *> eds = gf->edges();
2721  edgesNotRecovered.clear();
2722  for (size_t i = 0; i < eds.size(); i++)
2723  {
2724  const std::size_t NN = eds[i]->lines.size() ? 1 : 0;
2725  for (size_t j = 0; j < NN; j++)
2726  {
2727  MVertex *v1 = eds[i]->lines[j]->getVertex(0);
2728  MVertex *v2 = eds[i]->lines[j]->getVertex(1);
2729  auto itp1 = recoverMultiMapInv.lower_bound(v1);
2730  auto itp2 = recoverMultiMapInv.lower_bound(v2);
2731  if (itp1 != recoverMultiMapInv.end() && itp2 != recoverMultiMapInv.end())
2732  edgesNotRecovered.insert(EdgeToRecover(itp1->second->iD, itp2->second->iD, eds[i]));
2733  }
2734  }
2735  }
2736  remeshUnrecoveredEdges(recoverMultiMapInv, edgesNotRecovered, gf->meshStatistics.refineAllEdges);
2737  gf->meshStatistics.refineAllEdges = false;
2738  }
2739  // delete the mesh
2740  // getchar();
2741  if (debug)
2742  {
2743  char name[245];
2744  sprintf(name, "surface%d-initial-real-afterITER%d.pos", gf->tag(), RECUR_ITER);
2745  outputScalarField(m->triangles, name, 0, gf);
2746  sprintf(name, "surface%d-initial-param-afterITER%d.pos", gf->tag(), RECUR_ITER);
2747  outputScalarField(m->triangles, name, 1, gf);
2748  }
2749  delete m;
2750  if (RECUR_ITER < 10)
2751  {
2752  return meshGeneratorPeriodic(gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, debug);
2753  }
2754  return false;
2755  }
2756 
2757  if (RECUR_ITER > 0)
2758  Msg::Info(":-) All edges recovered after %d iteration%s", RECUR_ITER, (RECUR_ITER > 1) ? "s" : "");
2759 
2760  // look for a triangle that has a negative node and recursively tag all
2761  // exterior triangles
2762  {
2763  auto itt = m->triangles.begin();
2764  while (itt != m->triangles.end())
2765  {
2766  (*itt)->g = nullptr;
2767  ++itt;
2768  }
2769  itt = m->triangles.begin();
2770  while (itt != m->triangles.end())
2771  {
2772  BDS_Face *t = *itt;
2773  if (t->deleted)
2774  {
2775  // If triangle is deleted, it won't have the correct neighbours
2776  // to tag recursively
2777  ++itt;
2778  continue;
2779  }
2780  BDS_Point *n[4];
2781  if (t->getNodes(n))
2782  {
2783  if (n[0]->iD < 0 || n[1]->iD < 0 || n[2]->iD < 0)
2784  {
2785  recur_tag(t, &CLASS_EXTERIOR);
2786  break;
2787  }
2788  }
2789  ++itt;
2790  }
2791  }
2792 
2793  // now find an edge that has belongs to one of the exterior triangles
2794  {
2795  auto ite = m->edges.begin();
2796  while (ite != m->edges.end())
2797  {
2798  BDS_Edge *edge = *ite;
2799  if (edge->g && edge->numfaces() == 2)
2800  {
2801  if (edge->faces(0)->g == &CLASS_EXTERIOR)
2802  {
2803  recur_tag(edge->faces(1), &CLASS_F);
2804  break;
2805  }
2806  else if (edge->faces(1)->g == &CLASS_EXTERIOR)
2807  {
2808  recur_tag(edge->faces(0), &CLASS_F);
2809  break;
2810  }
2811  }
2812  ++ite;
2813  }
2814  auto itt = m->triangles.begin();
2815  while (itt != m->triangles.end())
2816  {
2817  if ((*itt)->g == &CLASS_EXTERIOR)
2818  (*itt)->g = nullptr;
2819  ++itt;
2820  }
2821  }
2822 
2823  // delete useless stuff
2824  {
2825  auto itt = m->triangles.begin();
2826  while (itt != m->triangles.end())
2827  {
2828  BDS_Face *t = *itt;
2829  if (!t->g)
2830  {
2831  m->del_face(t);
2832  }
2833  ++itt;
2834  }
2835  }
2836 
2837  m->cleanup();
2838 
2839  {
2840  auto ite = m->edges.begin();
2841  while (ite != m->edges.end())
2842  {
2843  BDS_Edge *edge = *ite;
2844  if (edge->numfaces() == 0)
2845  m->del_edge(edge);
2846  else
2847  {
2848  if (!edge->g)
2849  edge->g = &CLASS_F;
2850  if (!edge->p1->g || edge->p1->g->classif_degree > edge->g->classif_degree)
2851  edge->p1->g = edge->g;
2852  if (!edge->p2->g || edge->p2->g->classif_degree > edge->g->classif_degree)
2853  edge->p2->g = edge->g;
2854  }
2855  ++ite;
2856  }
2857  }
2858  m->cleanup();
2859  m->del_point(m->find_point(-1));
2860  m->del_point(m->find_point(-2));
2861  m->del_point(m->find_point(-3));
2862  m->del_point(m->find_point(-4));
2863 
2864  if (debug)
2865  {
2866  char name[245];
2867  sprintf(name, "surface%d-recovered-real.pos", gf->tag());
2868  outputScalarField(m->triangles, name, 0, gf);
2869  sprintf(name, "surface%d-recovered-param.pos", gf->tag());
2870  outputScalarField(m->triangles, name, 1, gf);
2871  }
2872 
2873  if (algoDelaunay2D(gf))
2874  {
2875  // Call this function to untangle elements in Cartesian space
2876  Msg::Debug("Delaunizing the initial mesh");
2877  int nb_swap;
2878  delaunayizeBDS(gf, *m, nb_swap);
2879  }
2880  else
2881  {
2882  // tag points that are degenerated
2883  modifyInitialMeshToRemoveDegeneracies(gf, *m, &recoverMap);
2884 
2885  Msg::Debug("Delaunizing the initial mesh");
2886  int nb_swap;
2887  delaunayizeBDS(gf, *m, nb_swap);
2888 
2889  refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true, nullptr, &recoverMap, &true_boundary);
2890  if (debug)
2891  {
2892  char name[245];
2893  sprintf(name, "surface%d-phase1-real.pos", gf->tag());
2894  outputScalarField(m->triangles, name, 0, gf);
2895  sprintf(name, "surface%d-phase1-param.pos", gf->tag());
2896  outputScalarField(m->triangles, name, 1, gf);
2897  }
2898  if (debug)
2899  {
2900  char name[245];
2901  sprintf(name, "surface%d-phase2-real.pos", gf->tag());
2902  outputScalarField(m->triangles, name, 0, gf);
2903  sprintf(name, "surface%d-phase2-param.pos", gf->tag());
2904  outputScalarField(m->triangles, name, 1, gf);
2905  }
2906  refineMeshBDS(gf, *m, -CTX::instance()->mesh.refineSteps, false, nullptr, &recoverMap, &true_boundary);
2907  if (debug)
2908  {
2909  char name[245];
2910  sprintf(name, "surface%d-phase3-real.pos", gf->tag());
2911  outputScalarField(m->triangles, name, 0);
2912  sprintf(name, "surface%d-phase3-param.pos", gf->tag());
2913  outputScalarField(m->triangles, name, 1, gf);
2914  }
2915  if (debug)
2916  {
2917  char name[245];
2918  sprintf(name, "surface%d-phase4-real.pos", gf->tag());
2919  outputScalarField(m->triangles, name, 0);
2920  sprintf(name, "surface%d-phase4-param.pos", gf->tag());
2921  outputScalarField(m->triangles, name, 1, gf);
2922  }
2923 
2924  if (gf->meshStatistics.status == GFace::FAILED)
2925  {
2926  // splitall
2927  gf->meshStatistics.status = GFace::PENDING;
2928  gf->meshStatistics.refineAllEdges = true;
2929  delete m;
2930  Msg::Info("Serializing surface %d and refining all its bounding edges", gf->tag());
2931  return true;
2932  }
2933  else
2934  {
2935  gf->meshStatistics.status = GFace::DONE;
2936  }
2937 
2938  if (debug)
2939  {
2940  char name[245];
2941  sprintf(name, "surface%d-just-real.pos", gf->tag());
2942  outputScalarField(m->triangles, name, 0, gf);
2943  }
2944  }
2945 
2946  // This is a structure that we need only for periodic cases. We will duplicate
2947  // the vertices (MVertex) that are on seams
2948 
2949  std::map<MVertex *, MVertex *> equivalence;
2950  std::map<MVertex *, SPoint2> parametricCoordinates;
2951  if (algoDelaunay2D(gf))
2952  {
2953  std::map<MVertex *, BDS_Point *> invertMap;
2954  auto it = recoverMap.begin();
2955  while (it != recoverMap.end())
2956  {
2957  // we have twice vertex MVertex with 2 different coordinates
2958  MVertex *mv1 = it->second;
2959  BDS_Point *bds = it->first;
2960  auto invIt = invertMap.find(mv1);
2961  if (invIt != invertMap.end())
2962  {
2963  // create a new "fake" vertex that will be destroyed afterwards
2964  MVertex *mv2 = nullptr;
2965  if (mv1->onWhat()->dim() == 1)
2966  {
2967  double t;
2968  mv1->getParameter(0, t);
2969  mv2 = new MEdgeVertex(mv1->x(), mv1->y(), mv1->z(), mv1->onWhat(), t, 0,
2970  ((MEdgeVertex *)mv1)->getLc());
2971  }
2972  else if (mv1->onWhat()->dim() == 0)
2973  {
2974  mv2 = new MVertex(mv1->x(), mv1->y(), mv1->z(), mv1->onWhat());
2975  }
2976  else
2977  Msg::Error("Could not reconstruct seam");
2978  if (mv2)
2979  {
2980  it->second = mv2;
2981  equivalence[mv2] = mv1;
2982  parametricCoordinates[mv2] = SPoint2(bds->u, bds->v);
2983  invertMap[mv2] = bds;
2984  }
2985  }
2986  else
2987  {
2988  parametricCoordinates[mv1] = SPoint2(bds->u, bds->v);
2989  invertMap[mv1] = bds;
2990  }
2991  ++it;
2992  }
2993  }
2994 
2995  // Msg::Info("%d points that are duplicated for Delaunay meshing",
2996  // equivalence.size());
2997 
2998  // fill the small gmsh structures
2999  BDS2GMSH(m, gf, recoverMap);
3000 
3001  if (debug)
3002  {
3003  char name[245];
3004  sprintf(name, "surface%d-final-real.pos", gf->tag());
3005  outputScalarField(m->triangles, name, 0, gf);
3006  sprintf(name, "surface%d-final-param.pos", gf->tag());
3007  outputScalarField(m->triangles, name, 1, gf);
3008  }
3009 
3010  bool infty = false;
3011  if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD)
3012  {
3013  infty = true;
3014  buildBackgroundMesh(gf, CTX::instance()->mesh.crossFieldClosestPoint, &equivalence, &parametricCoordinates);
3015  }
3016  else if (gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS || gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
3017  {
3018  infty = true;
3019  /* New version of PACK / QUADQS use a different background mesh */
3020  // buildBackgroundMesh(gf, false, &equivalence, &parametricCoordinates);
3021  }
3022 
3023  bool onlyInitialMesh = (gf->getMeshingAlgo() == ALGO_2D_INITIAL_ONLY);
3024 
3025  if (!onlyInitialMesh)
3026  {
3027  // boundary layer
3028  std::vector<MQuadrangle *> blQuads;
3029  std::vector<MTriangle *> blTris;
3030  std::set<MVertex *> verts;
3031  modifyInitialMeshForBoundaryLayers(gf, blQuads, blTris, verts, debug);
3032  gf->quadrangles.insert(gf->quadrangles.begin(), blQuads.begin(), blQuads.end());
3033  gf->triangles.insert(gf->triangles.begin(), blTris.begin(), blTris.end());
3034  gf->mesh_vertices.insert(gf->mesh_vertices.begin(), verts.begin(), verts.end());
3035  }
3036 
3037  if (algoDelaunay2D(gf) && !onlyInitialMesh)
3038  {
3039  if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL)
3040  bowyerWatsonFrontal(gf, &equivalence, &parametricCoordinates, &true_boundary);
3041  else if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD)
3042  bowyerWatsonFrontalLayers(gf, true, &equivalence, &parametricCoordinates);
3043  else if (gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
3044  bowyerWatsonParallelograms(gf, &equivalence, &parametricCoordinates);
3045  else if (gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
3046  {
3047  Msg::Error("ALGO_2D_PACK_PRLGRMS_CSTR deprecated");
3048  // bowyerWatsonParallelogramsConstrained(
3049  // gf, gf->constr_vertices, &equivalence, &parametricCoordinates);
3050  }
3051  else if (gf->getMeshingAlgo() == ALGO_2D_DELAUNAY || gf->getMeshingAlgo() == ALGO_2D_AUTO)
3052  bowyerWatson(gf, 1000000000, &equivalence, &parametricCoordinates);
3053  else
3054  meshGFaceBamg(gf);
3055  if (!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
3056  laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, infty);
3057  }
3058 
3059  if (debug)
3060  {
3061  char name[256];
3062  sprintf(name, "real%d.pos", gf->tag());
3063  outputScalarField(m->triangles, name, 0, gf);
3064  sprintf(name, "param%d.pos", gf->tag());
3065  outputScalarField(m->triangles, name, 1, gf);
3066  }
3067 
3068  // remove fake duplicate nodes on seams
3069  for (auto eq : equivalence)
3070  delete eq.first;
3071 
3072  // delete the mesh
3073  delete m;
3074 
3075  if ((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
3076  (CTX::instance()->mesh.algoRecombine <= 1 || CTX::instance()->mesh.algoRecombine == 4))
3077  {
3078 
3079  if (CTX::instance()->mesh.algoRecombine == 4)
3080  {
3081  meshGFaceQuadrangulateBipartiteLabelling(gf->tag());
3082  }
3083  else
3084  {
3085  bool blossom = (CTX::instance()->mesh.algoRecombine == 1);
3086  int topo = CTX::instance()->mesh.recombineOptimizeTopology;
3087  int repos = CTX::instance()->mesh.recombineNodeRepositioning;
3088  if (CTX::instance()->mesh.algo2d == ALGO_2D_QUAD_QUASI_STRUCT)
3089  repos = false;
3090  double minqual = CTX::instance()->mesh.recombineMinimumQuality;
3091  recombineIntoQuads(gf, blossom, topo, repos, minqual);
3092  }
3093  }
3094 
3095  computeElementShapes(gf, gf->meshStatistics.worst_element_shape, gf->meshStatistics.average_element_shape,
3096  gf->meshStatistics.best_element_shape, gf->meshStatistics.nbTriangle,
3097  gf->meshStatistics.nbGoodQuality);
3098  gf->meshStatistics.status = GFace::DONE;
3099 
3100  // Remove unused vertices, generated e.g. during background mesh
3101  deleteUnusedVertices(gf);
3102 
3103  return true;
3104 }
3105 
3106 void deMeshGFace::operator()(GFace *gf)
3107 {
3108  if (gf->isFullyDiscrete())
3109  return;
3110  gf->deleteMesh();
3111  gf->meshStatistics.status = GFace::PENDING;
3112  gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
3113 }
3114 
3115 static double TRIANGLE_VALIDITY(GFace *gf, MTriangle *t)
3116 {
3117  SPoint2 p1, p2, p3;
3118  reparamMeshVertexOnFace(t->getVertex(0), gf, p1);
3119  reparamMeshVertexOnFace(t->getVertex(1), gf, p2);
3120  reparamMeshVertexOnFace(t->getVertex(2), gf, p3);
3121  SVector3 N1 = gf->normal(p1);
3122  SVector3 N2 = gf->normal(p2);
3123  SVector3 N3 = gf->normal(p3);
3124  SVector3 N = N1 + N2 + N3;
3125  N.normalize();
3126  SVector3 d1(t->getVertex(1)->x() - t->getVertex(0)->x(), t->getVertex(1)->y() - t->getVertex(0)->y(),
3127  t->getVertex(1)->z() - t->getVertex(0)->z());
3128  SVector3 d2(t->getVertex(2)->x() - t->getVertex(0)->x(), t->getVertex(2)->y() - t->getVertex(0)->y(),
3129  t->getVertex(2)->z() - t->getVertex(0)->z());
3130  SVector3 c = crossprod(d1, d2);
3131  return dot(N, c);
3132 }
3133 
3134 static bool isMeshValid(GFace *gf)
3135 {
3136  size_t invalid = 0;
3137  for (size_t i = 0; i < gf->triangles.size(); i++)
3138  {
3139  double v = TRIANGLE_VALIDITY(gf, gf->triangles[i]);
3140  if (v < 0)
3141  invalid++;
3142  }
3143  if (invalid == 0 || invalid == gf->triangles.size())
3144  return true;
3145 
3146  return false;
3147 }
3148 
3149 // for debugging, change value from -1 to -100;
3150 int debugSurface = -1; //-100;
3151 
3152 void meshGFace::operator()(GFace *gf, bool print)
3153 {
3154  gf->model()->setCurrentMeshEntity(gf);
3155 
3156  if (gf->meshAttributes.method == MESH_NONE)
3157  return;
3158  if (CTX::instance()->mesh.meshOnlyVisible && !gf->getVisibility())
3159  return;
3160  if (CTX::instance()->mesh.meshOnlyEmpty && gf->getNumMeshElements())
3161  return;
3162 
3163  // destroy the mesh if it exists
3164  deMeshGFace dem;
3165  dem(gf);
3166 
3167  if (MeshTransfiniteSurface(gf))
3168  return;
3169  if (MeshExtrudedSurface(gf))
3170  return;
3171  if (gf->getMeshMaster() != gf)
3172  {
3173  GFace *gff = dynamic_cast<GFace *>(gf->getMeshMaster());
3174  if (gff)
3175  {
3176  if (gff->meshStatistics.status != GFace::DONE)
3177  {
3178  gf->meshStatistics.status = GFace::PENDING;
3179  return;
3180  }
3181  Msg::Info("Meshing surface %d (%s) as a copy of surface %d", gf->tag(), gf->getTypeString().c_str(),
3182  gf->getMeshMaster()->tag());
3183  copyMesh(gff, gf);
3184  gf->meshStatistics.status = GFace::DONE;
3185  return;
3186  }
3187  else
3188  Msg::Warning("Unknown mesh master surface %d", gf->getMeshMaster()->tag());
3189  }
3190 
3191  /* The ALGO_2D_QUAD_QUASI_STRUCT is using ALGO_2D_PACK_PRLGRMS
3192  * to generate a initial quad-dominant mesh */
3193  bool quadqs = false;
3194  if (CTX::instance()->mesh.algo2d == ALGO_2D_QUAD_QUASI_STRUCT)
3195  {
3196  quadqs = true;
3197  CTX::instance()->mesh.algo2d = ALGO_2D_PACK_PRLGRMS;
3198  }
3199 
3200  const char *algo = "Unknown";
3201 
3202  switch (gf->getMeshingAlgo())
3203  {
3204  case ALGO_2D_MESHADAPT:
3205  algo = "MeshAdapt";
3206  break;
3207  case ALGO_2D_AUTO:
3208  algo = (gf->geomType() == GEntity::Plane) ? "Delaunay" : "MeshAdapt";
3209  break;
3210  case ALGO_2D_INITIAL_ONLY:
3211  algo = "Initial Mesh Only";
3212  break;
3213  case ALGO_2D_DELAUNAY:
3214  algo = "Delaunay";
3215  break;
3216  case ALGO_2D_FRONTAL:
3217  algo = "Frontal-Delaunay";
3218  break;
3219  case ALGO_2D_BAMG:
3220  algo = "Bamg";
3221  break;
3222  case ALGO_2D_FRONTAL_QUAD:
3223  algo = "Frontal-Delaunay for Quads";
3224  break;
3225  case ALGO_2D_PACK_PRLGRMS:
3226  algo = "Packing of Parallelograms";
3227  break;
3228  case ALGO_2D_PACK_PRLGRMS_CSTR:
3229  algo = "Packing of Parallelograms Constrained";
3230  break;
3231  }
3232 
3233  if (print)
3234  Msg::Info("Meshing surface %d (%s, %s)", gf->tag(), gf->getTypeString().c_str(), algo);
3235 
3236  bool singularEdges = false;
3237  auto ite = gf->edges().begin();
3238  while (ite != gf->edges().end())
3239  {
3240  if ((*ite)->isSeam(gf))
3241  singularEdges = true;
3242  ite++;
3243  }
3244 
3245  bool periodic =
3246  (gf->getNativeType() != GEntity::GmshModel) && (gf->periodic(0) || gf->periodic(1) || singularEdges);
3247 
3248  quadMeshRemoveHalfOfOneDMesh halfmesh(gf, periodic);
3249 
3250  if (periodic)
3251  {
3252  if (!meshGeneratorPeriodic(gf, 0, repairSelfIntersecting1dMesh, debugSurface >= 0 || debugSurface == -100))
3253  {
3254  Msg::Error("Impossible to mesh periodic surface %d", gf->tag());
3255  gf->meshStatistics.status = GFace::FAILED;
3256  }
3257  }
3258  else
3259  {
3260  meshGenerator(gf, 0, repairSelfIntersecting1dMesh, (gf->getMeshingAlgo() == ALGO_2D_INITIAL_ONLY) ? 1 : 0,
3261  (debugSurface >= 0 || debugSurface == -100), NULL);
3262  }
3263 
3264  Msg::Debug("Type %d %d triangles generated, %d internal nodes", gf->geomType(), gf->triangles.size(),
3265  gf->mesh_vertices.size());
3266 
3267  halfmesh.finish();
3268 
3269  if (gf->getNumMeshElements() == 0 && gf->meshStatistics.status == GFace::DONE)
3270  {
3271  Msg::Warning("Surface %d consists of no elements", gf->tag());
3272  }
3273 
3274  // test validity for non-Gmsh models (currently we cannot reliably evaluate
3275  // the normal on the boundary of surfaces with the Gmsh kernel)
3276  if (CTX::instance()->mesh.algoSwitchOnFailure && gf->getNativeType() != GEntity::GmshModel &&
3277  gf->geomType() != GEntity::Plane && algoDelaunay2D(gf) && !isMeshValid(gf))
3278  {
3279  Msg::Debug("Delaunay-based mesher failed on surface %d -> moving to MeshAdapt", gf->tag());
3280  deMeshGFace killer;
3281  killer(gf);
3282  gf->setMeshingAlgo(1);
3283  (*this)(gf, print);
3284  gf->unsetMeshingAlgo();
3285  }
3286 
3287  if (quadqs)
3288  CTX::instance()->mesh.algo2d = ALGO_2D_QUAD_QUASI_STRUCT;
3289 }
3290 
3291 static bool getGFaceNormalFromVert(GFace *gf, MElement *el, SVector3 &nf)
3292 {
3293  bool found = false;
3294  for (std::size_t iElV = 0; iElV < el->getNumVertices(); iElV++)
3295  {
3296  MVertex *v = el->getVertex(iElV);
3297  SPoint2 param;
3298  if (v->onWhat() == gf && v->getParameter(0, param[0]) && v->getParameter(1, param[1]))
3299  {
3300  nf = gf->normal(param);
3301  found = true;
3302  break;
3303  }
3304  }
3305  return found;
3306 }
3307 
3308 static bool getGFaceNormalFromBary(GFace *gf, MElement *el, SVector3 &nf)
3309 {
3310  SPoint2 param(0., 0.);
3311  bool ok = true;
3312  for (std::size_t j = 0; j < el->getNumVertices(); j++)
3313  {
3314  SPoint2 p;
3315  // FIXME: use inexact reparam because some vertices might not be exactly on
3316  // the surface after the 3D Delaunay
3317  ok = reparamMeshVertexOnFace(el->getVertex(j), gf, p, false);
3318  if (!ok)
3319  break;
3320  param += p;
3321  }
3322  if (ok)
3323  {
3324  param *= 1.0 / el->getNumVertices();
3325  nf = gf->normal(param);
3326  }
3327  return ok;
3328 }
3329 
3330 static void getGFaceOrientation(GFace *gf, BoundaryLayerColumns *blc, bool existBL, bool fromVert, int &orientNonBL,
3331  int &orientBL)
3332 {
3333  for (std::size_t iEl = 0; iEl < gf->getNumMeshElements(); iEl++)
3334  {
3335  MElement *e = gf->getMeshElement(iEl);
3336  const bool isBLEl = existBL && (blc->_toFirst.find(e) != blc->_toFirst.end());
3337  SVector3 nf;
3338  // Check only if orientation of BL/non-BL el. not already known
3339  if ((!isBLEl && orientNonBL == 0) || (isBLEl && orientBL == 0))
3340  {
3341  const bool found = fromVert ? getGFaceNormalFromVert(gf, e, nf) : getGFaceNormalFromBary(gf, e, nf);
3342  if (found)
3343  {
3344  SVector3 ne = e->getFace(0).normal();
3345  const int orient = (dot(ne, nf) > 0.) ? 1 : -1;
3346  if (isBLEl)
3347  orientBL = orient;
3348  else
3349  orientNonBL = orient;
3350  }
3351  }
3352  // Stop when orientation found for non-BL and BL el.
3353  if ((orientNonBL != 0) && (orientBL != 0))
3354  break;
3355  }
3356 }
3357 
3358 void orientMeshGFace::operator()(GFace *gf)
3359 {
3360  if (!gf->getNumMeshElements())
3361  return;
3362 
3363  gf->model()->setCurrentMeshEntity(gf);
3364 
3365  if (gf->getMeshMaster() != gf)
3366  {
3367  // It's not clear if periodic meshes should be orientated according to the
3368  // orientation of the underlying CAD surface. Since we don't reorient
3369  // periodic curve meshes, let's also not reorient surface meshes for
3370  // now. This has implications for high-order periodic meshes: see comment in
3371  // FixPeriodicMesh().
3372  }
3373  else if (gf->isFullyDiscrete() || gf->geomType() == GEntity::BoundaryLayerSurface)
3374  {
3375  // Don't do anything
3376  }
3377  else
3378  {
3379  // In old versions we checked the orientation by comparing the orientation
3380  // of a line element on the boundary w.r.t. its connected surface
3381  // element. This is probably better than what follows, but
3382  // * it failed when the 3D Delaunay changes the 1D mesh (since we don't
3383  // recover it yet)
3384  // * it failed with OpenCASCADE geometries, where surface orientions do not
3385  // seem to be consistent with the orientation of the bounding edges
3386 
3387  // Now: orient surface elements w.r.t. normal to geometric model.
3388  // Assumes that originally, orientation is consistent among boundary layer
3389  // (BL) elements, and orientation is consistent among non-BL elements, but
3390  // BL and non-BL elements can be oriented differently
3391 
3392  // Determine whether there is a boundary layer (BL)
3393  BoundaryLayerColumns *blc = gf->getColumns();
3394  const bool existBL = !blc->_toFirst.empty();
3395 
3396  // Get orientation of BL and non-BL elements.
3397  // First, try to get normal to GFace from vertices.
3398  // If it fails, try to get normal to GFace from element barycenter
3399  int orientNonBL = 0, orientBL = existBL ? 0 : 1;
3400  getGFaceOrientation(gf, blc, existBL, true, orientNonBL, orientBL);
3401  if ((orientNonBL == 0) || (orientBL == 0))
3402  getGFaceOrientation(gf, blc, existBL, false, orientNonBL, orientBL);
3403 
3404  // Exit if could not determine orientation of both non-BL el. and BL el.
3405  if ((orientNonBL == 0) && (orientBL == 0))
3406  {
3407  Msg::Warning("Could not orient mesh in surface %d", gf->tag());
3408  return;
3409  }
3410 
3411  // Reverse BL and non-BL elements if needed
3412  if (existBL)
3413  { // If there is a BL, test BL/non-BL elements
3414  if ((orientNonBL == -1) || (orientBL == -1))
3415  {
3416  for (std::size_t iEl = 0; iEl < gf->getNumMeshElements(); iEl++)
3417  {
3418  MElement *e = gf->getMeshElement(iEl);
3419  // If el. outside of BL...
3420  if (blc->_toFirst.find(e) == blc->_toFirst.end())
3421  {
3422  // ... reverse if needed
3423  if (orientNonBL == -1)
3424  e->reverse();
3425  }
3426  else
3427  { // If el. in BL ... reverse if needed
3428  if (orientBL == -1)
3429  e->reverse();
3430  }
3431  }
3432  }
3433  }
3434  else
3435  { // If no BL, reverse all elements if needed
3436  if (orientNonBL == -1)
3437  {
3438  for (std::size_t iEl = 0; iEl < gf->getNumMeshElements(); iEl++)
3439  gf->getMeshElement(iEl)->reverse();
3440  }
3441  }
3442  }
3443 
3444  // Apply user-specified mesh orientation constraints
3445  if (gf->meshAttributes.reverseMesh)
3446  {
3447  for (std::size_t k = 0; k < gf->getNumMeshElements(); k++)
3448  gf->getMeshElement(k)->reverse();
3449  }
3450 }
parametricCoordinates
static void parametricCoordinates(MElement *t, GFace *gf, double u[4], double v[4], MVertex *close=nullptr)
Definition: meshGFaceOptimize.cpp:472
edgeColumn::_c2
const BoundaryLayerData & _c2
Definition: boundaryLayersData.h:44
BDS_Face::deleted
bool deleted
Definition: BDS.h:222
crossprod
SVector3 crossprod(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:150
edgeColumn
Definition: boundaryLayersData.h:43
filterElements.h
MTriangle.h
qualityMeasures.h
GFace::worst_element_shape
double worst_element_shape
Definition: GFace.h:397
BDS_Edge::p1
BDS_Point * p1
Definition: BDS.h:160
meshGFaceBipartiteLabelling.h
BoundaryLayerColumns::getColumns
edgeColumn getColumns(MVertex *v1, MVertex *v2, int side)
Definition: boundaryLayersData.cpp:23
GEntity::BoundaryLayerSurface
@ BoundaryLayerSurface
Definition: GEntity.h:116
GEntity::affineTransform
std::vector< double > affineTransform
Definition: GEntity.h:403
BoundaryLayerColumns::getColumn
const BoundaryLayerData & getColumn(MVertex *v, MEdge e) const
Definition: boundaryLayersData.h:101
meshGeneratorPeriodic
static bool meshGeneratorPeriodic(GFace *gf, int RECUR_ITER, bool repairSelfIntersecting1dMesh, bool debug=true)
Definition: meshGFace.cpp:2266
GFace::edges
virtual std::vector< GEdge * > const & edges() const
Definition: GFace.h:121
GFace::status
GEntity::MeshGenerationStatus status
Definition: GFace.h:395
GFace::firstDer
virtual Pair< SVector3, SVector3 > firstDer(const SPoint2 &param) const =0
GFace::nbTriangle
int nbTriangle
Definition: GFace.h:399
MEdge
Definition: MEdge.h:14
GFace::getMeshingAlgo
int getMeshingAlgo() const
Definition: GFace.cpp:63
meshGFace::operator()
void operator()(GFace *, bool print=true)
Definition: meshGFace.cpp:3152
GEdge::reparamOnFace
virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const
Definition: GEdge.cpp:482
copyMesh
static void copyMesh(GFace *source, GFace *target)
Definition: meshGFace.cpp:351
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
DocRecord::numPoints
int numPoints
Definition: DivideAndConquer.h:93
splitElementsInBoundaryLayerIfNeeded
void splitElementsInBoundaryLayerIfNeeded(GFace *gf)
Definition: meshGFaceOptimize.cpp:1457
BDS_Mesh
Definition: BDS.h:339
Msg::GetVerbosity
static int GetVerbosity()
Definition: GmshMessage.cpp:254
modifyInitialMeshForBoundaryLayers
static void modifyInitialMeshForBoundaryLayers(GFace *gf, std::vector< MQuadrangle * > &blQuads, std::vector< MTriangle * > &blTris, std::set< MVertex * > &verts, bool debug)
Definition: meshGFace.cpp:873
GPoint::y
double y() const
Definition: GPoint.h:22
GFace::recombine
int recombine
Definition: GFace.h:344
GFace.h
improved_translate
static bool improved_translate(GFace *gf, MVertex *vertex, SVector3 &v1, SVector3 &v2)
Definition: meshGFace.cpp:1078
bowyerWatsonParallelograms
void bowyerWatsonParallelograms(GFace *gf, std::map< MVertex *, MVertex * > *equivalence, std::map< MVertex *, SPoint2 > *parametricCoordinates)
Definition: meshGFaceDelaunayInsertion.cpp:1621
MLine::getVertex
virtual MVertex * getVertex(int num)
Definition: MLine.h:45
bowyerWatsonFrontalLayers
void bowyerWatsonFrontalLayers(GFace *gf, bool quad, std::map< MVertex *, MVertex * > *equivalence, std::map< MVertex *, SPoint2 > *parametricCoordinates)
Definition: meshGFaceDelaunayInsertion.cpp:1497
GFace::getEmbeddedVertices
std::vector< GVertex * > getEmbeddedVertices(bool force=false) const
Definition: GFace.cpp:355
GFace
Definition: GFace.h:33
GFace::triangles
std::vector< MTriangle * > triangles
Definition: GFace.h:428
Msg::Info
static void Info(const char *fmt,...)
Definition: GmshMessage.cpp:587
GEdge::lines
std::vector< MLine * > lines
Definition: GEdge.h:42
GEntity::model
GModel * model() const
Definition: GEntity.h:277
angle
double angle(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:157
SPoint2
Definition: SPoint2.h:12
OS.h
MElement::setPartition
virtual void setPartition(int num)
Definition: MElement.h:93
PolyMesh::HalfEdge::v
Vertex * v
Definition: meshPolyMesh.h:38
Msg::Debug
static void Debug(const char *fmt,...)
Definition: GmshMessage.cpp:752
BDS_Mesh::find_point
BDS_Point * find_point(int num)
Definition: BDS.cpp:277
robustPredicates.h
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
MESH_TRANSFINITE
#define MESH_TRANSFINITE
Definition: GmshDefines.h:259
discreteEdge
Definition: discreteEdge.h:12
MVertex
Definition: MVertex.h:24
buildAdditionalPoints2D
bool buildAdditionalPoints2D(GFace *gf)
Definition: boundaryLayersData.cpp:19
MElementOctree.h
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
BackgroundMesh.h
BDS_Edge::g
BDS_GeomEntity * g
Definition: BDS.h:161
GFace::method
char method
Definition: GFace.h:348
edgeColumn::_c1
const BoundaryLayerData & _c1
Definition: boundaryLayersData.h:44
BDS_Mesh::get_geom
BDS_GeomEntity * get_geom(int p1, int p2)
Definition: BDS.cpp:597
quadMeshRemoveHalfOfOneDMesh::_restore
void _restore()
Definition: meshGFace.cpp:264
reparamMeshVertexOnFace
bool reparamMeshVertexOnFace(MVertex const *v, const GFace *gf, SPoint2 &param, bool onSurface, bool failOnSeam, int dir)
Definition: MVertex.cpp:522
isMeshValid
static bool isMeshValid(GFace *gf)
Definition: meshGFace.cpp:3134
PointRecord::where
DPoint where
Definition: DivideAndConquer.h:28
SPoint3
Definition: SPoint3.h:14
addOrRemove
static void addOrRemove(MVertex *v1, MVertex *v2, std::set< MEdge, MEdgeLessThan > &bedges, std::set< MEdge, MEdgeLessThan > &removed)
Definition: meshGFace.cpp:855
BoundaryLayerColumns::endf
iterf endf()
Definition: boundaryLayersData.h:84
contextMeshOptions::recombineOptimizeTopology
int recombineOptimizeTopology
Definition: Context.h:30
MVertex::getNum
std::size_t getNum() const
Definition: MVertex.h:86
quadMeshRemoveHalfOfOneDMesh::_middle
std::map< MEdge, MVertex *, MEdgeLessThan > _middle
Definition: meshGFace.cpp:139
ALGO_2D_PACK_PRLGRMS
#define ALGO_2D_PACK_PRLGRMS
Definition: GmshDefines.h:245
MFace::normal
SVector3 normal() const
Definition: MFace.cpp:204
GFace::meshAttributes
struct GFace::@18 meshAttributes
SBoundingBox3d::min
SPoint3 min() const
Definition: SBoundingBox3d.h:90
Range::high
T high() const
Definition: Range.h:20
GFace::point
virtual GPoint point(double par1, double par2) const =0
SVector3
Definition: SVector3.h:16
deMeshGFace::operator()
void operator()(GFace *)
Definition: meshGFace.cpp:3106
meshGFaceOptimize.h
contextMeshOptions::recombineNodeRepositioning
int recombineNodeRepositioning
Definition: Context.h:31
meshGFaceBamg
void meshGFaceBamg(GFace *gf)
Definition: meshGFaceBamg.cpp:242
remeshUnrecoveredEdges
static void remeshUnrecoveredEdges(std::multimap< MVertex *, BDS_Point * > &recoverMultiMapInv, std::set< EdgeToRecover > &edgesNotRecovered, bool all)
Definition: meshGFace.cpp:534
BDS_Point::g
BDS_GeomEntity * g
Definition: BDS.h:62
GEdgeLoop::begin
iter begin()
Definition: GEdgeLoop.h:48
GEntity::getNativeType
virtual ModelType getNativeType() const
Definition: GEntity.h:268
meshGEdge.h
quadMeshRemoveHalfOfOneDMesh::quadMeshRemoveHalfOfOneDMesh
quadMeshRemoveHalfOfOneDMesh(GFace *gf, bool periodic)
Definition: meshGFace.cpp:281
MTriangle::getVertex
virtual MVertex * getVertex(int num)
Definition: MTriangle.h:62
BDS_Mesh::del_face
void del_face(BDS_Face *t)
Definition: BDS.cpp:515
MVertex::onWhat
GEntity * onWhat() const
Definition: MVertex.h:82
DocRecord::points
PointRecord * points
Definition: DivideAndConquer.h:95
MVertex::point
SPoint3 point() const
Definition: MVertex.h:67
BoundaryLayerColumns
Definition: boundaryLayersData.h:51
BDS_Mesh::triangles
std::vector< BDS_Face * > triangles
Definition: BDS.h:349
DivideAndConquer.h
GmshMessage.h
MLine.h
edges
static int edges[6][2]
Definition: meshGRegionLocalMeshMod.cpp:23
GEntity::GmshModel
@ GmshModel
Definition: GEntity.h:81
MAX_LC
#define MAX_LC
Definition: GEntity.h:19
BoundaryLayerColumns::_elemColumns
std::map< MElement *, std::vector< MElement * > > _elemColumns
Definition: boundaryLayersData.h:57
BDS_GeomEntity::classif_degree
int classif_degree
Definition: BDS.h:34
DocRecord
Definition: DivideAndConquer.h:64
GEntity
Definition: GEntity.h:31
GPoint
Definition: GPoint.h:13
BDS_Point::lcBGM
double & lcBGM()
Definition: BDS.h:65
GFace::closestPoint
virtual GPoint closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const
Definition: GFace.cpp:1323
GEntity::getMeshMaster
GEntity * getMeshMaster() const
Definition: GEntity.h:291
debugSurface
int debugSurface
Definition: meshGFace.cpp:3150
PolyMesh::HalfEdge::next
HalfEdge * next
Definition: meshPolyMesh.h:41
MElement::getNumVertices
virtual std::size_t getNumVertices() const =0
GFace::quadrangles
std::vector< MQuadrangle * > quadrangles
Definition: GFace.h:429
BDS_Edge::oppositeof
void oppositeof(BDS_Point *oface[2]) const
Definition: BDS.cpp:572
BDS_Edge::p2
BDS_Point * p2
Definition: BDS.h:160
Triangle::c
PointNumero c
Definition: DivideAndConquer.h:61
GFace::storage1
std::vector< SPoint3 > storage1
Definition: GFace.h:458
Fopen
FILE * Fopen(const char *f, const char *mode)
Definition: OS.cpp:273
PointRecord::adjacent
DListPeek adjacent
Definition: DivideAndConquer.h:29
GEntity::Plane
@ Plane
Definition: GEntity.h:105
MLine
Definition: MLine.h:21
GPoint::z
double z() const
Definition: GPoint.h:23
GFace::vertices
virtual std::vector< GVertex * > vertices() const
Definition: GFace.cpp:390
MEdgeVertex::getLc
double getLc() const
Definition: MVertex.h:165
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
Pair::first
L first() const
Definition: Pair.h:22
ALGO_2D_FRONTAL_QUAD
#define ALGO_2D_FRONTAL_QUAD
Definition: GmshDefines.h:244
MElement::getVertex
virtual const MVertex * getVertex(int num) const =0
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
GEntity::periodic
virtual bool periodic(int dim) const
Definition: GEntity.h:244
DocRecord::numTriangles
int numTriangles
Definition: DivideAndConquer.h:96
Triangle::a
PointNumero a
Definition: DivideAndConquer.h:61
GEdge.h
GEdgeLoop
Definition: GEdgeLoop.h:36
GFace::best_element_shape
double best_element_shape
Definition: GFace.h:397
BoundaryLayerColumns::beginf
iterf beginf()
Definition: boundaryLayersData.h:83
orientMeshGFace::operator()
void operator()(GFace *)
Definition: meshGFace.cpp:3358
Range
Definition: Range.h:10
algoDelaunay2D
static bool algoDelaunay2D(GFace *gf)
Definition: meshGFace.cpp:762
MEdge::getVertex
MVertex * getVertex(std::size_t i) const
Definition: MEdge.h:39
Pair::second
R second() const
Definition: Pair.h:23
GFace::edgeLoops
std::vector< GEdgeLoop > edgeLoops
Definition: GFace.h:47
GFace::meshStatistics
struct GFace::@19 meshStatistics
EdgeToRecover
Definition: BDS.h:316
GPoint::u
double u() const
Definition: GPoint.h:27
GEntity::mesh_vertices
std::vector< MVertex * > mesh_vertices
Definition: GEntity.h:56
contextMeshOptions::recombineMinimumQuality
double recombineMinimumQuality
Definition: Context.h:32
boundaryLayersData.h
MEdgeVertex
Definition: MVertex.h:137
HighOrder.h
GVertex
Definition: GVertex.h:23
BoundaryLayerData
Definition: boundaryLayersData.h:24
BDS_Mesh::points
std::set< BDS_Point *, PointLessThan > points
Definition: BDS.h:347
GEntity::getTypeString
virtual std::string getTypeString()
Definition: GEntity.h:133
ALGO_2D_FRONTAL
#define ALGO_2D_FRONTAL
Definition: GmshDefines.h:242
BoundaryLayerColumns::_normals
std::multimap< MEdge, SVector3, MEdgeLessThan > _normals
Definition: boundaryLayersData.h:64
BDS_GeomEntity
Definition: BDS.h:31
MVertex.h
recur_tag
void recur_tag(BDS_Face *t, BDS_GeomEntity *g)
Definition: BDS.cpp:605
qmTriangle::gamma
static double gamma(MTriangle *f)
Definition: qualityMeasures.cpp:146
GFace::embeddedVertices
std::set< GVertex *, GEntityPtrLessThan > & embeddedVertices()
Definition: GFace.h:160
norm
void norm(const double *vec, double *norm)
Definition: gmshLevelset.cpp:202
Numeric.h
BDS_Point::v
double v
Definition: BDS.h:57
MElement::getFace
virtual MFace getFace(int num) const =0
GEntity::DONE
@ DONE
Definition: GEntity.h:130
GFace::isFullyDiscrete
virtual bool isFullyDiscrete()
Definition: GFace.cpp:2945
PolyMesh::Vertex
Definition: meshPolyMesh.h:21
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
bowyerWatson
void bowyerWatson(GFace *gf, int MAXPNT, std::map< MVertex *, MVertex * > *equivalence, std::map< MVertex *, SPoint2 > *parametricCoordinates)
Definition: meshGFaceDelaunayInsertion.cpp:1033
GEdge::getBeginVertex
virtual GVertex * getBeginVertex() const
Definition: GEdge.h:63
meshGFaceDelaunayInsertion.h
GFace::parFromPoint
virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true, bool convTestXYZ=false) const
Definition: GFace.cpp:1254
BDS_Mesh::cleanup
void cleanup()
Definition: BDS.cpp:642
MESH_NONE
#define MESH_NONE
Definition: GmshDefines.h:258
MElement::reverse
virtual void reverse()
Definition: MElement.h:308
DPoint::h
double h
Definition: DivideAndConquer.h:24
GEntity::geomType
virtual GeomType geomType() const
Definition: GEntity.h:238
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
GEntity::compound
std::vector< GEntity * > compound
Definition: GEntity.h:59
PolyMesh::faces
std::vector< Face * > faces
Definition: meshPolyMesh.h:61
CTX::mesh
contextMeshOptions mesh
Definition: Context.h:313
GEdge::isSeam
virtual bool isSeam(const GFace *face) const
Definition: GEdge.h:101
GFace::storage2
std::vector< SVector3 > storage2
Definition: GFace.h:459
ALGO_2D_MESHADAPT
#define ALGO_2D_MESHADAPT
Definition: GmshDefines.h:238
deMeshGFace
Definition: meshGFace.h:30
MElement
Definition: MElement.h:30
getGFaceOrientation
static void getGFaceOrientation(GFace *gf, BoundaryLayerColumns *blc, bool existBL, bool fromVert, int &orientNonBL, int &orientBL)
Definition: meshGFace.cpp:3330
BDS.h
GEntity::tag
int tag() const
Definition: GEntity.h:280
MFaceVertex
Definition: MVertex.h:168
meshGenerator
static bool meshGenerator(GFace *, int, bool, int, bool, std::vector< GEdge * > *)
Definition: meshGFace.cpp:1220
discreteFace.h
BDS_Point::iD
int iD
Definition: BDS.h:61
GEdge::parBoundsOnFace
virtual Range< double > parBoundsOnFace(GFace *face=nullptr) const
Definition: GEdge.h:223
BoundaryLayerData::_column
std::vector< MVertex * > _column
Definition: boundaryLayersData.h:26
SVector3::x
double x(void) const
Definition: SVector3.h:30
SBoundingBox3d::makeCube
void makeCube()
Definition: SBoundingBox3d.h:94
element
Definition: shapeFunctions.h:12
PointRecord::data
void * data
Definition: DivideAndConquer.h:30
GFace::deleteMesh
virtual void deleteMesh()
Definition: GFace.cpp:160
BDS_Edge::faces
BDS_Face * faces(std::size_t const i) const
Definition: BDS.h:103
deleteUnusedVertices
static void deleteUnusedVertices(GFace *gf)
Definition: meshGFace.cpp:1190
GEntity::parBounds
virtual Range< double > parBounds(int i) const
Definition: GEntity.h:259
PolyMesh::Vertex::data
int data
Definition: meshPolyMesh.h:29
filterOverlappingElements
static void filterOverlappingElements(std::vector< MLine * > &lines, std::vector< MElement * > &els, std::map< MElement *, std::vector< MElement * > > &_elemColumns, std::map< MElement *, MElement * > &_toFirst)
Definition: filterElements.cpp:215
GFace::average_element_shape
double average_element_shape
Definition: GFace.h:397
refineMeshBDS
void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT, const bool computeNodalSizeField, std::map< MVertex *, BDS_Point * > *recoverMapInv, std::map< BDS_Point *, MVertex *, PointLessThan > *recoverMap, std::vector< SPoint2 > *true_boundary)
Definition: meshGFaceBDS.cpp:726
BoundaryLayerColumns::getNbColumns
int getNbColumns(MVertex *v) const
Definition: boundaryLayersData.h:118
BDS_Edge::numfaces
int numfaces() const
Definition: BDS.h:110
GFace::getMeshElement
MElement * getMeshElement(std::size_t index) const
Definition: GFace.cpp:234
MTriangle
Definition: MTriangle.h:26
GFace::getColumns
BoundaryLayerColumns * getColumns()
Definition: GFace.h:453
GModel::writeMSH
int writeMSH(const std::string &name, double version=2.2, bool binary=false, bool saveAll=false, bool saveParametric=false, double scalingFactor=1.0, int elementStartNum=0, int saveSinglePartition=0, bool append=false)
Definition: GModelIO_MSH.cpp:76
BDS_Face
Definition: BDS.h:164
computeElementShapes
static void computeElementShapes(GFace *gf, double &worst, double &avg, double &best, int &nT, int &greaterThan)
Definition: meshGFace.cpp:114
BDS_Mesh::add_geom
void add_geom(int degree, int tag)
Definition: BDS.cpp:539
GFaceInitialMesh
PolyMesh * GFaceInitialMesh(int faceTag, int recover, std::vector< double > *additional)
Definition: meshTriangulation.cpp:753
SPoint3::distance
double distance(const SPoint3 &p) const
Definition: SPoint3.h:176
GFace::getEmbeddedEdges
std::vector< GEdge * > getEmbeddedEdges(bool force=false) const
Definition: GFace.cpp:362
meshGFaceQuadrangulateBipartiteLabelling
void meshGFaceQuadrangulateBipartiteLabelling(int faceTag)
Definition: meshGFaceBipartiteLabelling.cpp:358
BDS_Face::getNodes
bool getNodes(BDS_Point *_n[4]) const
Definition: BDS.h:201
GFace::setMeshingAlgo
void setMeshingAlgo(int val)
Definition: GFace.h:371
recombineIntoQuads
void recombineIntoQuads(GFace *gf, bool blossom, int topologicalOptiPasses, bool nodeRepositioning, double minqual)
Definition: meshGFaceOptimize.cpp:1311
DocRecord::triangles
Triangle * triangles
Definition: DivideAndConquer.h:97
BoundaryLayerColumns::_toFirst
std::map< MElement *, MElement * > _toFirst
Definition: boundaryLayersData.h:56
ALGO_3D_RTREE
#define ALGO_3D_RTREE
Definition: GmshDefines.h:254
MVertex::getParameter
virtual bool getParameter(int i, double &par) const
Definition: MVertex.h:97
DPoint::v
double v
Definition: DivideAndConquer.h:23
meshGFace.h
outputScalarField
void outputScalarField(std::vector< BDS_Face * > &t, const char *iii, int param, GFace *gf)
Definition: BDS.cpp:46
SPoint2::x
double x(void) const
Definition: SPoint2.h:86
discreteEdge.h
CTX::debugSurface
int debugSurface
Definition: Context.h:138
Context.h
SPoint2::distance
double distance(const SPoint2 &p) const
Definition: SPoint2.h:124
getGFaceNormalFromBary
static bool getGFaceNormalFromBary(GFace *gf, MElement *el, SVector3 &nf)
Definition: meshGFace.cpp:3308
GFace::storage3
std::vector< SVector3 > storage3
Definition: GFace.h:460
backgroundMesh::unset
static void unset()
Definition: BackgroundMesh.cpp:60
meshTriangulation.h
SVector3::y
double y(void) const
Definition: SVector3.h:31
laplaceSmoothing
void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
Definition: meshGFaceOptimize.cpp:978
GVertex.h
buildBackgroundMesh
void buildBackgroundMesh(GFace *gf, bool crossFieldClosestPoint, std::map< MVertex *, MVertex * > *equivalence, std::map< MVertex *, SPoint2 > *parametricCoordinates)
Definition: meshGFaceDelaunayInsertion.cpp:1453
GFace::edgeCounterparts
std::map< GEdge *, std::pair< GEdge *, int > > edgeCounterparts
Definition: GFace.h:50
GFace::embeddedEdges
std::vector< GEdge * > & embeddedEdges()
Definition: GFace.h:156
pointInsideParametricDomain
bool pointInsideParametricDomain(std::vector< SPoint2 > &bnd, SPoint2 &p, SPoint2 &out, int &N)
Definition: meshGFace.cpp:43
MQuadrangle.h
contextMeshOptions::lcFactor
double lcFactor
Definition: Context.h:21
BDS_Mesh::recover_edge
BDS_Edge * recover_edge(int p1, int p2, bool &_fatal, std::set< EdgeToRecover > *e2r=nullptr, std::set< EdgeToRecover > *not_recovered=nullptr)
Definition: BDS.cpp:347
meshGFaceBDS.h
contextMeshOptions::algo2d
int algo2d
Definition: Context.h:29
BDS_Mesh::add_point
BDS_Point * add_point(int num, double x, double y, double z)
Definition: BDS.cpp:257
ALGO_2D_INITIAL_ONLY
#define ALGO_2D_INITIAL_ONLY
Definition: GmshDefines.h:240
quadMeshRemoveHalfOfOneDMesh::finish
void finish()
Definition: meshGFace.cpp:333
quadMeshRemoveHalfOfOneDMesh::_gf
GFace * _gf
Definition: meshGFace.cpp:137
BDS_Mesh::del_edge
void del_edge(BDS_Edge *e)
Definition: BDS.cpp:525
GEdge::point
virtual GPoint point(double p) const =0
Pair
Definition: Pair.h:10
BDS_Mesh::del_point
void del_point(BDS_Point *p)
Definition: BDS.cpp:533
GEdge
Definition: GEdge.h:26
quadMeshRemoveHalfOfOneDMesh::_subdivide
void _subdivide()
Definition: meshGFace.cpp:140
GFace::nbGoodQuality
int nbGoodQuality
Definition: GFace.h:400
modifyInitialMeshToRemoveDegeneracies
void modifyInitialMeshToRemoveDegeneracies(GFace *gf, BDS_Mesh &m, std::map< BDS_Point *, MVertex *, PointLessThan > *recoverMap)
Definition: meshGFaceBDS.cpp:690
directions_storage
static void directions_storage(GFace *gf)
Definition: meshGFace.cpp:1111
GFace::refineAllEdges
bool refineAllEdges
Definition: GFace.h:396
backgroundMesh::current
static backgroundMesh * current()
Definition: BackgroundMesh.cpp:68
ALGO_2D_AUTO
#define ALGO_2D_AUTO
Definition: GmshDefines.h:239
Msg::GetNumThreads
static int GetNumThreads()
Definition: GmshMessage.cpp:1635
getGFaceNormalFromVert
static bool getGFaceNormalFromVert(GFace *gf, MElement *el, SVector3 &nf)
Definition: meshGFace.cpp:3291
GEntity::vertexCounterparts
std::map< GVertex *, GVertex * > vertexCounterparts
Definition: GEntity.h:62
trueBoundary
static void trueBoundary(GFace *gf, std::vector< SPoint2 > &bnd, int debug)
Definition: meshGFace.cpp:66
GPoint::v
double v() const
Definition: GPoint.h:28
BDS_Edge
Definition: BDS.h:85
BDS_Point::u
double u
Definition: BDS.h:57
GPoint.h
GModel.h
BDS_Point
Definition: BDS.h:49
BDS_Mesh::edges
std::vector< BDS_Edge * > edges
Definition: BDS.h:348
MVertex::y
double y() const
Definition: MVertex.h:61
GFace::normal
virtual SVector3 normal(const SPoint2 &param) const
Definition: GFace.cpp:1416
contextMeshOptions::algoRecombine
int algoRecombine
Definition: Context.h:30
GEdge::getEndVertex
virtual GVertex * getEndVertex() const
Definition: GEdge.h:64
BDS_Point::lc
double & lc()
Definition: BDS.h:66
backgroundMesh::getAngle
double getAngle(double u, double v, double w) const
Definition: BackgroundMesh.cpp:603
Range::low
T low() const
Definition: Range.h:18
MeshTransfiniteSurface
int MeshTransfiniteSurface(GFace *gf)
Definition: meshGFaceTransfinite.cpp:166
GFace::getNumMeshElements
std::size_t getNumMeshElements() const
Definition: GFace.cpp:181
Triangle::b
PointNumero b
Definition: DivideAndConquer.h:61
GFace::nbEdge
int nbEdge
Definition: GFace.h:399
GEntity::FAILED
@ FAILED
Definition: GEntity.h:130
SBoundingBox3d::max
SPoint3 max() const
Definition: SBoundingBox3d.h:91
ALGO_2D_DELAUNAY
#define ALGO_2D_DELAUNAY
Definition: GmshDefines.h:241
PolyMesh
Definition: meshPolyMesh.h:15
GEntity::dim
virtual int dim() const
Definition: GEntity.h:196
meshGFace::repairSelfIntersecting1dMesh
const bool repairSelfIntersecting1dMesh
Definition: meshGFace.h:22
GEdgeLoop::end
iter end()
Definition: GEdgeLoop.h:49
delaunayizeBDS
void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
Definition: meshGFaceBDS.cpp:236
GPoint::x
double x() const
Definition: GPoint.h:21
GEntity::getVisibility
virtual char getVisibility()
Definition: GEntity.cpp:35
GEntity::correspondingVertices
std::map< MVertex *, MVertex * > correspondingVertices
Definition: GEntity.h:406
BDS_Face::g
BDS_GeomEntity * g
Definition: BDS.h:224
GEntity::PENDING
@ PENDING
Definition: GEntity.h:130
DocRecord::MakeMeshWithPoints
void MakeMeshWithPoints()
Definition: DivideAndConquer.cpp:855
c1
const double c1
Definition: GaussQuadratureHex.cpp:16
MeshExtrudedSurface
int MeshExtrudedSurface(GFace *gf, std::set< std::pair< MVertex *, MVertex * > > *constrainedEdges=nullptr)
Definition: meshGFaceExtruded.cpp:291
PolyMesh::Vertex::position
SVector3 position
Definition: meshPolyMesh.h:27
MQuadrangle
Definition: MQuadrangle.h:26
robustPredicates::orient2d
REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
Definition: robustPredicates.cpp:1633
ALGO_2D_QUAD_QUASI_STRUCT
#define ALGO_2D_QUAD_QUASI_STRUCT
Definition: GmshDefines.h:247
SBoundingBox3d
Definition: SBoundingBox3d.h:21
ALGO_2D_BAMG
#define ALGO_2D_BAMG
Definition: GmshDefines.h:243
quadMeshRemoveHalfOfOneDMesh
Definition: meshGFace.cpp:135
GFace::storage4
std::vector< double > storage4
Definition: GFace.h:461
GModel::setCurrentMeshEntity
void setCurrentMeshEntity(GEntity *e)
Definition: GModel.cpp:2202
PolyMesh::vertices
std::vector< Vertex * > vertices
Definition: meshPolyMesh.h:59
TRIANGLE_VALIDITY
static double TRIANGLE_VALIDITY(GFace *gf, MTriangle *t)
Definition: meshGFace.cpp:3115
meshGFaceBamg.h
BDS_Mesh::MAXPOINTNUMBER
int MAXPOINTNUMBER
Definition: BDS.h:341
GFace::reverseMesh
bool reverseMesh
Definition: GFace.h:359
BDS_Mesh::find_edge
BDS_Edge * find_edge(int p1, int p2)
Definition: BDS.cpp:301
BDS_Mesh::add_triangle
BDS_Face * add_triangle(int p1, int p2, int p3)
Definition: BDS.cpp:496
GFace::set
void set(const std::vector< GEdge * > &f)
Definition: GFace.h:125
getGEdge
static GEdge * getGEdge(GFace *gf, MVertex *v1, MVertex *v2)
Definition: meshGFace.cpp:2248
recoverEdge
static bool recoverEdge(BDS_Mesh *m, GFace *gf, GEdge *ge, std::map< MVertex *, BDS_Point * > &recoverMapInv, std::set< EdgeToRecover > *e2r, std::set< EdgeToRecover > *notRecovered, int pass)
Definition: meshGFace.cpp:779
MVertex::x
double x() const
Definition: MVertex.h:60
PolyMesh::HalfEdge
Definition: meshPolyMesh.h:32
bowyerWatsonFrontal
void bowyerWatsonFrontal(GFace *gf, std::map< MVertex *, MVertex * > *equivalence, std::map< MVertex *, SPoint2 > *parametricCoordinates, std::vector< SPoint2 > *true_boundary)
Definition: meshGFaceDelaunayInsertion.cpp:1289
SVector3::normalize
double normalize()
Definition: SVector3.h:38
BDS2GMSH
static void BDS2GMSH(BDS_Mesh *m, GFace *gf, std::map< BDS_Point *, MVertex *, PointLessThan > &recoverMap)
Definition: meshGFace.cpp:1150
buildConsecutiveListOfVertices
static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, std::vector< BDS_Point * > &result, SBoundingBox3d &bbox, BDS_Mesh *m, std::map< BDS_Point *, MVertex *, PointLessThan > &recoverMap, int &count, int countTot, double tol, bool seam_the_first=false)
Definition: meshGFace.cpp:1940
ALGO_2D_PACK_PRLGRMS_CSTR
#define ALGO_2D_PACK_PRLGRMS_CSTR
Definition: GmshDefines.h:246
quadMeshRemoveHalfOfOneDMesh::_backup
std::map< GEdge *, std::vector< MLine * > > _backup
Definition: meshGFace.cpp:138
SPoint2::y
double y(void) const
Definition: SPoint2.h:88
backgroundMesh::set
static void set(GFace *)
Definition: BackgroundMesh.cpp:40
contextMeshOptions::randFactor
double randFactor
Definition: Context.h:21