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 <stdio.h>
7 #include <string.h>
8 #include <assert.h>
9 #include <sstream>
10 #include "GmshConfig.h"
11 #include "meshGRegionBoundaryRecovery.h"
12
13 #if defined(HAVE_TETGENBR)
14
15 #include "meshGRegion.h"
16 #include "meshGRegionDelaunayInsertion.h"
17 #include "robustPredicates.h"
18 #include "GModel.h"
19 #include "GRegion.h"
20 #include "GFace.h"
21 #include "MVertex.h"
22 #include "MLine.h"
23 #include "MPoint.h"
24 #include "MTriangle.h"
25 #include "MQuadrangle.h"
26 #include "MTetrahedron.h"
27 #include "Context.h"
28 #include "OS.h"
29 #if !defined(HAVE_NO_STDINT_H)
30 #include <stdint.h>
31 #elif defined(HAVE_NO_INTPTR_T)
32 typedef unsigned long intptr_t;
33 #endif
34
35 #if defined(HAVE_POST)
36 #include "PView.h"
37 #endif
38
39 #if defined(HAVE_FLTK)
40 #include "FlGui.h"
41 #include "drawContext.h"
42 #endif
43
44 #if 0
45 static int computeTetGenVersion2(uint32_t v1, uint32_t* v2Choices,
46 const int iface2)
47 {
48 int i;
49 for (i = 0; i < 3; i++) {
50 if(v1 == v2Choices[i]){
51 break;
52 }
53 }
54
55 if(i == 3)
56 Msg::Error("Should never happen (file:%s line:%d)", __FILE__, __LINE__);
57
58 // version%4 : corresponding face in adjacent tet
59 // version/4 : which of the 3 rotation of the facet the tetrahedra has...
60 return 4 * i + iface2;
61 }
62 #endif
63
64 namespace tetgenBR {
65
66 #define REAL double
67
68 struct brdata {
69 GRegion *gr;
70 splitQuadRecovery *sqr;
71 };
72
73 // dummy tetgenio class
74 class tetgenio {
75 public:
76 int firstnumber;
77 int numberofpointattributes;
78 int numberoftetrahedronattributes;
79 int numberofsegmentconstraints;
80 REAL *segmentconstraintlist;
81 int numberoffacetconstraints;
82 REAL *facetconstraintlist;
83 int numberofpoints;
84 int *pointlist;
85 int *pointattributelist;
86 int numberofpointmakers;
87 int *pointmarkerlist;
88 int numberofpointmtrs;
89 int *pointmtrlist;
90 int numberofedges;
91 int *edgelist;
92 int *edgemarkerlist;
93 int numberofholes;
94 REAL *holelist;
95 int numberofregions;
96 REAL *regionlist;
97 int mesh_dim;
tetgenio()98 tetgenio()
99 {
100 firstnumber = 1;
101 numberofpointattributes = 0;
102 numberoftetrahedronattributes = 0;
103 numberofsegmentconstraints = 0;
104 segmentconstraintlist = nullptr;
105 numberoffacetconstraints = 0;
106 facetconstraintlist = nullptr;
107 numberofpoints = 0;
108 pointlist = nullptr;
109 pointattributelist = nullptr;
110 numberofpointmakers = 0;
111 pointmarkerlist = nullptr;
112 numberofpointmtrs = 0;
113 pointmtrlist = nullptr;
114 numberofedges = 0;
115 edgelist = nullptr;
116 edgemarkerlist = nullptr;
117 numberofholes = 0;
118 holelist = nullptr;
119 numberofregions = 0;
120 regionlist = nullptr;
121 mesh_dim = 0;
122 }
123 };
124
125 // redefinition of predicates using our own
126 #define orient3d robustPredicates::orient3d
127 #define insphere robustPredicates::insphere
orient4d(double *,double *,double *,double *,double *,double,double,double,double,double)128 static double orient4d(double *, double *, double *, double *, double *,
129 double, double, double, double, double)
130 {
131 return 0.;
132 }
clock()133 static int clock() { return 0; }
134 #define clock_t int
135 #if !defined(TETLIBRARY)
136 #define TETLIBRARY
137 #endif
138 #define printf Msg::Auto
139 #include "tetgenBR.h"
140 #include "tetgenBR.cxx"
141 #undef printf
142
reconstructmesh(void * p)143 bool tetgenmesh::reconstructmesh(void *p)
144 {
145 GRegion *_gr = ((brdata *)p)->gr;
146 splitQuadRecovery *_sqr = ((brdata *)p)->sqr;
147
148 char opts[128];
149 sprintf(opts, "YpeQT%gp/%g", CTX::instance()->mesh.toleranceInitialDelaunay,
150 CTX::instance()->mesh.angleToleranceFacetOverlap);
151 b->parse_commandline(opts);
152
153 double t_start = Cpu(), w_start = TimeOfDay();
154 std::vector<MVertex *> _vertices;
155 std::map<int, MVertex *> _extras;
156 // Get the set of vertices from GRegion.
157 {
158 std::set<MVertex *, MVertexPtrLessThan> all;
159 std::vector<GFace *> const &f = _gr->faces();
160 for(auto it = f.begin(); it != f.end(); ++it) {
161 GFace *gf = *it;
162 for(std::size_t i = 0; i < gf->triangles.size(); i++) {
163 MVertex *v0 = gf->triangles[i]->getVertex(0);
164 MVertex *v1 = gf->triangles[i]->getVertex(1);
165 MVertex *v2 = gf->triangles[i]->getVertex(2);
166 all.insert(v0);
167 all.insert(v1);
168 all.insert(v2);
169 }
170 if(_sqr) {
171 for(std::size_t i = 0; i < gf->quadrangles.size(); i++) {
172 MVertex *v0 = gf->quadrangles[i]->getVertex(0);
173 MVertex *v1 = gf->quadrangles[i]->getVertex(1);
174 MVertex *v2 = gf->quadrangles[i]->getVertex(2);
175 MVertex *v3 = gf->quadrangles[i]->getVertex(3);
176 MVertex *newv =
177 new MVertex((v0->x() + v1->x() + v2->x() + v3->x()) * 0.25,
178 (v0->y() + v1->y() + v2->y() + v3->y()) * 0.25,
179 (v0->z() + v1->z() + v2->z() + v3->z()) * 0.25, gf);
180 // the extra vertex will be added in a GRegion (and reclassified
181 // correctly on that GRegion) when the pyramid is generated
182 MFace mf = gf->quadrangles[i]->getFace(0);
183 _sqr->add(mf, newv, gf);
184 all.insert(v0);
185 all.insert(v1);
186 all.insert(v2);
187 all.insert(v3);
188 all.insert(newv);
189 }
190 }
191 }
192 std::vector<GEdge *> const &e = _gr->embeddedEdges();
193 for(auto it = e.begin(); it != e.end(); ++it) {
194 GEdge *ge = *it;
195 for(std::size_t i = 0; i < ge->lines.size(); i++) {
196 all.insert(ge->lines[i]->getVertex(0));
197 all.insert(ge->lines[i]->getVertex(1));
198 }
199 }
200 std::vector<GVertex *> const &v = _gr->embeddedVertices();
201 for(auto it = v.begin(); it != v.end(); ++it) {
202 GVertex *gv = *it;
203 for(std::size_t i = 0; i < gv->points.size(); i++) {
204 all.insert(gv->points[i]->getVertex(0));
205 }
206 }
207 all.insert(_gr->mesh_vertices.begin(), _gr->mesh_vertices.end());
208
209 _vertices.insert(_vertices.begin(), all.begin(), all.end());
210 }
211
212 initializepools();
213
214 // Store all coordinates of the vertices as these will be pertubated in
215 // function delaunayTriangulation
216 std::map<MVertex *, SPoint3> originalCoordinates;
217 for(std::size_t i = 0; i < _vertices.size(); i++) {
218 MVertex *v = _vertices[i];
219 originalCoordinates[v] = v->point();
220 }
221
222 std::vector<MTetrahedron *> tets;
223
224 delaunayMeshIn3D(_vertices,
225 tets); // will add 8 MVertices at the end of _vertices
226 if(Msg::GetErrorCount()) return false;
227
228 Msg::Debug("Points have been tetrahedralized");
229
230 {
231 point pointloop;
232 REAL x, y, z;
233
234 // Read the points.
235 for(std::size_t i = 0; i < _vertices.size(); i++) {
236 makepoint(&pointloop, UNUSEDVERTEX);
237 // Read the point coordinates.
238 x = pointloop[0] = _vertices[i]->x();
239 y = pointloop[1] = _vertices[i]->y();
240 z = pointloop[2] = _vertices[i]->z();
241 // Determine the smallest and largest x, y and z coordinates.
242 if(i == 0) {
243 xmin = xmax = x;
244 ymin = ymax = y;
245 zmin = zmax = z;
246 }
247 else {
248 xmin = (x < xmin) ? x : xmin;
249 xmax = (x > xmax) ? x : xmax;
250 ymin = (y < ymin) ? y : ymin;
251 ymax = (y > ymax) ? y : ymax;
252 zmin = (z < zmin) ? z : zmin;
253 zmax = (z > zmax) ? z : zmax;
254 }
255 }
256
257 // 'longest' is the largest possible edge length formed by input vertices.
258 x = xmax - xmin;
259 y = ymax - ymin;
260 z = zmax - zmin;
261 longest = sqrt(x * x + y * y + z * z);
262 if(longest == 0.0) {
263 Msg::Warning("The point set is trivial");
264 return true;
265 }
266
267 // Two identical points are distinguished by 'lengthlimit'.
268 if(minedgelength == 0.0) { minedgelength = longest * b->epsilon; }
269 }
270
271 point *idx2verlist;
272
273 // Create a map from indices to vertices.
274 makeindex2pointmap(idx2verlist);
275 // 'idx2verlist' has length 'in->numberofpoints + 1'.
276 idx2verlist[0] = dummypoint; // Let 0th-entry be dummypoint.
277 // Index the vertices, starting at 1 (vertex index 0 is used as special code
278 // in tetgenBR in case of failure)
279 for(std::size_t i = 0; i < _vertices.size(); i++) {
280 _vertices[i]->setIndex(i + 1);
281 }
282
283 {
284 tetrahedron *ver2tetarray;
285 triface tetloop, checktet, prevchktet;
286 triface hulltet, face1, face2;
287 tetrahedron tptr;
288 point p[4], q[3];
289 REAL ori; //, attrib, volume;
290 int bondflag;
291 int t1ver;
292 int idx, k;
293
294 Msg::Info("Reconstructing mesh...");
295
296 // Allocate an array that maps each vertex to its adjacent tets.
297 ver2tetarray = new tetrahedron[_vertices.size() + in->firstnumber];
298 for(std::size_t i = 0; i < _vertices.size() + in->firstnumber; i++) {
299 setpointtype(idx2verlist[i], VOLVERTEX); // initial type.
300 ver2tetarray[i] = nullptr;
301 }
302
303 #if 0
304 /* N E W V E R S I O N */
305 std::vector<triface> ts( tets.size() );
306 for(std::size_t i = 0; i < tets.size(); i++) {
307 point p[4];
308 // index tetrahedra in order to have access to neighbors ids.
309 tets[i]->tet()->forceNum(i+1);
310 p[0] = idx2verlist[tets[i]->getVertex(0)->getIndex()];
311 p[1] = idx2verlist[tets[i]->getVertex(1)->getIndex()];
312 p[2] = idx2verlist[tets[i]->getVertex(2)->getIndex()];
313 p[3] = idx2verlist[tets[i]->getVertex(3)->getIndex()];
314 setvertices(ts[i], p[0], p[1], p[2], p[3]);
315 }
316 // we can make this in parallel, iterations are totally independent
317 for (uint64_t i = 0; i < tets.size(); i++) {
318 triface tf1 = ts[i];
319
320 for (tf1.ver=0; tf1.ver<4; tf1.ver++){
321 uint64_t neigh = tets[i]->getNeigh(tf1.ver)->tet()->getNum() - 1;
322 triface tf2 = ts[neigh];
323 int iface2 = tf1.ver;
324
325 int face2[3] = {
326 tets[i]->getVertex(faces_tetra(tf1.ver),0)->getIndex(),
327 tets[i]->getVertex(faces_tetra(tf1.ver),1)->getIndex(),
328 tets[i]->getVertex(faces_tetra(tf1.ver),2)->getIndex()};
329
330 tf2.ver = computeTetGenVersion2(faces2[0], face2, iface2);
331 bond(tf1,tf2);
332 }
333 }
334
335 #else
336
337 /* N E W V E R S I O N */
338
339 // Create the tetrahedra and connect those that share a common face.
340 for(std::size_t i = 0; i < tets.size(); i++) {
341 // Get the four vertices.
342 for(int j = 0; j < 4; j++) {
343 p[j] = idx2verlist[tets[i]->getVertex(j)->getIndex()];
344 }
345 // Check the orientation.
346 ori = orient3d(p[0], p[1], p[2], p[3]);
347 if(ori > 0.0) {
348 // Swap the first two vertices.
349 q[0] = p[0];
350 p[0] = p[1];
351 p[1] = q[0];
352 }
353 else if(ori == 0.0) {
354 if(!b->quiet) {
355 Msg::Warning("Tet #%d is degenerated", i + in->firstnumber);
356 }
357 }
358 // Create a new tetrahedron.
359 maketetrahedron(&tetloop); // tetloop.ver = 11.
360 setvertices(tetloop, p[0], p[1], p[2], p[3]);
361 // Try connecting this tet to others that share the common faces.
362 for(tetloop.ver = 0; tetloop.ver < 4; tetloop.ver++) {
363 p[3] = oppo(tetloop);
364 // Look for other tets having this vertex.
365 idx = pointmark(p[3]) - in->firstnumber;
366 tptr = ver2tetarray[idx];
367 // Link the current tet to the next one in the stack.
368 tetloop.tet[8 + tetloop.ver] = tptr;
369 // Push the current tet onto the stack.
370 ver2tetarray[idx] = encode(tetloop);
371 decode(tptr, checktet);
372 if(checktet.tet != nullptr) {
373 p[0] = org(tetloop); // a
374 p[1] = dest(tetloop); // b
375 p[2] = apex(tetloop); // c
376 prevchktet = tetloop;
377 do {
378 q[0] = org(checktet); // a'
379 q[1] = dest(checktet); // b'
380 q[2] = apex(checktet); // c'
381 // Check the three faces at 'd' in 'checktet'.
382 bondflag = 0;
383 for(int j = 0; j < 3; j++) {
384 // Go to the face [b',a',d], or [c',b',d], or [a',c',d].
385 esym(checktet, face2);
386 if(face2.tet[face2.ver & 3] == nullptr) {
387 k = ((j + 1) % 3);
388 if(q[k] == p[0]) { // b', c', a' = a
389 if(q[j] == p[1]) { // a', b', c' = b
390 // [#,#,d] is matched to [b,a,d].
391 esym(tetloop, face1);
392 bond(face1, face2);
393 bondflag++;
394 }
395 }
396 if(q[k] == p[1]) { // b',c',a' = b
397 if(q[j] == p[2]) { // a',b',c' = c
398 // [#,#,d] is matched to [c,b,d].
399 enext(tetloop, face1);
400 esymself(face1);
401 bond(face1, face2);
402 bondflag++;
403 }
404 }
405 if(q[k] == p[2]) { // b',c',a' = c
406 if(q[j] == p[0]) { // a',b',c' = a
407 // [#,#,d] is matched to [a,c,d].
408 eprev(tetloop, face1);
409 esymself(face1);
410 bond(face1, face2);
411 bondflag++;
412 }
413 }
414 }
415 else {
416 bondflag++;
417 }
418 enextself(checktet);
419 } // j
420 // Go to the next tet in the link.
421 tptr = checktet.tet[8 + checktet.ver];
422 if(bondflag == 3) {
423 // All three faces at d in 'checktet' have been connected.
424 // It can be removed from the link.
425 prevchktet.tet[8 + prevchktet.ver] = tptr;
426 }
427 else {
428 // Bakup the previous tet in the link.
429 prevchktet = checktet;
430 }
431 decode(tptr, checktet);
432 } while(checktet.tet != nullptr);
433 } // if(checktet.tet != nullptr)
434 } // for(tetloop.ver = 0; ...
435 } // i
436
437 // Remember a tet of the mesh.
438 recenttet = tetloop;
439
440 // Create hull tets, create the point-to-tet map, and clean up the
441 // temporary spaces used in each tet.
442 hullsize = tetrahedrons->items;
443
444 tetrahedrons->traversalinit();
445 tetloop.tet = tetrahedrontraverse();
446 while(tetloop.tet != (tetrahedron *)nullptr) {
447 tptr = encode(tetloop);
448 for(tetloop.ver = 0; tetloop.ver < 4; tetloop.ver++) {
449 if(tetloop.tet[tetloop.ver] == nullptr) {
450 // Create a hull tet.
451 maketetrahedron(&hulltet);
452 p[0] = org(tetloop);
453 p[1] = dest(tetloop);
454 p[2] = apex(tetloop);
455 setvertices(hulltet, p[1], p[0], p[2], dummypoint);
456 bond(tetloop, hulltet);
457 // Try connecting this to others that share common hull edges.
458 for(int j = 0; j < 3; j++) {
459 fsym(hulltet, face2);
460 while(1) {
461 if(face2.tet == nullptr) break;
462 esymself(face2);
463 if(apex(face2) == dummypoint) break;
464 fsymself(face2);
465 }
466 if(face2.tet != nullptr) {
467 // Found an adjacent hull tet.
468 assert(face2.tet[face2.ver & 3] == nullptr);
469 esym(hulltet, face1);
470 bond(face1, face2);
471 }
472 enextself(hulltet);
473 }
474 // hullsize++;
475 }
476 // Create the point-to-tet map.
477 setpoint2tet((point)(tetloop.tet[4 + tetloop.ver]), tptr);
478 // Clean the temporary used space.
479 tetloop.tet[8 + tetloop.ver] = nullptr;
480 }
481 tetloop.tet = tetrahedrontraverse();
482 }
483
484 hullsize = tetrahedrons->items - hullsize;
485
486 delete[] ver2tetarray;
487 for(std::size_t i = 0; i < tets.size(); i++) delete tets[i];
488 tets.clear(); // Release all memory in this vector.
489 }
490 #endif
491
492 std::vector<GFace *> const &f_list = _gr->faces();
493 std::vector<GEdge *> const &e_list = _gr->embeddedEdges();
494
495 {
496 Msg::Info(" - Creating surface mesh");
497 face newsh;
498 face newseg;
499 point p[4];
500 int idx;
501
502 for(auto it = f_list.begin(); it != f_list.end(); ++it) {
503 GFace *gf = *it;
504 for(std::size_t i = 0; i < gf->triangles.size(); i++) {
505 for(int j = 0; j < 3; j++) {
506 p[j] = idx2verlist[gf->triangles[i]->getVertex(j)->getIndex()];
507 if(pointtype(p[j]) == VOLVERTEX) {
508 setpointtype(p[j], FACETVERTEX);
509 }
510 }
511 makeshellface(subfaces, &newsh);
512 setshvertices(newsh, p[0], p[1], p[2]);
513 setshellmark(newsh, gf->tag()); // the GFace's tag.
514 recentsh = newsh;
515 for(int j = 0; j < 3; j++) {
516 makeshellface(subsegs, &newseg);
517 setshvertices(newseg, sorg(newsh), sdest(newsh), nullptr);
518 // Set the default segment marker '-1'.
519 setshellmark(newseg, -1);
520 ssbond(newsh, newseg);
521 senextself(newsh);
522 }
523 }
524 } // it
525
526 if(_sqr) {
527 std::map<MFace, GFace *, MFaceLessThan> f = _sqr->getTri();
528 for(auto it = f.begin(); it != f.end(); it++) {
529 const MFace &mf = it->first;
530 for(int j = 0; j < 3; j++) {
531 p[j] = idx2verlist[mf.getVertex(j)->getIndex()];
532 if(pointtype(p[j]) == VOLVERTEX) {
533 setpointtype(p[j], FACETVERTEX);
534 }
535 }
536 makeshellface(subfaces, &newsh);
537 setshvertices(newsh, p[0], p[1], p[2]);
538 setshellmark(newsh, it->second->tag());
539 recentsh = newsh;
540 for(int j = 0; j < 3; j++) {
541 makeshellface(subsegs, &newseg);
542 setshvertices(newseg, sorg(newsh), sdest(newsh), nullptr);
543 // Set the default segment marker '-1'.
544 setshellmark(newseg, -1);
545 ssbond(newsh, newseg);
546 senextself(newsh);
547 }
548 }
549 }
550
551 // Connecting triangles, removing redundant segments.
552 unifysegments();
553
554 Msg::Info(" - Identifying boundary edges");
555
556 face *shperverlist;
557 int *idx2shlist;
558 face searchsh, neighsh;
559 face segloop, checkseg;
560 point checkpt;
561
562 // Construct a map from points to subfaces.
563 makepoint2submap(subfaces, idx2shlist, shperverlist);
564
565 // Process the set of PSC edges.
566 // Remeber that all segments have default marker '-1'.
567 // int COUNTER = 0;
568 for(auto it = e_list.begin(); it != e_list.end(); ++it) {
569 GEdge *ge = *it;
570 for(std::size_t i = 0; i < ge->lines.size(); i++) {
571 for(int j = 0; j < 2; j++) {
572 p[j] = idx2verlist[ge->lines[i]->getVertex(j)->getIndex()];
573 setpointtype(p[j], RIDGEVERTEX);
574 }
575 if(p[0] == p[1]) {
576 // This is a potential problem in surface mesh.
577 continue; // Skip this edge.
578 }
579 // Find a face contains the edge p[0], p[1].
580 newseg.sh = nullptr;
581 searchsh.sh = nullptr;
582 idx = pointmark(p[0]) - in->firstnumber;
583 for(int j = idx2shlist[idx]; j < idx2shlist[idx + 1]; j++) {
584 checkpt = sdest(shperverlist[j]);
585 if(checkpt == p[1]) {
586 searchsh = shperverlist[j];
587 break; // Found.
588 }
589 else {
590 checkpt = sapex(shperverlist[j]);
591 if(checkpt == p[1]) {
592 senext2(shperverlist[j], searchsh);
593 sesymself(searchsh);
594 break;
595 }
596 }
597 } // j
598 if(searchsh.sh != nullptr) {
599 // Check if this edge is already a segment of the mesh.
600 sspivot(searchsh, checkseg);
601 if(checkseg.sh != nullptr) {
602 // This segment already exist.
603 newseg = checkseg;
604 }
605 else {
606 // Create a new segment at this edge.
607 makeshellface(subsegs, &newseg);
608 setshvertices(newseg, p[0], p[1], nullptr);
609 ssbond(searchsh, newseg);
610 spivot(searchsh, neighsh);
611 if(neighsh.sh != nullptr) { ssbond(neighsh, newseg); }
612 }
613 }
614 else {
615 // It is a dangling segment (not belong to any facets).
616 // Check if segment [p[0],p[1]] already exists.
617 // TODO: Change the brute-force search. Slow!
618 /* point *ppt;
619 subsegs->traversalinit();
620 segloop.sh = shellfacetraverse(subsegs);
621 while (segloop.sh != nullptr){
622 ppt = (point *) &(segloop.sh[3]);
623 if(((ppt[0] == p[0]) && (ppt[1] == p[1])) ||
624 ((ppt[0] == p[1]) && (ppt[1] == p[0]))){
625 // Found!
626 newseg = segloop;
627 break;
628 }
629 segloop.sh = shellfacetraverse(subsegs);
630 }*/
631 if(newseg.sh == nullptr) {
632 makeshellface(subsegs, &newseg);
633 setshvertices(newseg, p[0], p[1], nullptr);
634 }
635 }
636 setshellmark(newseg, ge->tag());
637 } // i
638 } // e_list
639
640 delete[] shperverlist;
641 delete[] idx2shlist;
642
643 Msg::Debug(" %ld (%ld) subfaces (segments)", subfaces->items,
644 subsegs->items);
645
646 // The total number of iunput segments.
647 insegments = subsegs->items;
648
649 if(0) { outmesh2medit("dump2"); }
650 }
651
652 delete[] idx2verlist;
653
654 // Boundary recovery.
655
656 clock_t t;
657 Msg::Info(" - Recovering boundary");
658 recoverboundary(t);
659
660 carveholes();
661
662 if(subvertstack->objects > 0l) { suppresssteinerpoints(); }
663
664 recoverdelaunay();
665
666 // let's try
667 optimizemesh();
668
669 if((dupverts > 0l) || (unuverts > 0l)) {
670 // Remove hanging nodes.
671 // cannot call this here due to 8 additional exterior vertices we
672 // inserted jettisonnodes();
673 }
674
675 long tetnumber, facenumber;
676
677 Msg::Debug("Statistics:\n");
678 Msg::Debug(" Input points: %ld", _vertices.size());
679 if(b->plc) {
680 Msg::Debug(" Input facets: %ld", f_list.size());
681 Msg::Debug(" Input segments: %ld", e_list.size());
682 }
683
684 tetnumber = tetrahedrons->items - hullsize;
685 facenumber = (tetnumber * 4l + hullsize) / 2l;
686
687 if(b->weighted) { // -w option
688 Msg::Debug(" Mesh points: %ld", points->items - nonregularcount);
689 }
690 else {
691 Msg::Debug(" Mesh points: %ld", points->items);
692 }
693 Msg::Debug(" Mesh tetrahedra: %ld", tetnumber);
694 Msg::Debug(" Mesh faces: %ld", facenumber);
695 if(meshedges > 0l) { Msg::Debug(" Mesh edges: %ld", meshedges); }
696 else {
697 if(!nonconvex) {
698 long vsize = points->items - dupverts - unuverts;
699 if(b->weighted) vsize -= nonregularcount;
700 meshedges = vsize + facenumber - tetnumber - 1;
701 Msg::Debug(" Mesh edges: %ld", meshedges);
702 }
703 }
704
705 if(b->plc || b->refine) {
706 Msg::Debug(" Mesh faces on facets: %ld", subfaces->items);
707 Msg::Debug(" Mesh edges on segments: %ld", subsegs->items);
708 if(st_volref_count > 0l) {
709 Msg::Debug(" Steiner points inside domain: %ld", st_volref_count);
710 }
711 if(st_facref_count > 0l) {
712 Msg::Debug(" Steiner points on facets: %ld", st_facref_count);
713 }
714 if(st_segref_count > 0l) {
715 Msg::Debug(" Steiner points on segments: %ld", st_segref_count);
716 }
717 }
718 else {
719 Msg::Debug(" Convex hull faces: %ld", hullsize);
720 if(meshhulledges > 0l) {
721 Msg::Debug(" Convex hull edges: %ld", meshhulledges);
722 }
723 }
724 if(b->weighted) { // -w option
725 Msg::Debug(" Skipped non-regular points: %ld", nonregularcount);
726 }
727
728 // Debug
729 if(0) { outmesh2medit("dump"); }
730
731 {
732 // Write mesh into to GRegion.
733
734 Msg::Debug("Writing to GRegion...");
735
736 point p[4];
737
738 // In some hard cases, the surface mesh may be modified.
739 // Find the list of GFaces, GEdges that have been modified.
740 std::set<int> l_faces, l_edges;
741
742 if(points->items > (int)_vertices.size()) {
743 face parentseg, parentsh, spinsh;
744 point pointloop;
745 // Create newly added mesh vertices.
746 // The new vertices must be added at the end of the point list.
747 points->traversalinit();
748 pointloop = pointtraverse();
749 while(pointloop != (point)nullptr) {
750 if(issteinerpoint(pointloop)) {
751 // Check if this Steiner point locates on boundary.
752 if(pointtype(pointloop) == FREESEGVERTEX) {
753 sdecode(point2sh(pointloop), parentseg);
754 assert(parentseg.sh != nullptr);
755 l_edges.insert(shellmark(parentseg));
756 // Get the GEdge containing this vertex.
757 GEdge *ge = nullptr;
758 GFace *gf = nullptr;
759 int etag = shellmark(parentseg);
760 for(auto it = e_list.begin(); it != e_list.end(); ++it) {
761 if((*it)->tag() == etag) {
762 ge = *it;
763 break;
764 }
765 }
766 if(ge != nullptr) {
767 MEdgeVertex *v = new MEdgeVertex(pointloop[0], pointloop[1],
768 pointloop[2], ge, 0);
769 double uu = 0;
770 if(reparamMeshVertexOnEdge(v, ge, uu)) {
771 v->setParameter(0, uu);
772 }
773 v->setIndex(pointmark(pointloop));
774 _gr->mesh_vertices.push_back(v);
775 _extras[pointmark(pointloop) - in->firstnumber] = v;
776 }
777 spivot(parentseg, parentsh);
778 if(parentsh.sh != nullptr) {
779 if(ge == nullptr) {
780 // We treat this vertex a facet vertex.
781 int ftag = shellmark(parentsh);
782 for(auto it = f_list.begin(); it != f_list.end(); ++it) {
783 if((*it)->tag() == ftag) {
784 gf = *it;
785 break;
786 }
787 }
788 if(gf != nullptr) {
789 MFaceVertex *v = new MFaceVertex(
790 pointloop[0], pointloop[1], pointloop[2], gf, 0, 0);
791 SPoint2 param;
792 if(reparamMeshVertexOnFace(v, gf, param)) {
793 v->setParameter(0, param.x());
794 v->setParameter(1, param.y());
795 }
796 v->setIndex(pointmark(pointloop));
797 _gr->mesh_vertices.push_back(v);
798 _extras[pointmark(pointloop) - in->firstnumber] = v;
799 }
800 }
801 // Record all the GFaces' tag at this segment.
802 spinsh = parentsh;
803 while(1) {
804 l_faces.insert(shellmark(spinsh));
805 spivotself(spinsh);
806 if(spinsh.sh == parentsh.sh) break;
807 }
808 }
809 if((ge == nullptr) && (gf == nullptr)) {
810 // Create an interior mesh vertex.
811 MVertex *v =
812 new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
813 v->setIndex(pointmark(pointloop));
814 _extras[pointmark(pointloop) - in->firstnumber] = v;
815 _gr->mesh_vertices.push_back(v);
816 }
817 }
818 else if(pointtype(pointloop) == FREEFACETVERTEX) {
819 sdecode(point2sh(pointloop), parentsh);
820 assert(parentsh.sh != nullptr);
821 l_faces.insert(shellmark(parentsh));
822 // Get the GFace containing this vertex.
823 GFace *gf = nullptr;
824 int ftag = shellmark(parentsh);
825 for(auto it = f_list.begin(); it != f_list.end(); ++it) {
826 if((*it)->tag() == ftag) {
827 gf = *it;
828 break;
829 }
830 }
831 if(gf != nullptr) {
832 MFaceVertex *v = new MFaceVertex(pointloop[0], pointloop[1],
833 pointloop[2], gf, 0, 0);
834 SPoint2 param;
835 if(reparamMeshVertexOnFace(v, gf, param)) {
836 v->setParameter(0, param.x());
837 v->setParameter(1, param.y());
838 }
839 v->setIndex(pointmark(pointloop));
840 _gr->mesh_vertices.push_back(v);
841 _extras[pointmark(pointloop) - in->firstnumber] = v;
842 }
843 else {
844 // Create a mesh vertex.
845 MVertex *v =
846 new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
847 v->setIndex(pointmark(pointloop));
848 _gr->mesh_vertices.push_back(v);
849 _extras[pointmark(pointloop) - in->firstnumber] = v;
850 }
851 }
852 else {
853 MVertex *v =
854 new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
855 v->setIndex(pointmark(pointloop));
856 _gr->mesh_vertices.push_back(v);
857 _extras[pointmark(pointloop) - in->firstnumber] = v;
858 }
859 }
860 pointloop = pointtraverse();
861 }
862 // assert((int)_vertices.size() == points->items);
863 }
864
865 if(!_extras.empty())
866 Msg::Info(" - Added %d Steiner point%s", _extras.size(),
867 (_extras.size() > 1) ? "s" : "");
868
869 if(l_edges.size() > 0) {
870 // There are Steiner points on segments!
871 face segloop;
872 // Re-create the segment mesh in the corresponding GEdges.
873 for(auto it = l_edges.begin(); it != l_edges.end(); ++it) {
874 // Find the GEdge with tag = *it.
875
876 int etag = *it;
877 GEdge *ge = nullptr;
878 for(auto it = e_list.begin(); it != e_list.end(); ++it) {
879 if((*it)->tag() == etag) {
880 ge = *it;
881 break;
882 }
883 }
884 if(ge != nullptr) {
885 Msg::Info(" - Steiner points exist on curve %d", ge->tag());
886 // Delete the old triangles.
887 for(std::size_t i = 0; i < ge->lines.size(); i++)
888 delete ge->lines[i];
889 ge->lines.clear();
890 ge->deleteVertexArrays();
891 // Create the new triangles.
892 segloop.shver = 0;
893 subsegs->traversalinit();
894 segloop.sh = shellfacetraverse(subsegs);
895 while(segloop.sh != nullptr) {
896 if(shellmark(segloop) == etag) {
897 p[0] = sorg(segloop);
898 p[1] = sdest(segloop);
899 int idx1 = pointmark(p[0]) - in->firstnumber;
900 MVertex *v1 = idx1 >= (int)_vertices.size() ? _extras[idx1] :
901 _vertices[idx1];
902 int idx2 = pointmark(p[1]) - in->firstnumber;
903 MVertex *v2 = idx2 >= (int)_vertices.size() ? _extras[idx2] :
904 _vertices[idx2];
905 MLine *t = new MLine(v1, v2);
906 ge->lines.push_back(t);
907 }
908 segloop.sh = shellfacetraverse(subsegs);
909 }
910 }
911 else {
912 Msg::Debug("Unknown curve %d with Steiner point(s)", etag);
913 }
914 } // it
915 }
916
917 if(l_faces.size() > 0) {
918 // There are Steiner points on facets!
919 face subloop;
920 // Re-create the surface mesh in the corresponding GFaces.
921 for(auto it = l_faces.begin(); it != l_faces.end(); ++it) {
922 // Find the GFace with tag = *it.
923
924 int ftag = *it;
925 GFace *gf = nullptr;
926 for(auto it = f_list.begin(); it != f_list.end(); ++it) {
927 if((*it)->tag() == ftag) {
928 gf = *it;
929 break;
930 }
931 }
932 if(gf != nullptr) {
933 // Delete the old triangles.
934 Msg::Info(" - Steiner points exist on surface %d", gf->tag());
935 for(std::size_t i = 0; i < gf->triangles.size(); i++)
936 delete gf->triangles[i];
937 gf->triangles.clear();
938 gf->deleteVertexArrays();
939
940 if(gf->quadrangles.size()) {
941 Msg::Warning("Steiner points not handled for quad surface mesh");
942 }
943
944 // Create the new triangles.
945 subloop.shver = 0;
946 subfaces->traversalinit();
947 subloop.sh = shellfacetraverse(subfaces);
948 while(subloop.sh != nullptr) {
949 if(shellmark(subloop) == ftag) {
950 p[0] = sorg(subloop);
951 p[1] = sdest(subloop);
952 p[2] = sapex(subloop);
953 int idx1 = pointmark(p[0]) - in->firstnumber;
954 MVertex *v1 = idx1 >= (int)_vertices.size() ? _extras[idx1] :
955 _vertices[idx1];
956 int idx2 = pointmark(p[1]) - in->firstnumber;
957 MVertex *v2 = idx2 >= (int)_vertices.size() ? _extras[idx2] :
958 _vertices[idx2];
959 int idx3 = pointmark(p[2]) - in->firstnumber;
960 MVertex *v3 = idx3 >= (int)_vertices.size() ? _extras[idx3] :
961 _vertices[idx3];
962 MTriangle *t = new MTriangle(v1, v2, v3);
963 gf->triangles.push_back(t);
964 }
965 subloop.sh = shellfacetraverse(subfaces);
966 }
967 }
968 else {
969 Msg::Debug("Unknown surface %d with Steiner point(s)", ftag);
970 }
971 } // it
972 }
973
974 triface tetloop;
975
976 tetloop.ver = 11;
977 tetrahedrons->traversalinit();
978 tetloop.tet = tetrahedrontraverse();
979
980 while(tetloop.tet != (tetrahedron *)nullptr) {
981 p[0] = org(tetloop);
982 p[1] = dest(tetloop);
983 p[2] = apex(tetloop);
984 p[3] = oppo(tetloop);
985
986 int idx1 = pointmark(p[0]) - in->firstnumber;
987 MVertex *v1 =
988 idx1 >= (int)_vertices.size() ? _extras[idx1] : _vertices[idx1];
989 int idx2 = pointmark(p[1]) - in->firstnumber;
990 MVertex *v2 =
991 idx2 >= (int)_vertices.size() ? _extras[idx2] : _vertices[idx2];
992 int idx3 = pointmark(p[2]) - in->firstnumber;
993 MVertex *v3 =
994 idx3 >= (int)_vertices.size() ? _extras[idx3] : _vertices[idx3];
995 int idx4 = pointmark(p[3]) - in->firstnumber;
996 MVertex *v4 =
997 idx4 >= (int)_vertices.size() ? _extras[idx4] : _vertices[idx4];
998 MTetrahedron *t = new MTetrahedron(v1, v2, v3, v4);
999 _gr->tetrahedra.push_back(t);
1000 tetloop.tet = tetrahedrontraverse();
1001 }
1002 } // mesh output
1003
1004 Msg::Info("Done reconstructing mesh (Wall %gs, CPU %gs)",
1005 TimeOfDay() - w_start, Cpu() - t_start);
1006
1007 // Put all coordinates back so they are not pertubated anymore
1008 // (pertubation done in delaunayTriangulation)
1009 for(auto vIter = originalCoordinates.begin();
1010 vIter != originalCoordinates.end(); ++vIter) {
1011 const SPoint3 &coordinates = vIter->second;
1012 vIter->first->setXYZ(coordinates.x(), coordinates.y(), coordinates.z());
1013 }
1014
1015 // delete 8 new enclosing box vertices added in delaunayMeshIn3d
1016 for(std::size_t i = _vertices.size() - 8; i < _vertices.size(); i++)
1017 delete _vertices[i];
1018
1019 return true;
1020 }
1021
1022 // Dump the input surface mesh.
1023 // 'mfilename' is a filename without suffix.
outsurfacemesh(const char * mfilename)1024 void tetgenmesh::outsurfacemesh(const char *mfilename)
1025 {
1026 FILE *outfile = nullptr;
1027 char sfilename[256];
1028 int firstindex;
1029
1030 point pointloop;
1031 int pointnumber;
1032 strcpy(sfilename, mfilename);
1033 strcat(sfilename, ".node");
1034 outfile = fopen(sfilename, "w");
1035 if(!b->quiet) { printf("Writing %s.\n", sfilename); }
1036 fprintf(outfile, "%ld 3 0 0\n", points->items);
1037 // Determine the first index (0 or 1).
1038 firstindex = b->zeroindex ? 0 : in->firstnumber;
1039 points->traversalinit();
1040 pointloop = pointtraverse();
1041 pointnumber = firstindex; // in->firstnumber;
1042 while(pointloop != (point)nullptr) {
1043 // Point number, x, y and z coordinates.
1044 fprintf(outfile, "%4d %.17g %.17g %.17g", pointnumber,
1045 pointloop[0], pointloop[1], pointloop[2]);
1046 fprintf(outfile, "\n");
1047 pointloop = pointtraverse();
1048 pointnumber++;
1049 }
1050 fclose(outfile);
1051
1052 face faceloop;
1053 point torg, tdest, tapex;
1054 strcpy(sfilename, mfilename);
1055 strcat(sfilename, ".smesh");
1056 outfile = fopen(sfilename, "w");
1057 if(!b->quiet) { printf("Writing %s.\n", sfilename); }
1058 int shift = 0; // Default no shiftment.
1059 if((in->firstnumber == 1) && (firstindex == 0)) {
1060 shift = 1; // Shift the output indices by 1.
1061 }
1062 fprintf(outfile, "0 3 0 0\n");
1063 fprintf(outfile, "%ld 1\n", subfaces->items);
1064 subfaces->traversalinit();
1065 faceloop.sh = shellfacetraverse(subfaces);
1066 while(faceloop.sh != (shellface *)nullptr) {
1067 torg = sorg(faceloop);
1068 tdest = sdest(faceloop);
1069 tapex = sapex(faceloop);
1070 fprintf(outfile, "3 %4d %4d %4d %d\n", pointmark(torg) - shift,
1071 pointmark(tdest) - shift, pointmark(tapex) - shift,
1072 shellmark(faceloop));
1073 faceloop.sh = shellfacetraverse(subfaces);
1074 }
1075 fprintf(outfile, "0\n");
1076 fprintf(outfile, "0\n");
1077 fclose(outfile);
1078
1079 face edgeloop;
1080 int edgenumber;
1081 strcpy(sfilename, mfilename);
1082 strcat(sfilename, ".edge");
1083 outfile = fopen(sfilename, "w");
1084 if(!b->quiet) { printf("Writing %s.\n", sfilename); }
1085 fprintf(outfile, "%ld 1\n", subsegs->items);
1086 subsegs->traversalinit();
1087 edgeloop.sh = shellfacetraverse(subsegs);
1088 edgenumber = firstindex; // in->firstnumber;
1089 while(edgeloop.sh != (shellface *)nullptr) {
1090 torg = sorg(edgeloop);
1091 tdest = sdest(edgeloop);
1092 fprintf(outfile, "%5d %4d %4d %d\n", edgenumber,
1093 pointmark(torg) - shift, pointmark(tdest) - shift,
1094 shellmark(edgeloop));
1095 edgenumber++;
1096 edgeloop.sh = shellfacetraverse(subsegs);
1097 }
1098 fclose(outfile);
1099 }
1100
outmesh2medit(const char * mfilename)1101 void tetgenmesh::outmesh2medit(const char *mfilename)
1102 {
1103 FILE *outfile;
1104 char mefilename[256];
1105 tetrahedron *tetptr;
1106 triface tface, tsymface;
1107 face segloop, checkmark;
1108 point ptloop, p1, p2, p3, p4;
1109 long ntets, faces;
1110 int shift = 0;
1111 int marker;
1112
1113 if(mfilename != (char *)nullptr && mfilename[0] != '\0') {
1114 strcpy(mefilename, mfilename);
1115 }
1116 else {
1117 strcpy(mefilename, "unnamed");
1118 }
1119 strcat(mefilename, ".mesh");
1120
1121 if(!b->quiet) { printf("Writing %s.\n", mefilename); }
1122 outfile = fopen(mefilename, "w");
1123 if(outfile == (FILE *)nullptr) {
1124 Msg::Error("Could not open file '%s'", mefilename);
1125 return;
1126 }
1127
1128 fprintf(outfile, "MeshVersionFormatted 1\n");
1129 fprintf(outfile, "\n");
1130 fprintf(outfile, "Dimension\n");
1131 fprintf(outfile, "3\n");
1132 fprintf(outfile, "\n");
1133
1134 fprintf(outfile, "\n# Set of mesh vertices\n");
1135 fprintf(outfile, "Vertices\n");
1136 fprintf(outfile, "%ld\n", points->items);
1137
1138 points->traversalinit();
1139 ptloop = pointtraverse();
1140 // pointnumber = 1;
1141 while(ptloop != (point)nullptr) {
1142 // Point coordinates.
1143 fprintf(outfile, "%.17g %.17g %.17g", ptloop[0], ptloop[1],
1144 ptloop[2]);
1145 fprintf(outfile, " 0\n");
1146 // setpointmark(ptloop, pointnumber);
1147 ptloop = pointtraverse();
1148 // pointnumber++;
1149 }
1150
1151 // Medit need start number form 1.
1152 if(in->firstnumber == 1) { shift = 0; }
1153 else {
1154 shift = 1;
1155 }
1156
1157 // Compute the number of faces.
1158 ntets = tetrahedrons->items - hullsize;
1159
1160 fprintf(outfile, "\n# Set of Tetrahedra\n");
1161 fprintf(outfile, "Tetrahedra\n");
1162 fprintf(outfile, "%ld\n", ntets);
1163
1164 tetrahedrons->traversalinit();
1165 tetptr = tetrahedrontraverse();
1166 while(tetptr != (tetrahedron *)nullptr) {
1167 if(!b->reversetetori) {
1168 p1 = (point)tetptr[4];
1169 p2 = (point)tetptr[5];
1170 }
1171 else {
1172 p1 = (point)tetptr[5];
1173 p2 = (point)tetptr[4];
1174 }
1175 p3 = (point)tetptr[6];
1176 p4 = (point)tetptr[7];
1177 fprintf(outfile, "%5d %5d %5d %5d", pointmark(p1) + shift,
1178 pointmark(p2) + shift, pointmark(p3) + shift,
1179 pointmark(p4) + shift);
1180 if(numelemattrib > 0) {
1181 fprintf(outfile, " %.17g", elemattribute(tetptr, 0));
1182 }
1183 else {
1184 fprintf(outfile, " 0");
1185 }
1186 fprintf(outfile, "\n");
1187 tetptr = tetrahedrontraverse();
1188 }
1189
1190 // faces = (ntets * 4l + hullsize) / 2l;
1191 faces = subfaces->items;
1192 face sface;
1193
1194 fprintf(outfile, "\n# Set of Triangles\n");
1195 fprintf(outfile, "Triangles\n");
1196 fprintf(outfile, "%ld\n", faces);
1197
1198 subfaces->traversalinit();
1199 sface.sh = shellfacetraverse(subfaces);
1200 while(sface.sh != nullptr) {
1201 p1 = sorg(sface);
1202 p2 = sdest(sface);
1203 p3 = sapex(sface);
1204 fprintf(outfile, "%5d %5d %5d", pointmark(p1) + shift,
1205 pointmark(p2) + shift, pointmark(p3) + shift);
1206 marker = shellmark(sface);
1207 fprintf(outfile, " %d\n", marker);
1208 sface.sh = shellfacetraverse(subfaces);
1209 }
1210
1211 fprintf(outfile, "\nEnd\n");
1212 fclose(outfile);
1213 }
1214
1215 } // namespace tetgenBR
1216
meshGRegionBoundaryRecovery(GRegion * gr,splitQuadRecovery * sqr)1217 bool meshGRegionBoundaryRecovery(GRegion *gr, splitQuadRecovery *sqr)
1218 {
1219 bool ret = false;
1220 try {
1221 tetgenBR::tetgenmesh *m = new tetgenBR::tetgenmesh();
1222 m->in = new tetgenBR::tetgenio();
1223 m->b = new tetgenBR::tetgenbehavior();
1224 tetgenBR::brdata data = {gr, sqr};
1225 ret = m->reconstructmesh((void *)&data);
1226 delete m->in;
1227 delete m->b;
1228 delete m;
1229 } catch(int err) {
1230 if(err == 1) {
1231 Msg::Error("Out of memory in boundary mesh recovery");
1232 ret = false;
1233 }
1234 else if(err == 3) {
1235 std::map<int, MVertex *> all;
1236 std::vector<GFace *> f = gr->faces();
1237 for(auto it = f.begin(); it != f.end(); ++it) {
1238 GFace *gf = *it;
1239 for(std::size_t i = 0; i < gf->triangles.size(); i++) {
1240 for(int j = 0; j < 3; j++) {
1241 MVertex *v = gf->triangles[i]->getVertex(j);
1242 all[v->getIndex()] = v;
1243 }
1244 }
1245 }
1246 std::vector<GEdge *> const &e = gr->embeddedEdges();
1247 for(auto it = e.begin(); it != e.end(); ++it) {
1248 GEdge *ge = *it;
1249 for(std::size_t i = 0; i < ge->lines.size(); i++) {
1250 for(int j = 0; j < 2; j++) {
1251 MVertex *v = ge->lines[i]->getVertex(j);
1252 all[v->getIndex()] = v;
1253 }
1254 }
1255 }
1256 std::vector<GVertex *> const &v = gr->embeddedVertices();
1257 for(auto it = v.begin(); it != v.end(); ++it) {
1258 GVertex *gv = *it;
1259 for(std::size_t i = 0; i < gv->points.size(); i++) {
1260 MVertex *v = gv->points[i]->getVertex(0);
1261 all[v->getIndex()] = v;
1262 }
1263 }
1264 for(std::size_t i = 0; i < gr->mesh_vertices.size(); i++) {
1265 MVertex *v = gr->mesh_vertices[i];
1266 all[v->getIndex()] = v;
1267 }
1268 std::string what;
1269 bool pnt = true;
1270 switch(tetgenBR::sevent.e_type) {
1271 case 1: what = "segment-segment intersection"; break;
1272 case 2: what = "segment-facet intersection"; break;
1273 case 3: what = "facet-facet intersection"; break;
1274 case 4:
1275 what = "overlapping segments";
1276 pnt = false;
1277 break;
1278 case 5:
1279 what = "segment in facet";
1280 pnt = false;
1281 break;
1282 case 6:
1283 what = "overlapping facets";
1284 pnt = false;
1285 break;
1286 case 7: what = "vertex in segment"; break;
1287 case 8: what = "vertex in facet"; break;
1288 default: what = "unknown"; break;
1289 }
1290 int vtags[2][3] = {
1291 {tetgenBR::sevent.f_vertices1[0], tetgenBR::sevent.f_vertices1[1],
1292 tetgenBR::sevent.f_vertices1[2]},
1293 {tetgenBR::sevent.f_vertices2[0], tetgenBR::sevent.f_vertices2[1],
1294 tetgenBR::sevent.f_vertices2[2]}};
1295 int ftags[2] = {tetgenBR::sevent.f_marker1, tetgenBR::sevent.f_marker2};
1296 int etags[2] = {tetgenBR::sevent.s_marker1, tetgenBR::sevent.s_marker2};
1297 std::ostringstream pb;
1298 std::vector<double> x, y, z, val;
1299 for(int f = 0; f < 2; f++) {
1300 if(ftags[f] > 0) {
1301 GFace *gf = gr->model()->getFaceByTag(ftags[f]);
1302 if(gf) {
1303 gr->model()->addLastMeshEntityError(gf);
1304 pb << " surface " << ftags[f];
1305 }
1306 }
1307 if(etags[f] > 0) {
1308 GEdge *ge = gr->model()->getEdgeByTag(etags[f]);
1309 if(ge) {
1310 gr->model()->addLastMeshEntityError(ge);
1311 pb << " curve " << etags[f];
1312 }
1313 }
1314 for(int i = 0; i < 3; i++) {
1315 MVertex *v = all[vtags[f][i]];
1316 if(v) {
1317 gr->model()->addLastMeshVertexError(v);
1318 x.push_back(v->x());
1319 y.push_back(v->y());
1320 z.push_back(v->z());
1321 val.push_back(f);
1322 }
1323 }
1324 }
1325 if(pnt) {
1326 double px = tetgenBR::sevent.int_point[0];
1327 double py = tetgenBR::sevent.int_point[1];
1328 double pz = tetgenBR::sevent.int_point[2];
1329 pb << ", intersection (" << px << "," << py << "," << pz << ")";
1330 x.push_back(px);
1331 y.push_back(py);
1332 z.push_back(pz);
1333 val.push_back(3.);
1334 }
1335 Msg::Error("Invalid boundary mesh (%s) on%s", what.c_str(),
1336 pb.str().c_str());
1337 #if defined(HAVE_POST)
1338 new PView("Boundary mesh issue", x, y, z, val);
1339 #if defined(HAVE_FLTK)
1340 if(FlGui::available()) FlGui::instance()->updateViews(true, true);
1341 drawContext::global()->draw();
1342 #endif
1343 #endif
1344 ret = false;
1345 }
1346 else {
1347 Msg::Error("Could not recover boundary mesh: error %d", err);
1348 ret = false;
1349 }
1350 }
1351 return ret;
1352 }
1353
1354 #else
1355
meshGRegionBoundaryRecovery(GRegion * gr,splitQuadRecovery * sqr)1356 bool meshGRegionBoundaryRecovery(GRegion *gr, splitQuadRecovery *sqr)
1357 {
1358 return false;
1359 }
1360
1361 #endif
1362