1 // Gmsh - Copyright (C) 1997-2021 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 <utility>
7 #include <list>
8 #include <map>
9 #include <unordered_map>
10 #include "Context.h"
11 #include "gmsh.h"
12 #include "GModel.h"
13 #include "GFace.h"
14 #include "GEdge.h"
15 #include "MLine.h"
16 #include "MVertex.h"
17 #include "MTriangle.h"
18 #include "MQuadrangle.h"
19 #include "meshTriangulation.h"
20 #include "SBoundingBox3d.h"
21 #include "robustPredicates.h"
22 #include "meshGFaceDelaunayInsertion.h"
23 #include "qualityMeasures.h"
24 #include "Numeric.h"
25 #include "SPoint3.h"
26
swap(double & a,double & b)27 void swap(double &a, double &b)
28 {
29 double temp = a;
30 a = b;
31 b = temp;
32 }
33
HilbertCoordinates(double x,double y,double x0,double y0,double xRed,double yRed,double xBlue,double yBlue)34 size_t HilbertCoordinates(double x, double y, double x0, double y0, double xRed,
35 double yRed, double xBlue, double yBlue)
36 {
37 size_t BIG = 1073741824;
38 size_t RESULT = 0;
39 for(int i = 0; i < 16; i++) {
40 double coordRed = (x - x0) * xRed + (y - y0) * yRed;
41 double coordBlue = (x - x0) * xBlue + (y - y0) * yBlue;
42 xRed /= 2;
43 yRed /= 2;
44 xBlue /= 2;
45 yBlue /= 2;
46 if(coordRed <= 0 && coordBlue <= 0) { // quadrant 0
47 x0 -= (xBlue + xRed);
48 y0 -= (yBlue + yRed);
49 swap(xRed, xBlue);
50 swap(yRed, yBlue);
51 }
52 else if(coordRed <= 0 && coordBlue >= 0) { // quadrant 1
53 RESULT += BIG;
54 x0 += (xBlue - xRed);
55 y0 += (yBlue - yRed);
56 }
57 else if(coordRed >= 0 && coordBlue >= 0) { // quadrant 2
58 RESULT += 2 * BIG;
59 x0 += (xBlue + xRed);
60 y0 += (yBlue + yRed);
61 }
62 else if(coordRed >= 0 && coordBlue <= 0) { // quadrant 3
63 x0 += (-xBlue + xRed);
64 y0 += (-yBlue + yRed);
65 swap(xRed, xBlue);
66 swap(yRed, yBlue);
67 xBlue = -xBlue;
68 yBlue = -yBlue;
69 xRed = -xRed;
70 yRed = -yRed;
71 RESULT += 3 * BIG;
72 }
73 else
74 Msg::Warning("Hilbert failed %d %d", coordRed, coordBlue);
75 BIG /= 4;
76 }
77 return RESULT;
78 }
79
80 struct pair_hash {
81 template <class T1, class T2>
operator ()pair_hash82 std::size_t operator()(const std::pair<T1, T2> &pair) const
83 {
84 return std::hash<T1>()(pair.first) ^ std::hash<T2>()(pair.second);
85 }
86 };
87
PolyMesh2GFace(PolyMesh * pm,int faceTag)88 int PolyMesh2GFace(PolyMesh *pm, int faceTag)
89 {
90 GFace *gf = GModel::current()->getFaceByTag(faceTag);
91
92 if(!gf){
93 Msg::Error("PolyMesh2GFace cannot find surface %d", faceTag);
94 return 0;
95 }
96
97 for(auto t : gf->triangles) delete t;
98 for(auto q : gf->quadrangles) delete q;
99 gf->triangles.clear();
100 gf->quadrangles.clear();
101
102 std::unordered_map<int, MVertex *> news;
103 std::unordered_map<PolyMesh::HalfEdge *, GPoint> hon;
104
105 if(!pm->high_order_nodes.empty()) {
106 for(size_t i = 0; i < pm->high_order_nodes.size(); i++) {
107 auto it = hon.find(pm->hedges[i]);
108 if(it == hon.end()) {
109 double uv[2] = {0, 0};
110 SVector3 p = pm->high_order_nodes[i];
111 GPoint gp = gf->closestPoint(SPoint3(p.x(), p.y(), p.z()), uv);
112 if(!gp.succeeded()) {
113 gp.x() = p.x();
114 gp.y() = p.y();
115 gp.z() = p.z();
116 }
117 hon[pm->hedges[i]] = gp;
118 if(pm->hedges[i]->opposite) hon[pm->hedges[i]] = gp;
119 }
120 }
121 }
122
123 std::map<MEdge, GPoint, MEdgeLessThan> hop;
124
125 for(auto f : pm->faces) {
126 if(f->data == faceTag) {
127 PolyMesh::Vertex *v[4] = {f->he->v, f->he->next->v, f->he->next->next->v,
128 f->he->next->next->next->v};
129 MVertex *v_gmsh[4];
130
131 for(int i = 0; i < pm->num_sides(f->he); i++) {
132 if(v[i]->data != -1) {
133 auto it = news.find(v[i]->data);
134 if(it == news.end())
135 v_gmsh[i] = GModel::current()->getMeshVertexByTag(v[i]->data);
136 else
137 v_gmsh[i] = it->second;
138 }
139 else {
140 double uv[2] = {0, 0};
141 GPoint gp = gf->closestPoint(
142 SPoint3(v[i]->position.x(), v[i]->position.y(), v[i]->position.z()),
143 uv);
144 if(gp.succeeded())
145 v_gmsh[i] =
146 new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
147 else
148 v_gmsh[i] = new MFaceVertex(v[i]->position.x(), v[i]->position.y(),
149 v[i]->position.z(), gf,
150 v[i]->position.x(), v[i]->position.y());
151 gf->mesh_vertices.push_back(v_gmsh[i]);
152 v[i]->data = v_gmsh[i]->getNum();
153 news[v[i]->data] = v_gmsh[i];
154 }
155 }
156
157 if(hon.size()) {
158 MEdge l01(v_gmsh[0], v_gmsh[1]);
159 hop[l01] = hon[f->he];
160 MEdge l12(v_gmsh[1], v_gmsh[2]);
161 hop[l12] = hon[f->he->next];
162 MEdge l20(v_gmsh[2], v_gmsh[0]);
163 hop[l20] = hon[f->he->next->next];
164 }
165
166 if(pm->num_sides(f->he) == 3)
167 gf->triangles.push_back(new MTriangle(v_gmsh[0], v_gmsh[1], v_gmsh[2]));
168 else if(pm->num_sides(f->he) == 4)
169 gf->quadrangles.push_back(
170 new MQuadrangle(v_gmsh[0], v_gmsh[1], v_gmsh[2], v_gmsh[3]));
171 }
172 }
173
174 if(!hon.empty()) {
175 GModel::current()->setOrderN(2, 0, 0);
176 #if 1
177 for(auto t : gf->triangles) {
178 for(int i = 0; i < 3; i++) {
179 MEdge li(t->getVertex(i), t->getVertex((i + 1) % 3));
180 GPoint gp = hop[li];
181 MVertex *vint = t->getVertex(i + 3);
182 vint->x() = gp.x();
183 vint->y() = gp.y();
184 vint->z() = gp.z();
185 }
186 }
187 #endif
188 }
189
190 CTX::instance()->mesh.changed = ENT_ALL;
191
192 return 0;
193 }
194
GFace2PolyMesh(int faceTag,PolyMesh ** pm)195 int GFace2PolyMesh(int faceTag, PolyMesh **pm)
196 {
197 // FIXME should probably not use the public API here (and certainly not
198 // initialize it!)
199 gmsh::initialize();
200 *pm = new PolyMesh;
201
202 std::unordered_map<size_t, size_t> nodeLabels;
203 std::unordered_map<std::pair<size_t, size_t>, PolyMesh::HalfEdge *, pair_hash>
204 opposites;
205
206 // FIXME should probably not use the public API here
207 std::vector<int> elementTypes;
208 std::vector<std::vector<std::size_t> > elementTags;
209 std::vector<std::vector<std::size_t> > nodeTags;
210 gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTags, 2,
211 faceTag);
212
213 for(size_t K = 0; K < elementTypes.size(); K++) {
214 int eT = elementTypes[K];
215 int nNod = 0;
216 if(eT == 2)
217 nNod = 3;
218 else if(eT == 3)
219 nNod = 4;
220 else {
221 Msg::Error("GFace2PolyMesh only support quads (element type 3) and "
222 "triangles (element type 2)");
223 return -1;
224 }
225 PolyMesh::Vertex *v[4] = {nullptr, nullptr, nullptr, nullptr};
226
227 for(size_t i = 0; i < elementTags[K].size(); i++) {
228 for(int j = 0; j < nNod; j++) {
229 size_t nodeTag = nodeTags[K][nNod * i + j];
230 auto it = nodeLabels.find(nodeTag);
231 if(it == nodeLabels.end()) {
232 // FIXME should probably not use the public API here
233 std::vector<double> coord(3), parametricCoord(3);
234 int entityDim, entityTag;
235 gmsh::model::mesh::getNode(nodeTag, coord, parametricCoord, entityDim,
236 entityTag);
237 v[j] = new PolyMesh::Vertex(coord[0], coord[1], coord[2], nodeTag);
238 nodeLabels[nodeTag] = (*pm)->vertices.size();
239 (*pm)->vertices.push_back(v[j]);
240 }
241 else
242 v[j] = (*pm)->vertices[it->second];
243 }
244
245 PolyMesh::HalfEdge *he[4];
246 for(int j = 0; j < nNod; j++) {
247 he[j] = new PolyMesh::HalfEdge(v[j]);
248 (*pm)->hedges.push_back(he[j]);
249 v[j]->he = he[j];
250 }
251
252 PolyMesh::Face *ff = new PolyMesh::Face(he[0]);
253 (*pm)->faces.push_back(ff);
254
255 for(int j = 0; j < nNod; j++) {
256 he[j]->next = he[(j + 1) % nNod];
257 he[j]->prev = he[(j - 1 + nNod) % nNod];
258 he[j]->f = ff;
259 // size_t n0 = v[j]->data;
260 // size_t n1 = v[(j+1)%nNod]->data;
261 // std::pair<size_t, size_t> pj =
262 // std::make_pair(std::min(n0,n1),std::max(n0,n1));
263 // auto itj = opposites.find(pj);
264 // if(itj == opposites.end()) opposites[pj] = he[j];
265 // else {
266 // he[j]->opposite = itj->second;
267 // itj->second->opposite = he[j];
268 // }
269 }
270 }
271 }
272
273 HalfEdgePtrLessThan compare;
274 std::sort((*pm)->hedges.begin(), (*pm)->hedges.end(), compare);
275
276 HalfEdgePtrEqual equal;
277 for(size_t i = 0; i < (*pm)->hedges.size() - 1; i++) {
278 PolyMesh::HalfEdge *h0 = (*pm)->hedges[i];
279 PolyMesh::HalfEdge *h1 = (*pm)->hedges[i + 1];
280 if(equal(h0, h1)) {
281 h0->opposite = h1;
282 h1->opposite = h0;
283 i++;
284 }
285 }
286 return 0;
287 }
288
delaunayEdgeCriterionPlaneIsotropic(PolyMesh::HalfEdge * he,void *)289 static int delaunayEdgeCriterionPlaneIsotropic(PolyMesh::HalfEdge *he, void *)
290 {
291 if(he->opposite == nullptr) return -1;
292 PolyMesh::Vertex *v0 = he->v;
293 PolyMesh::Vertex *v1 = he->next->v;
294 PolyMesh::Vertex *v2 = he->next->next->v;
295 PolyMesh::Vertex *v = he->opposite->next->next->v;
296
297 // FIXME : should be oriented anyway !
298 double result = -robustPredicates::incircle(v0->position, v1->position,
299 v2->position, v->position);
300
301 return (result > 0) ? 1 : 0;
302 }
303
faceCircumCenter(PolyMesh::HalfEdge * he,GFace * gf,double * res,double * uv)304 static void faceCircumCenter(PolyMesh::HalfEdge *he, GFace *gf, double *res,
305 double *uv)
306 {
307 PolyMesh::Vertex *v0 = he->v;
308 PolyMesh::Vertex *v1 = he->next->v;
309 PolyMesh::Vertex *v2 = he->next->next->v;
310 GPoint p0 = gf->point(v0->position.x(), v0->position.y());
311 GPoint p1 = gf->point(v1->position.x(), v1->position.y());
312 GPoint p2 = gf->point(v2->position.x(), v2->position.y());
313 double q0[3] = {p0.x(), p0.y(), p0.z()};
314 double q1[3] = {p1.x(), p1.y(), p1.z()};
315 double q2[3] = {p2.x(), p2.y(), p2.z()};
316 circumCenterXYZ(q0, q1, q2, res, uv);
317 }
318
faceQuality(PolyMesh::HalfEdge * he,GFace * gf)319 static double faceQuality(PolyMesh::HalfEdge *he, GFace *gf)
320 {
321 PolyMesh::Vertex *v0 = he->v;
322 PolyMesh::Vertex *v1 = he->next->v;
323 PolyMesh::Vertex *v2 = he->next->next->v;
324 GPoint p0 = gf->point(v0->position.x(), v0->position.y());
325 GPoint p1 = gf->point(v1->position.x(), v1->position.y());
326 GPoint p2 = gf->point(v2->position.x(), v2->position.y());
327 return qmTriangle::gamma(p0.x(), p0.y(), p0.z(), p1.x(), p1.y(), p1.z(),
328 p2.x(), p2.y(), p2.z());
329 }
330
331 /*
332 static int qualityCriterion3D(PolyMesh::HalfEdge *he, void *p){
333 if (he->data > 0) return -1;
334 if (he->opposite == nullptr) return -1;
335 if (p == nullptr) return -1;
336
337 GFace *gf = (GFace*)p;
338
339 PolyMesh::Vertex *v0 = he->v;
340 PolyMesh::Vertex *v1 = he->next->v;
341 PolyMesh::Vertex *v2 = he->next->next->v;
342 PolyMesh::Vertex *v3 = he->opposite->next->next->v;
343
344 GPoint p0 = gf->point (v0->position.x(),v0->position.y());
345 GPoint p1 = gf->point (v1->position.x(),v1->position.y());
346 GPoint p2 = gf->point (v2->position.x(),v2->position.y());
347 GPoint p3 = gf->point (v3->position.x(),v3->position.y());
348
349 double q1 = qmTriangle::gamma
350 (p0.x(),p0.y(),p0.z(),p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z()); double q2 =
351 qmTriangle::gamma
352 (p2.x(),p2.y(),p2.z(),p3.x(),p3.y(),p3.z(),p0.x(),p0.y(),p0.z());
353
354 double o1 = qmTriangle::gamma
355 (p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z(),p3.x(),p3.y(),p3.z()); double o2 =
356 qmTriangle::gamma
357 (p3.x(),p3.y(),p3.z(),p0.x(),p0.y(),p0.z(),p1.x(),p1.y(),p1.z());
358
359 return std::max(fabs(q1),fabs(q2)) > std::max(fabs(o1),fabs(o2)) ? 0 : 1;
360 }
361 */
362
Walk(PolyMesh::Face * f,double x,double y)363 static PolyMesh::Face *Walk(PolyMesh::Face *f, double x, double y)
364 {
365 double POS[2] = {x, y};
366 PolyMesh::HalfEdge *he = f->he;
367
368 while(1) {
369 PolyMesh::Vertex *v0 = he->v;
370 PolyMesh::Vertex *v1 = he->next->v;
371 PolyMesh::Vertex *v2 = he->next->next->v;
372
373 double s0 = -robustPredicates::orient2d(v0->position, v1->position, POS);
374 double s1 = -robustPredicates::orient2d(v1->position, v2->position, POS);
375 double s2 = -robustPredicates::orient2d(v2->position, v0->position, POS);
376
377 if(s0 >= 0 && s1 >= 0 && s2 >= 0) {
378 /* printf("Face %g %g %g / %g %g %g / %g %g %g \n",
379 v0->position.x(), v0->position.y(), v0->position.z(),
380 v1->position.x(), v1->position.y(), v1->position.z(),
381 v2->position.x(), v2->position.y(), v2->position.z());
382 printf("point %g %g CURRENT FACE %p %g %g %g\n", x,y,he->f,
383 s0,s1,s2); */
384 // getchar();
385 return he->f;
386 }
387 else if(s0 <= 0 && s1 >= 0 && s2 >= 0)
388 he = he->opposite;
389 else if(s1 <= 0 && s0 >= 0 && s2 >= 0)
390 he = he->next->opposite;
391 else if(s2 <= 0 && s0 >= 0 && s1 >= 0)
392 he = he->next->next->opposite;
393 else if(s0 <= 0 && s1 <= 0)
394 he = s0 > s1 ? he->opposite : he->next->opposite;
395 else if(s0 <= 0 && s2 <= 0)
396 he = s0 > s2 ? he->opposite : he->next->next->opposite;
397 else if(s1 <= 0 && s2 <= 0)
398 he = s1 > s2 ? he->next->opposite : he->next->next->opposite;
399 else {
400 Msg::Error("Could not find half-edge in walk for point %g %g on "
401 "face %g %g %g / %g %g %g / %g %g %g "
402 "(orientation tests %g %g %g)", x, y,
403 v0->position.x(), v0->position.y(), v0->position.z(),
404 v1->position.x(), v1->position.y(), v1->position.z(),
405 v2->position.x(), v2->position.y(), v2->position.z(),
406 s0, s1, s2);
407 }
408 if(he == nullptr) break;
409 }
410 // should only come here wether the triangulated domain is not convex
411 return nullptr;
412 }
413
414 // recover an edge that goes from v_start --> v_end
415 // ----------------------------------- assume it's internal !!!
416
intersect(PolyMesh::Vertex * v0,PolyMesh::Vertex * v1,PolyMesh::Vertex * b0,PolyMesh::Vertex * b1)417 static int intersect(PolyMesh::Vertex *v0, PolyMesh::Vertex *v1,
418 PolyMesh::Vertex *b0, PolyMesh::Vertex *b1)
419 {
420 double s0 =
421 robustPredicates::orient2d(v0->position, v1->position, b0->position);
422 double s1 =
423 robustPredicates::orient2d(v0->position, v1->position, b1->position);
424 if(s0 * s1 >= 0) return 0;
425 double t0 =
426 robustPredicates::orient2d(b0->position, b1->position, v0->position);
427 double t1 =
428 robustPredicates::orient2d(b0->position, b1->position, v1->position);
429 if(t0 * t1 >= 0) return 0;
430 return 1;
431 }
432
recover_edge(PolyMesh * pm,PolyMesh::Vertex * v_start,PolyMesh::Vertex * v_end)433 static int recover_edge(PolyMesh *pm, PolyMesh::Vertex *v_start,
434 PolyMesh::Vertex *v_end)
435 {
436 PolyMesh::HalfEdge *he = v_start->he;
437 std::list<PolyMesh::HalfEdge *> _list;
438
439 do {
440 PolyMesh::Vertex *v1 = he->next->v;
441 if(v1 == v_end) {
442 return 0; // edge exists
443 }
444 PolyMesh::Vertex *v2 = he->next->next->v;
445 if(v2 == v_end) {
446 return 0; // edge exists
447 }
448
449 if(intersect(v_start, v_end, v1, v2)) {
450 // printf("INTERSECTION WITH %d %d\n", v1->data, v2->data);
451 _list.push_back(he->next);
452 break;
453 }
454 he = he->next->next->opposite;
455 } while(he != v_start->he);
456
457 if(_list.empty()) { return -1; }
458
459 // find all intersections
460 while(1) {
461 he = _list.back();
462 he = he->opposite;
463 if(!he) return -2;
464 he = he->next;
465 PolyMesh::Vertex *v1 = he->v;
466 PolyMesh::Vertex *v2 = he->next->v;
467 if(v2 == v_end) {
468 // printf("END FOUND %d\n", v2->data);
469 break;
470 }
471 if(intersect(v_start, v_end, v1, v2)) {
472 // printf("INTESECTION %d %d\n", v1->data, v2->data);
473 _list.push_back(he);
474 }
475 else {
476 he = he->next;
477 v1 = he->v;
478 v2 = he->next->v;
479 if(v2 == v_end) {
480 // printf("END FOUND %d\n", v2->data);
481 break;
482 }
483 if(intersect(v_start, v_end, v1, v2)) {
484 // printf("INTESECTION %d %d\n", v1->data, v2->data);
485 _list.push_back(he);
486 }
487 else {
488 return -3;
489 }
490 }
491 }
492
493 int nbIntersection = _list.size();
494 // printf("%d intersections\n", nbIntersection);
495 // int K = 100;
496 while(!_list.empty()) {
497 he = *_list.begin();
498 _list.erase(_list.begin());
499 // ensure that swap is allowed (convex quad)
500 if(intersect(he->v, he->next->v, he->next->next->v,
501 he->opposite->next->next->v)) {
502 // ensure that swap removes one intersection
503 int still_intersect = intersect(v_start, v_end, he->next->next->v,
504 he->opposite->next->next->v);
505 // printf("swapping %d %d\n", he->v->data, he->next->v->data);
506 pm->swap_edge(he);
507 // pm->print4debug(K++);
508 if(still_intersect) _list.push_back(he);
509 }
510 else
511 _list.push_back(he);
512 }
513 return nbIntersection;
514 }
515
Color(PolyMesh::HalfEdge * he,int color)516 static PolyMesh::HalfEdge *Color(PolyMesh::HalfEdge *he, int color)
517 {
518 std::stack<PolyMesh::Face *> _stack;
519 _stack.push(he->f);
520
521 PolyMesh::HalfEdge *other_side = nullptr;
522
523 while(!_stack.empty()) {
524 PolyMesh::Face *f = _stack.top();
525 _stack.pop();
526 f->data = color;
527 he = f->he;
528 for(int i = 0; i < 3; i++) {
529 if(he->data == -1 && he->opposite != nullptr &&
530 he->opposite->f->data == -1) {
531 _stack.push(he->opposite->f);
532 }
533 else if(he->data != -1 && he->opposite != nullptr) {
534 other_side = he->opposite;
535 }
536 he = he->next;
537 }
538 }
539 return other_side;
540 }
541
GFaceDelaunayRefinement(size_t faceTag)542 void GFaceDelaunayRefinement(size_t faceTag)
543 {
544 GFace *gf = GModel::current()->getFaceByTag(faceTag);
545
546 PolyMesh *pm = GFaceInitialMesh(faceTag, 1);
547
548 std::list<PolyMesh::HalfEdge *> _list;
549 double _limit = 0.7;
550 for(auto f : pm->faces) {
551 double q = faceQuality(f->he, gf);
552 if(q < _limit && f->data == gf->tag()) _list.push_back(f->he);
553 }
554 // int I = 1;
555 while(!_list.empty()) {
556 PolyMesh::HalfEdge *he = *_list.begin();
557 _list.erase(_list.begin());
558 double q = faceQuality(he, gf);
559 if(q < _limit) {
560 double uv[2];
561 SPoint3 cc;
562 faceCircumCenter(he, gf, cc, uv);
563 GPoint gp = gf->closestPoint(cc, uv);
564 if(gp.succeeded()) {
565 PolyMesh::Face *f = he->f;
566 f = Walk(f, gp.u(), gp.v());
567 if(f && f->data == (int)faceTag) {
568 std::vector<PolyMesh::HalfEdge *> _touched;
569 pm->split_triangle(-1, gp.u(), gp.v(), 0, f,
570 delaunayEdgeCriterionPlaneIsotropic, gf,
571 &_touched);
572 if(_touched.size() == 3) {
573 // we should unsplit ...
574 // pm->undo_split(_touched);
575 }
576 else {
577 std::vector<PolyMesh::Face *> _f;
578 for(auto h : _touched)
579 if(std::find(_f.begin(), _f.end(), h->f) == _f.end())
580 _f.push_back(h->f);
581
582 // printf("step %d %lu touched : ", I, _f.size());
583 for(auto pf : _f) {
584 q = faceQuality(pf->he, gf);
585 // printf("%12.5E ", q);
586 if(q < _limit && pf->data == gf->tag()) _list.push_back(pf->he);
587 }
588 // printf("\n");
589 }
590 // pm->print4debug(100000 + I++);
591 }
592 }
593 }
594 }
595 }
596
GFaceDelaunayRefinementOldMesher(int faceTag)597 void GFaceDelaunayRefinementOldMesher(int faceTag)
598 {
599 PolyMesh *pm = GFaceInitialMesh(faceTag);
600
601 GFace *gf = GModel::current()->getFaceByTag(faceTag);
602
603 // use old code ---
604
605 for(auto f : pm->faces) {
606 if(f->data == faceTag) {
607 size_t n0 = f->he->v->data;
608 size_t n1 = f->he->next->v->data;
609 size_t n2 = f->he->next->next->v->data;
610 MVertex *v0 = GModel::current()->getMeshVertexByTag(n0);
611 MVertex *v1 = GModel::current()->getMeshVertexByTag(n1);
612 MVertex *v2 = GModel::current()->getMeshVertexByTag(n2);
613 gf->triangles.push_back(new MTriangle(v0, v1, v2));
614 }
615 }
616 delete pm;
617 // bowyerWatsonFrontal(gf);
618 }
619
620 struct nodeCopies {
621 MVertex *mv;
622 size_t nbCopies;
623 double u[8], v[8]; // max 8 copies -- reduced to 4
624 size_t id[8];
nodeCopiesnodeCopies625 nodeCopies(MVertex *_mv, double _u, double _v) : mv(_mv), nbCopies(1)
626 {
627 u[0] = _u;
628 v[0] = _v;
629 }
addCopynodeCopies630 void addCopy(double _u, double _v)
631 {
632 for(size_t i = 0; i < nbCopies; i++) {
633 if(fabs(u[i] - _u) < 1.e-9 && fabs(v[i] - _v) < 1.e-9) return;
634 }
635 u[nbCopies] = _u;
636 v[nbCopies] = _v;
637 nbCopies++;
638 }
closestnodeCopies639 size_t closest(double _u, double _v)
640 {
641 double minD = 1.e22;
642 size_t I = 0;
643 for(size_t i = 0; i < nbCopies; i++) {
644 double dist = sqrt((_u - u[i]) * (_u - u[i]) + (_v - v[i]) * (_v - v[i]));
645 if(dist < minD) {
646 minD = dist;
647 I = i;
648 }
649 }
650 return id[I];
651 }
652 };
653
654 // INITIAL MESH --------- colored
655
getNodeCopies(GFace * gf,std::unordered_map<size_t,nodeCopies> & copies)656 static void getNodeCopies(GFace *gf,
657 std::unordered_map<size_t, nodeCopies> &copies)
658 {
659 std::vector<GEdge *> edges = gf->edges();
660 std::vector<GEdge *> emb_edges = gf->getEmbeddedEdges();
661 edges.insert(edges.end(), emb_edges.begin(), emb_edges.end());
662 std::set<GEdge *> touched;
663
664 if(edges.empty())
665 edges.insert(edges.end(), gf->model()->firstEdge(),
666 gf->model()->lastEdge());
667
668 for(auto e : edges) {
669 if(!e->isMeshDegenerated()) {
670 std::set<MVertex *, MVertexPtrLessThan> e_vertices;
671 for(std::size_t i = 0; i < e->lines.size(); i++) {
672 MVertex *v1 = e->lines[i]->getVertex(0);
673 MVertex *v2 = e->lines[i]->getVertex(1);
674 e_vertices.insert(v1);
675 e_vertices.insert(v2);
676 }
677 int direction = -1;
678 if(e->isSeam(gf)) {
679 direction = 0;
680 if(touched.find(e) == touched.end())
681 touched.insert(e);
682 else
683 direction = 1;
684 }
685 // printf("model edge %lu %lu vertices\n", e->tag(), e_vertices.size());
686 for(auto v : e_vertices) {
687 SPoint2 param;
688 if(direction != -1) {
689 double t = 0;
690 if(v->onWhat()->dim() == 0)
691 reparamMeshVertexOnEdge(v, e, t);
692 else if(v->onWhat()->dim() == 1)
693 v->getParameter(0, t);
694 else
695 Msg::Error("a seam edge without CAD ?");
696 param = e->reparamOnFace(gf, t, direction);
697 }
698 else {
699 // Hmm...
700 if(!gf->haveParametrization() &&
701 gf->geomType() == GEntity::DiscreteSurface) {
702 param = SPoint2(v->x(), v->y());
703 }
704 else
705 reparamMeshVertexOnFace(v, gf, param);
706 }
707 std::unordered_map<size_t, nodeCopies>::iterator it =
708 copies.find(v->getNum());
709 if(it == copies.end()) {
710 nodeCopies c(v, param.x(), param.y());
711 copies.insert(std::make_pair(v->getNum(), c));
712 }
713 else {
714 it->second.addCopy(param.x(), param.y());
715 }
716 }
717 }
718 }
719
720 std::vector<GVertex *> emb_vertx = gf->getEmbeddedVertices();
721 for(auto v : emb_vertx) {
722 SPoint2 param;
723 reparamMeshVertexOnFace(v->mesh_vertices[0], gf, param);
724 nodeCopies c(v->mesh_vertices[0], param.x(), param.y());
725 copies.insert(std::make_pair(v->mesh_vertices[0]->getNum(), c));
726 }
727 }
728
addPoints(PolyMesh * pm,std::vector<double> & pts,SBoundingBox3d & bb)729 void addPoints(PolyMesh *pm, std::vector<double> &pts, SBoundingBox3d &bb)
730 {
731 const size_t N = pts.size() / 2;
732 std::vector<double> X(N), Y(N);
733 std::vector<size_t> HC(N), IND(N);
734 PolyMesh::Face *f = pm->faces[0];
735 for(size_t i = 0; i < N; i++) {
736 X[i] = pts[2 * i];
737 Y[i] = pts[2 * i + 1];
738 HC[i] = HilbertCoordinates(X[i], Y[i], bb.center().x(), bb.center().y(),
739 bb.max().x() - bb.center().x(), 0, 0,
740 bb.max().y() - bb.center().y());
741 IND[i] = i;
742 }
743 std::sort(IND.begin(), IND.end(),
744 [&](size_t i, size_t j) { return HC[i] < HC[j]; });
745
746 for(size_t i = 0; i < N; i++) {
747 size_t I = IND[i];
748 f = Walk(f, X[I], Y[I]);
749 pm->split_triangle(i, X[I], Y[I], 0, f, delaunayEdgeCriterionPlaneIsotropic,
750 nullptr);
751 }
752 }
753
GFaceInitialMesh(int faceTag,int recover,std::vector<double> * additional)754 PolyMesh *GFaceInitialMesh(int faceTag, int recover,
755 std::vector<double> *additional)
756 {
757 GFace *gf = GModel::current()->getFaceByTag(faceTag);
758
759 if(!gf) Msg::Error("GFaceInitialMesh: no face with tag %d", faceTag);
760
761 PolyMesh *pm = new PolyMesh;
762
763 std::unordered_map<size_t, nodeCopies> copies;
764 getNodeCopies(gf, copies);
765
766 SBoundingBox3d bb;
767 for(auto c : copies) {
768 for(size_t i = 0; i < c.second.nbCopies; i++)
769 bb += SPoint3(c.second.u[i], c.second.v[i], 0);
770 }
771 bb *= 1.1;
772 pm->initialize_rectangle(bb.min().x(), bb.max().x(), bb.min().y(),
773 bb.max().y());
774 PolyMesh::Face *f = pm->faces[0];
775 for(std::unordered_map<size_t, nodeCopies>::iterator it = copies.begin();
776 it != copies.end(); ++it) {
777 for(size_t i = 0; i < it->second.nbCopies; i++) {
778 double x = it->second.u[i];
779 double y = it->second.v[i];
780 // find face in which lies x,y
781 f = Walk(f, x, y);
782 // split f and then swap edges to recover delaunayness
783 pm->split_triangle(-1, x, y, 0, f, delaunayEdgeCriterionPlaneIsotropic,
784 nullptr);
785 // remember node tags
786 it->second.id[i] = pm->vertices.size() - 1;
787 pm->vertices[pm->vertices.size() - 1]->data = it->first;
788 }
789 }
790
791 //pm->print4debug(faceTag);
792
793 if(recover) {
794 std::vector<GEdge *> edges = gf->edges();
795 std::vector<GEdge *> emb_edges = gf->getEmbeddedEdges();
796 edges.insert(edges.end(), emb_edges.begin(), emb_edges.end());
797 if(edges.empty())
798 edges.insert(edges.end(), gf->model()->firstEdge(),
799 gf->model()->lastEdge());
800
801 for(auto e : edges) {
802 if(!e->isMeshDegenerated()) {
803 for(auto l : e->lines) {
804 auto c0 = copies.find(l->getVertex(0)->getNum());
805 auto c1 = copies.find(l->getVertex(1)->getNum());
806 if(c0 == copies.end() || c1 == copies.end())
807 Msg::Error("unable to find %lu %lu %d %d",
808 l->getVertex(0)->getNum(), l->getVertex(1)->getNum(),
809 c0 == copies.end(), c1 == copies.end());
810 if(c0->second.nbCopies > c1->second.nbCopies) {
811 auto cc = c0;
812 c0 = c1;
813 c1 = cc;
814 }
815 for(size_t j = 0; j < c0->second.nbCopies; j++) {
816 PolyMesh::Vertex *v0 = pm->vertices[c0->second.id[j]];
817 PolyMesh::Vertex *v1 = pm->vertices[c1->second.closest(
818 c0->second.u[j], c0->second.v[j])];
819 int result = recover_edge(pm, v0, v1);
820 if(result < 0) {
821 Msg::Warning("Impossible to recover edge %lu %lu (error tag %d)",
822 l->getVertex(0)->getNum(), l->getVertex(0)->getNum(),
823 result);
824 }
825 else {
826 PolyMesh::HalfEdge *he = pm->getEdge(v0, v1);
827 if(he) {
828 if(he->opposite) he->opposite->data = e->tag();
829 he->data = e->tag();
830 }
831 }
832 }
833 }
834 }
835 }
836
837 // color all PolyMesh::Faces
838 // the first 4 vertices are "infinite vertices" --> color them with tag -2
839 // meaning exterior
840 PolyMesh::HalfEdge *other_side = Color(pm->vertices[0]->he, -2);
841 // other_side is inthernal to the face --> color them with tag faceTag
842 other_side = Color(other_side, faceTag);
843 // holes will be tagged -1
844
845 // flip edges that have been scrambled
846 int iter = 0;
847 while(iter++ < 100) {
848 int count = 0;
849 for(auto he : pm->hedges) {
850 if(he->opposite && he->f->data == faceTag &&
851 he->opposite->f->data == faceTag) {
852 if(delaunayEdgeCriterionPlaneIsotropic(he, nullptr)) {
853 if(intersect(he->v, he->next->v, he->next->next->v,
854 he->opposite->next->next->v)) {
855 pm->swap_edge(he);
856 count++;
857 }
858 }
859 }
860 }
861 if(!count) break;
862 }
863 }
864
865 if(additional) addPoints(pm, *additional, bb);
866
867 return pm;
868 }
869