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 // Contributor(s):
7 //   Gaetan Bricteux
8 
9 #include <stdlib.h>
10 #include <sstream>
11 #include "GmshConfig.h"
12 #include "GModel.h"
13 #include "MElement.h"
14 #include "MElementCut.h"
15 #include "gmshLevelset.h"
16 #include "MQuadrangle.h"
17 
18 #if defined(HAVE_DINTEGRATION)
19 #include "Integration3D.h"
20 #endif
21 
22 //---------------------------------------- MPolyhedron
23 //----------------------------
24 
_init()25 void MPolyhedron::_init()
26 {
27   if(_parts.size() == 0) return;
28 
29   for(std::size_t i = 0; i < _parts.size(); i++) {
30     if(_parts[i]->getVolume() * _parts[0]->getVolume() < 0.)
31       _parts[i]->reverse();
32     for(int j = 0; j < 4; j++) {
33       int k;
34       for(k = _faces.size() - 1; k >= 0; k--)
35         if(_parts[i]->getFace(j) == _faces[k]) break;
36       if(k < 0)
37         _faces.push_back(_parts[i]->getFace(j));
38       else
39         _faces.erase(_faces.begin() + k);
40     }
41     if(_parts.size() < 4) { // No interior vertex
42       for(int j = 0; j < 6; j++) {
43         int k;
44         for(k = _edges.size() - 1; k >= 0; k--)
45           if(_parts[i]->getEdge(j) == _edges[k]) break;
46         if(k < 0) _edges.push_back(_parts[i]->getEdge(j));
47       }
48       for(int j = 0; j < 4; j++) {
49         int k;
50         for(k = _vertices.size() - 1; k >= 0; k--)
51           if(_parts[i]->getVertex(j) == _vertices[k]) break;
52         if(k < 0) _vertices.push_back(_parts[i]->getVertex(j));
53       }
54     }
55   }
56 
57   if(_parts.size() >= 4) {
58     for(std::size_t i = 0; i < _faces.size(); i++) {
59       for(int j = 0; j < 3; j++) {
60         int k;
61         for(k = _edges.size() - 1; k >= 0; k--)
62           if(_faces[i].getEdge(j) == _edges[k]) break;
63         if(k < 0) _edges.push_back(_faces[i].getEdge(j));
64       }
65       for(int j = 0; j < 3; j++) {
66         int k;
67         for(k = _vertices.size() - 1; k >= 0; k--)
68           if(_faces[i].getVertex(j) == _vertices[k]) break;
69         if(k < 0) _vertices.push_back(_faces[i].getVertex(j));
70       }
71     }
72   }
73 
74   // innerVertices
75   for(std::size_t i = 0; i < _parts.size(); i++) {
76     for(int j = 0; j < 4; j++) {
77       if(std::find(_vertices.begin(), _vertices.end(),
78                    _parts[i]->getVertex(j)) == _vertices.end())
79         _innerVertices.push_back(_parts[i]->getVertex(j));
80     }
81   }
82 }
83 
isInside(double u,double v,double w) const84 bool MPolyhedron::isInside(double u, double v, double w) const
85 {
86   if(!_orig) return false;
87   double uvw[3] = {u, v, w};
88   for(std::size_t i = 0; i < _parts.size(); i++) {
89     double verts[4][3];
90     for(int j = 0; j < 4; j++) {
91       MVertex *vij = _parts[i]->getVertex(j);
92       double v_xyz[3] = {vij->x(), vij->y(), vij->z()};
93       _orig->xyz2uvw(v_xyz, verts[j]);
94     }
95     MVertex v0(verts[0][0], verts[0][1], verts[0][2]);
96     MVertex v1(verts[1][0], verts[1][1], verts[1][2]);
97     MVertex v2(verts[2][0], verts[2][1], verts[2][2]);
98     MVertex v3(verts[3][0], verts[3][1], verts[3][2]);
99     MTetrahedron t(&v0, &v1, &v2, &v3);
100     double ksi[3];
101     t.xyz2uvw(uvw, ksi);
102     if(t.isInside(ksi[0], ksi[1], ksi[2])) return true;
103   }
104   return false;
105 }
106 
getIntegrationPoints(int pOrder,int * npts,IntPt ** pts)107 void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
108 {
109   *npts = 0;
110   if(_intpt) delete[] _intpt;
111   if(!_orig) return;
112   double jac[3][3];
113   _intpt = new IntPt[getNGQTetPts(pOrder) * _parts.size()];
114   for(std::size_t i = 0; i < _parts.size(); i++) {
115     int nptsi;
116     IntPt *ptsi;
117     _parts[i]->getIntegrationPoints(pOrder, &nptsi, &ptsi);
118     double uvw[4][3];
119     for(int j = 0; j < 4; j++) {
120       double xyz[3] = {_parts[i]->getVertex(j)->x(),
121                        _parts[i]->getVertex(j)->y(),
122                        _parts[i]->getVertex(j)->z()};
123       _orig->xyz2uvw(xyz, uvw[j]);
124     }
125     MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
126     MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
127     MVertex v2(uvw[2][0], uvw[2][1], uvw[2][2]);
128     MVertex v3(uvw[3][0], uvw[3][1], uvw[3][2]);
129     MTetrahedron tt(&v0, &v1, &v2, &v3);
130 
131     for(int ip = 0; ip < nptsi; ip++) {
132       const double u = ptsi[ip].pt[0];
133       const double v = ptsi[ip].pt[1];
134       const double w = ptsi[ip].pt[2];
135       SPoint3 p;
136       tt.pnt(u, v, w, p);
137       _intpt[*npts + ip].pt[0] = p.x();
138       _intpt[*npts + ip].pt[1] = p.y();
139       _intpt[*npts + ip].pt[2] = p.z();
140       double partJac = _parts[i]->getJacobian(p.x(), p.y(), p.z(), jac);
141       double Jac = getJacobian(p.x(), p.y(), p.z(), jac);
142       _intpt[*npts + ip].weight = ptsi[ip].weight * partJac / Jac;
143     }
144     *npts += nptsi;
145   }
146   *pts = _intpt;
147 }
148 
149 //------------------------------------------- MPolygon
150 //------------------------------
151 
_initVertices()152 void MPolygon::_initVertices()
153 {
154   if(_parts.size() == 0) return;
155 
156   // reorient the parts
157   SVector3 n;
158   if(_orig)
159     n = _orig->getFace(0).normal();
160   else
161     n = _parts[0]->getFace(0).normal();
162   for(std::size_t i = 0; i < _parts.size(); i++) {
163     SVector3 ni = _parts[i]->getFace(0).normal();
164     if(dot(n, ni) < 0.) _parts[i]->reverse();
165   }
166   // select only outer edges in vector edg
167   std::vector<MEdge> edg;
168   std::vector<MEdge> multiEdges;
169   edg.reserve(_parts[0]->getNumEdges());
170   for(int j = 0; j < _parts[0]->getNumEdges(); j++)
171     edg.push_back(_parts[0]->getEdge(j));
172   for(std::size_t i = 1; i < _parts.size(); i++) {
173     for(int j = 0; j < _parts[i]->getNumEdges(); j++) {
174       bool found = false;
175       MEdge ed = _parts[i]->getEdge(j);
176       int k;
177       for(k = edg.size() - 1; k >= 0; k--) {
178         if(ed == edg[k]) {
179           edg.erase(edg.begin() + k);
180           found = true;
181           break;
182         }
183       }
184       if(!found) {
185         for(k = 0; k < (int)multiEdges.size(); k++)
186           if(multiEdges[k].isInside(ed.getVertex(0)) &&
187              multiEdges[k].isInside(ed.getVertex(1))) {
188             found = true;
189             break;
190           }
191       }
192       if(!found) {
193         for(k = edg.size() - 1; k >= 0; k--) {
194           if(edg[k].isInside(ed.getVertex(0)) &&
195              edg[k].isInside(ed.getVertex(1))) {
196             multiEdges.push_back(edg[k]);
197             edg.erase(edg.begin() + k);
198             found = true;
199             break;
200           }
201         }
202       }
203       if(!found) {
204         for(k = edg.size() - 1; k >= 0; k--) {
205           if(ed.isInside(edg[k].getVertex(0)) &&
206              ed.isInside(edg[k].getVertex(1))) {
207             edg.erase(edg.begin() + k);
208             int nbME = multiEdges.size();
209             if(nbME == 0 || multiEdges[nbME - 1] != ed)
210               multiEdges.push_back(ed);
211             found = true;
212           }
213         }
214       }
215       if(!found) edg.push_back(ed);
216     }
217   }
218 
219   // organize edges to get vertices in rotating order
220   _edges.push_back(edg[0]);
221   edg.erase(edg.begin());
222   while(edg.size()) {
223     MVertex *v = _edges[_edges.size() - 1].getVertex(1);
224     for(std::size_t i = 0; i < edg.size(); i++) {
225       if(edg[i].getVertex(0) == v) {
226         _edges.push_back(edg[i]);
227         edg.erase(edg.begin() + i);
228         break;
229       }
230       if(edg[i].getVertex(1) == v) {
231         _edges.push_back(MEdge(edg[i].getVertex(1), edg[i].getVertex(0)));
232         edg.erase(edg.begin() + i);
233         break;
234       }
235       if(i == edg.size() - 1) {
236         _edges.push_back(edg[0]);
237         edg.erase(edg.begin());
238         break;
239       }
240     }
241   }
242 
243   for(std::size_t i = 0; i < _edges.size(); i++) {
244     _vertices.push_back(_edges[i].getVertex(0));
245   }
246 
247   // innerVertices
248   for(std::size_t i = 0; i < _parts.size(); i++) {
249     for(int j = 0; j < 3; j++) {
250       if(std::find(_vertices.begin(), _vertices.end(),
251                    _parts[i]->getVertex(j)) == _vertices.end())
252         _innerVertices.push_back(_parts[i]->getVertex(j));
253     }
254   }
255 }
256 
isInside(double u,double v,double w) const257 bool MPolygon::isInside(double u, double v, double w) const
258 {
259   if(!getParent()) return false;
260   double uvw[3] = {u, v, w};
261   for(std::size_t i = 0; i < _parts.size(); i++) {
262     double v_uvw[3][3];
263     for(int j = 0; j < 3; j++) {
264       MVertex *vij = _parts[i]->getVertex(j);
265       double v_xyz[3] = {vij->x(), vij->y(), vij->z()};
266       getParent()->xyz2uvw(v_xyz, v_uvw[j]);
267     }
268     MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
269     MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
270     MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
271     MTriangle t(&v0, &v1, &v2);
272     double ksi[3];
273     t.xyz2uvw(uvw, ksi);
274     if(t.isInside(ksi[0], ksi[1], ksi[2])) return true;
275   }
276   return false;
277 }
278 
getIntegrationPoints(int pOrder,int * npts,IntPt ** pts)279 void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
280 {
281   *npts = 0;
282   if(_intpt) delete[] _intpt;
283   if(!getParent()) return;
284   double jac[3][3];
285   _intpt = new IntPt[getNGQTPts(pOrder) * _parts.size()];
286   for(std::size_t i = 0; i < _parts.size(); i++) {
287     int nptsi;
288     IntPt *ptsi;
289     _parts[i]->getIntegrationPoints(pOrder, &nptsi, &ptsi);
290     double uvw[3][3];
291     for(int j = 0; j < 3; j++) {
292       double xyz[3] = {_parts[i]->getVertex(j)->x(),
293                        _parts[i]->getVertex(j)->y(),
294                        _parts[i]->getVertex(j)->z()};
295       getParent()->xyz2uvw(xyz, uvw[j]);
296     }
297     MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
298     MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
299     MVertex v2(uvw[2][0], uvw[2][1], uvw[2][2]);
300     MTriangle tt(&v0, &v1, &v2);
301     for(int ip = 0; ip < nptsi; ip++) {
302       const double u = ptsi[ip].pt[0];
303       const double v = ptsi[ip].pt[1];
304       const double w = ptsi[ip].pt[2];
305       SPoint3 p;
306       tt.pnt(u, v, w, p);
307       _intpt[*npts + ip].pt[0] = p.x();
308       _intpt[*npts + ip].pt[1] = p.y();
309       _intpt[*npts + ip].pt[2] = p.z();
310       double partJac = _parts[i]->getJacobian(p.x(), p.y(), p.z(), jac);
311       double Jac = getJacobian(p.x(), p.y(), p.z(), jac);
312       _intpt[*npts + ip].weight = ptsi[ip].weight * partJac / Jac;
313     }
314     *npts += nptsi;
315   }
316   *pts = _intpt;
317 }
318 
319 //----------------------------------- MLineChild ------------------------------
320 
isInside(double u,double v,double w) const321 bool MLineChild::isInside(double u, double v, double w) const
322 {
323   if(!_orig) return false;
324   double uvw[3] = {u, v, w};
325   double v_uvw[2][3];
326   for(int i = 0; i < 2; i++) {
327     const MVertex *vi = getVertex(i);
328     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
329     _orig->xyz2uvw(v_xyz, v_uvw[i]);
330   }
331   MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
332   MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
333   MLine l(&v0, &v1);
334   double ksi[3];
335   l.xyz2uvw(uvw, ksi);
336   if(l.isInside(ksi[0], ksi[1], ksi[2])) return true;
337   return false;
338 }
339 
getIntegrationPoints(int pOrder,int * npts,IntPt ** pts)340 void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
341 {
342   *npts = 0;
343   if(_intpt) delete[] _intpt;
344   if(!_orig) return;
345   _intpt = new IntPt[getNGQLPts(pOrder)];
346   int nptsi;
347   IntPt *ptsi;
348   double v_uvw[2][3];
349   for(int i = 0; i < 2; i++) {
350     MVertex *vi = getVertex(i);
351     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
352     _orig->xyz2uvw(v_xyz, v_uvw[i]);
353   }
354   MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
355   MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
356   MLine l(&v0, &v1);
357   l.getIntegrationPoints(pOrder, &nptsi, &ptsi);
358   for(int ip = 0; ip < nptsi; ip++) {
359     const double u = ptsi[ip].pt[0];
360     const double v = ptsi[ip].pt[1];
361     const double w = ptsi[ip].pt[2];
362     SPoint3 p;
363     l.pnt(u, v, w, p);
364     _intpt[*npts + ip].pt[0] = p.x();
365     _intpt[*npts + ip].pt[1] = p.y();
366     _intpt[*npts + ip].pt[2] = p.z();
367     _intpt[*npts + ip].weight = ptsi[ip].weight;
368   }
369   *npts = nptsi;
370   *pts = _intpt;
371 }
372 
373 //----------------------------------- MTriangleBorder
374 //------------------------------
375 
isInside(double u,double v,double w) const376 bool MTriangleBorder::isInside(double u, double v, double w) const
377 {
378   if(!getParent()) return false;
379   double uvw[3] = {u, v, w};
380   double v_uvw[3][3];
381   for(int i = 0; i < 3; i++) {
382     const MVertex *vi = getVertex(i);
383     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
384     getParent()->xyz2uvw(v_xyz, v_uvw[i]);
385   }
386   MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
387   MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
388   MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
389   MTriangle t(&v0, &v1, &v2);
390   double ksi[3];
391   t.xyz2uvw(uvw, ksi);
392   if(t.isInside(ksi[0], ksi[1], ksi[2])) return true;
393   return false;
394 }
395 
getIntegrationPoints(int pOrder,int * npts,IntPt ** pts)396 void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
397 {
398   *npts = 0;
399   if(_intpt) delete[] _intpt;
400   if(!getParent()) return;
401   _intpt = new IntPt[getNGQTPts(pOrder)];
402   int nptsi;
403   IntPt *ptsi;
404 
405   double uvw[3][3];
406   for(int j = 0; j < 3; j++) {
407     double xyz[3] = {_v[j]->x(), _v[j]->y(), _v[j]->z()};
408     getParent()->xyz2uvw(xyz, uvw[j]);
409   }
410   MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
411   MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
412   MVertex v2(uvw[2][0], uvw[2][1], uvw[2][2]);
413   MTriangle tt(&v0, &v1, &v2);
414   tt.getIntegrationPoints(pOrder, &nptsi, &ptsi);
415   double jac[3][3];
416   for(int ip = 0; ip < nptsi; ip++) {
417     const double u = ptsi[ip].pt[0];
418     const double v = ptsi[ip].pt[1];
419     const double w = ptsi[ip].pt[2];
420     tt.getJacobian(u, v, w, jac);
421     SPoint3 p;
422     tt.pnt(u, v, w, p);
423     _intpt[ip].pt[0] = p.x();
424     _intpt[ip].pt[1] = p.y();
425     _intpt[ip].pt[2] = p.z();
426     _intpt[ip].weight = ptsi[ip].weight;
427   }
428   *npts = nptsi;
429   *pts = _intpt;
430 }
431 
432 //-------------------------------------- MLineBorder
433 //------------------------------
434 
isInside(double u,double v,double w) const435 bool MLineBorder::isInside(double u, double v, double w) const
436 {
437   if(!getParent()) return false;
438   double uvw[3] = {u, v, w};
439   double v_uvw[2][3];
440   for(int i = 0; i < 2; i++) {
441     const MVertex *vi = getVertex(i);
442     double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
443     getParent()->xyz2uvw(v_xyz, v_uvw[i]);
444   }
445   MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
446   MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
447   MLine l(&v0, &v1);
448   double ksi[3];
449   l.xyz2uvw(uvw, ksi);
450   if(l.isInside(ksi[0], ksi[1], ksi[2])) return true;
451   return false;
452 }
453 
getIntegrationPoints(int pOrder,int * npts,IntPt ** pts)454 void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
455 {
456   *npts = 0;
457   if(_intpt) delete[] _intpt;
458   if(!getParent()) return;
459   _intpt = new IntPt[getNGQLPts(pOrder)];
460   int nptsi;
461   IntPt *ptsi;
462 
463   double uvw[2][3];
464   for(int j = 0; j < 2; j++) {
465     double xyz[3] = {_v[j]->x(), _v[j]->y(), _v[j]->z()};
466     getParent()->xyz2uvw(xyz, uvw[j]);
467   }
468   MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
469   MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
470   MLine ll(&v0, &v1);
471   ll.getIntegrationPoints(pOrder, &nptsi, &ptsi);
472   for(int ip = 0; ip < nptsi; ip++) {
473     const double u = ptsi[ip].pt[0];
474     const double v = ptsi[ip].pt[1];
475     const double w = ptsi[ip].pt[2];
476     SPoint3 p;
477     ll.pnt(u, v, w, p);
478     _intpt[ip].pt[0] = p.x();
479     _intpt[ip].pt[1] = p.y();
480     _intpt[ip].pt[2] = p.z();
481     _intpt[ip].weight = ptsi[ip].weight;
482   }
483   *npts = nptsi;
484   *pts = _intpt;
485 }
486 
487 //---------------------------------------- CutMesh
488 //----------------------------
489 
490 #if defined(HAVE_DINTEGRATION)
491 
492 static void
assignPhysicals(GModel * GM,std::vector<int> & gePhysicals,int reg,int dim,std::map<int,std::map<int,std::string>> physicals[4],std::map<int,int> & newPhysTags,int lsTag)493 assignPhysicals(GModel *GM, std::vector<int> &gePhysicals, int reg, int dim,
494                 std::map<int, std::map<int, std::string> > physicals[4],
495                 std::map<int, int> &newPhysTags, int lsTag)
496 {
497   for(std::size_t i = 0; i < gePhysicals.size(); i++) {
498     int phys = gePhysicals[i];
499 
500     if(lsTag > 0 && newPhysTags.count(phys)) {
501       int phys2 = newPhysTags[phys];
502       if(phys2 &&
503          (!physicals[dim].count(reg) || !physicals[dim][reg].count(phys2))) {
504         std::string name = GM->getPhysicalName(dim, phys);
505         if(name != "" && newPhysTags.count(-phys)) {
506           auto it = physicals[dim].begin();
507           for(; it != physicals[dim].end(); it++) {
508             auto it2 = it->second.begin();
509             for(; it2 != it->second.end(); it2++)
510               if(it2->second == name)
511                 physicals[dim][it->first][it2->first] = name + "_out";
512           }
513           name += "_in";
514         }
515         physicals[dim][reg][phys2] = name;
516       }
517     }
518     else if(lsTag < 0 && newPhysTags.count(-phys)) {
519       int phys2 = newPhysTags[-phys];
520       if(phys2 &&
521          (!physicals[dim].count(reg) || !physicals[dim][reg].count(phys2))) {
522         std::string name = GM->getPhysicalName(dim, phys);
523         if(name != "" && newPhysTags.count(phys)) {
524           auto it = physicals[dim].begin();
525           for(; it != physicals[dim].end(); it++) {
526             auto it2 = it->second.begin();
527             for(; it2 != it->second.end(); it2++)
528               if(it2->second == name)
529                 physicals[dim][it->first][it2->first] = name + "_in";
530           }
531           name += "_out";
532         }
533         physicals[dim][reg][phys2] = name;
534       }
535     }
536   }
537 }
538 
539 static void
assignLsPhysical(GModel * GM,int reg,int dim,std::map<int,std::map<int,std::string>> physicals[4],int physTag,int lsTag)540 assignLsPhysical(GModel *GM, int reg, int dim,
541                  std::map<int, std::map<int, std::string> > physicals[4],
542                  int physTag, int lsTag)
543 {
544   if(!physicals[dim][reg].count(physTag)) {
545     std::stringstream strs;
546     strs << lsTag;
547     std::string sdim = (dim == 2) ? "S" : "L";
548     physicals[dim][reg][physTag] = "levelset_" + sdim + strs.str();
549     if(physTag != lsTag)
550       Msg::Info("Levelset %d -> physical %d", lsTag, physTag);
551   }
552 }
553 
getElementaryTag(int lsTag,int elementary,std::map<int,int> & newElemTags)554 static int getElementaryTag(int lsTag, int elementary,
555                             std::map<int, int> &newElemTags)
556 {
557   if(lsTag < 0) {
558     if(newElemTags.count(elementary))
559       return newElemTags[elementary];
560     else {
561       int reg = ++newElemTags[0];
562       newElemTags[elementary] = reg;
563       return reg;
564     }
565   }
566   return elementary;
567 }
getPhysicalTag(int lsTag,const std::vector<int> & physicals,std::map<int,int> & newPhysTags)568 static void getPhysicalTag(int lsTag, const std::vector<int> &physicals,
569                            std::map<int, int> &newPhysTags)
570 {
571   for(std::size_t i = 0; i < physicals.size(); i++) {
572     int phys = physicals[i];
573     if(lsTag < 0) {
574       if(!newPhysTags.count(-phys)) newPhysTags[-phys] = ++newPhysTags[0];
575       phys = newPhysTags[-phys];
576     }
577     else if(!newPhysTags.count(phys))
578       newPhysTags[phys] = phys;
579   }
580 }
581 
getBorderTag(int lsTag,int count,int & maxTag,std::map<int,int> & borderTags)582 static int getBorderTag(int lsTag, int count, int &maxTag,
583                         std::map<int, int> &borderTags)
584 {
585   if(borderTags.count(lsTag)) return borderTags[lsTag];
586   if(count || lsTag < 0) {
587     int tag = ++maxTag;
588     borderTags[lsTag] = tag;
589     return tag;
590   }
591   maxTag = std::max(maxTag, lsTag);
592   borderTags[lsTag] = lsTag;
593   return lsTag;
594 }
595 
elementSplitMesh(MElement * e,std::vector<gLevelset * > & RPN,fullMatrix<double> & verticesLs,GEntity * ge,GModel * GM,std::size_t & numEle,std::map<int,MVertex * > & vertexMap,std::map<MElement *,MElement * > & newParents,std::map<MElement *,MElement * > & newDomains,std::map<int,std::vector<MElement * >> elements[10],std::map<int,std::map<int,std::string>> physicals[4],std::map<int,int> newElemTags[4],std::map<int,int> newPhysTags[4],std::map<int,int> borderElemTags[2],std::map<int,int> borderPhysTags[2])596 static void elementSplitMesh(
597   MElement *e, std::vector<gLevelset *> &RPN, fullMatrix<double> &verticesLs,
598   GEntity *ge, GModel *GM, std::size_t &numEle,
599   std::map<int, MVertex *> &vertexMap,
600   std::map<MElement *, MElement *> &newParents,
601   std::map<MElement *, MElement *> &newDomains,
602   std::map<int, std::vector<MElement *> > elements[10],
603   std::map<int, std::map<int, std::string> > physicals[4],
604   std::map<int, int> newElemTags[4], std::map<int, int> newPhysTags[4],
605   std::map<int, int> borderElemTags[2], std::map<int, int> borderPhysTags[2])
606 {
607   int elementary = ge->tag();
608   int eType = e->getType();
609   std::vector<int> gePhysicals = ge->physicals;
610   int iLs = (verticesLs.size1() == 1) ? 0 : verticesLs.size1() - 1;
611   int gLsTag = RPN[RPN.size() - 1]->getTag();
612 
613   MElement *copy = e->copy(vertexMap, newParents, newDomains);
614   double eps = 1.e-10;
615   bool splitElem = false;
616 
617   // split acording to center of gravity
618   // double lsMean = 0.0;
619   // int lsTag = 1;
620   // for(int k = 0; k < e->getNumVertices(); k++)
621   //   lsMean += verticesLs(iLs, e->getVertex(k)->getIndex()) - eps;
622   // if (lsMean > 0) lsTag = -1;
623 
624   // EMI : better for embedded dirichlet with smoothed properties
625   // split according to values of vertices (keep +)
626   int lsTag = (verticesLs(iLs, e->getVertex(0)->getIndex()) > eps) ? -1 : 1;
627   int ils = 0;
628   for(std::size_t k = 1; k < e->getNumVertices(); k++) {
629     int lsTag2 = (verticesLs(iLs, e->getVertex(k)->getIndex()) > eps) ? -1 : 1;
630     if(lsTag * lsTag2 < 0) {
631       lsTag = -1;
632       splitElem = true;
633       if(RPN.size() > 1) {
634         for(; ils < verticesLs.size1() - 1; ils++) {
635           int lsT1 =
636             (verticesLs(ils, e->getVertex(0)->getIndex()) > eps) ? -1 : 1;
637           int lsT2 =
638             (verticesLs(ils, e->getVertex(k)->getIndex()) > eps) ? -1 : 1;
639           if(lsT1 * lsT2 < 0) break;
640         }
641         for(int i = 0; i <= ils; i++)
642           if(!RPN[i]->isPrimitive()) ils++;
643         gLsTag = RPN[ils]->getTag();
644       }
645       break;
646     }
647   }
648 
649   // invert tag
650   // lsTag = 1; //negative ls
651   // for(int k = 0; k < e->getNumVertices(); k++){
652   //  double val = verticesLs(iLs, e->getVertex(k)->getIndex());
653   //  if (val > 0.0) { lsTag = -1; break; }
654   //}
655 
656   switch(eType) {
657   case TYPE_TET:
658   case TYPE_HEX:
659   case TYPE_PYR:
660   case TYPE_PRI:
661   case TYPE_POLYH: {
662     int reg = getElementaryTag(lsTag, elementary, newElemTags[3]);
663     getPhysicalTag(lsTag, gePhysicals, newPhysTags[3]);
664     if(eType == TYPE_TET)
665       elements[4][reg].push_back(copy);
666     else if(eType == TYPE_HEX)
667       elements[5][reg].push_back(copy);
668     else if(eType == TYPE_PRI)
669       elements[6][reg].push_back(copy);
670     else if(eType == TYPE_PYR)
671       elements[7][reg].push_back(copy);
672     else if(eType == TYPE_POLYH)
673       elements[9][reg].push_back(copy);
674     assignPhysicals(GM, gePhysicals, reg, 3, physicals, newPhysTags[3], lsTag);
675   } break;
676   case TYPE_TRI:
677   case TYPE_QUA:
678   case TYPE_POLYG: {
679     int reg = getElementaryTag(lsTag, elementary, newElemTags[2]);
680     getPhysicalTag(lsTag, gePhysicals, newPhysTags[2]);
681     if(eType == TYPE_TRI)
682       elements[2][reg].push_back(copy);
683     else if(eType == TYPE_QUA)
684       elements[3][reg].push_back(copy);
685     else if(eType == TYPE_POLYG)
686       elements[8][reg].push_back(copy);
687     assignPhysicals(GM, gePhysicals, reg, 2, physicals, newPhysTags[2], lsTag);
688   } break;
689   case TYPE_LIN: {
690     int reg = getElementaryTag(lsTag, elementary, newElemTags[1]);
691     getPhysicalTag(lsTag, gePhysicals, newPhysTags[1]);
692     elements[1][reg].push_back(copy);
693     assignPhysicals(GM, gePhysicals, reg, 1, physicals, newPhysTags[1], lsTag);
694   } break;
695   case TYPE_PNT: {
696     int reg = getElementaryTag(lsTag, elementary, newElemTags[0]);
697     getPhysicalTag(lsTag, gePhysicals, newPhysTags[0]);
698     elements[0][reg].push_back(copy);
699     assignPhysicals(GM, gePhysicals, reg, 0, physicals, newPhysTags[0], lsTag);
700   } break;
701   default: Msg::Error("This type of element cannot be split."); return;
702   }
703 
704   // create level set interface (pt in 1D, line in 2D or face in 3D)
705   if(splitElem && e->getDim() == 2) {
706     for(int k = 0; k < e->getNumEdges(); k++) {
707       MEdge me = e->getEdge(k);
708       MVertex *v0 = me.getVertex(0);
709       MVertex *v1 = me.getVertex(1);
710       MVertex *v0N = vertexMap[v0->getNum()];
711       MVertex *v1N = vertexMap[v1->getNum()];
712       double val0 = (verticesLs(iLs, v0->getIndex()) > eps) ? 1 : -1;
713       double val1 = (verticesLs(iLs, v1->getIndex()) > eps) ? 1 : -1;
714       if(val0 * val1 > 0.0 && val0 < 0.0) {
715         getPhysicalTag(-1, gePhysicals, newPhysTags[1]);
716         int c = elements[1].count(gLsTag);
717         int reg = getBorderTag(gLsTag, c, newElemTags[1][0], borderElemTags[0]);
718         int physTag =
719           (!gePhysicals.size()) ?
720             0 :
721             getBorderTag(gLsTag, c, newPhysTags[1][0], borderPhysTags[0]);
722         int i;
723         for(i = elements[1][reg].size() - 1; i >= 0; i--) {
724           MElement *el = elements[1][reg][i];
725           if((el->getVertex(0) == v0N && el->getVertex(1) == v1N) ||
726              (el->getVertex(0) == v1N && el->getVertex(1) == v0N))
727             break;
728         }
729         if(i < 0) {
730           MLine *lin = new MLine(v0N, v1N);
731           elements[1][reg].push_back(lin);
732           if(physTag) assignLsPhysical(GM, reg, 1, physicals, physTag, gLsTag);
733         }
734         else {
735           MElement *el = elements[1][reg][i];
736           elements[1][reg].erase(elements[1][reg].begin() + i);
737           delete el;
738         }
739       }
740     }
741   }
742   else if(splitElem && e->getDim() == 3) {
743     for(int k = 0; k < e->getNumFaces(); k++) {
744       MFace mf = e->getFace(k);
745       bool sameSign = true;
746       double val0 =
747         (verticesLs(iLs, mf.getVertex(0)->getIndex()) > eps) ? 1 : -1;
748       for(std::size_t j = 1; j < mf.getNumVertices(); j++) {
749         double valj =
750           (verticesLs(iLs, mf.getVertex(j)->getIndex()) > eps) ? 1 : -1;
751         if(val0 * valj < 0.0) {
752           sameSign = false;
753           break;
754         }
755       }
756       if(sameSign && val0 < 0.0) {
757         getPhysicalTag(-1, gePhysicals, newPhysTags[2]);
758         int c = elements[2].count(gLsTag) + elements[3].count(gLsTag) +
759                 elements[8].count(gLsTag);
760         int reg = getBorderTag(gLsTag, c, newElemTags[2][0], borderElemTags[1]);
761         int physTag =
762           (!gePhysicals.size()) ?
763             0 :
764             getBorderTag(gLsTag, c, newPhysTags[2][0], borderPhysTags[1]);
765         if(mf.getNumVertices() == 3) {
766           MFace f1 = MFace(vertexMap[mf.getVertex(0)->getNum()],
767                            vertexMap[mf.getVertex(1)->getNum()],
768                            vertexMap[mf.getVertex(2)->getNum()]);
769           int i = elements[2][reg].size() - 1;
770           for(; i >= 0; i--) {
771             MElement *el = elements[2][reg][i];
772             MFace f2 =
773               MFace(el->getVertex(0), el->getVertex(1), el->getVertex(2));
774             if(f1 == f2) break;
775           }
776           if(i < 0) {
777             MTriangle *tri =
778               new MTriangle(f1.getVertex(0), f1.getVertex(1), f1.getVertex(2));
779             elements[2][reg].push_back(tri);
780           }
781           else {
782             MElement *el = elements[2][reg][i];
783             elements[2][reg].erase(elements[2][reg].begin() + i);
784             delete el;
785           }
786         }
787         else if(mf.getNumVertices() == 4) {
788           MFace f1 = MFace(vertexMap[mf.getVertex(0)->getNum()],
789                            vertexMap[mf.getVertex(1)->getNum()],
790                            vertexMap[mf.getVertex(2)->getNum()],
791                            vertexMap[mf.getVertex(3)->getNum()]);
792           int i;
793           for(i = elements[3][reg].size() - 1; i >= 0; i--) {
794             MElement *el = elements[3][reg][i];
795             MFace f2 = MFace(el->getVertex(0), el->getVertex(1),
796                              el->getVertex(2), el->getVertex(3));
797             if(f1 == f2) break;
798           }
799           if(i < 0) {
800             MQuadrangle *quad =
801               new MQuadrangle(f1.getVertex(0), f1.getVertex(1), f1.getVertex(2),
802                               f1.getVertex(2));
803             elements[3][reg].push_back(quad);
804           }
805           else {
806             MElement *el = elements[3][reg][i];
807             elements[3][reg].erase(elements[3][reg].begin() + i);
808             delete el;
809           }
810         }
811         if(physTag) assignLsPhysical(GM, reg, 2, physicals, physTag, gLsTag);
812       }
813     }
814   }
815 }
816 
equalV(MVertex * v,const DI_Point * p)817 static bool equalV(MVertex *v, const DI_Point *p)
818 {
819   return (fabs(v->x() - p->x()) < 1.e-15 && fabs(v->y() - p->y()) < 1.e-15 &&
820           fabs(v->z() - p->z()) < 1.e-15);
821 }
822 
getElementVertexNum(DI_Point * p,MElement * e)823 static int getElementVertexNum(DI_Point *p, MElement *e)
824 {
825   for(std::size_t i = 0; i < e->getNumVertices(); i++)
826     if(equalV(e->getVertex(i), p)) return e->getVertex(i)->getNum();
827   return -1;
828 }
829 
830 typedef std::set<MVertex *, MVertexPtrLessThanLexicographic>
831   newVerticesContainer;
832 
elementCutMesh(MElement * e,std::vector<gLevelset * > & RPN,fullMatrix<double> & verticesLs,GEntity * ge,GModel * GM,std::size_t & numEle,std::map<int,MVertex * > & vertexMap,newVerticesContainer & newVertices,std::map<MElement *,MElement * > & newParents,std::map<MElement *,MElement * > & newDomains,std::multimap<MElement *,MElement * > borders[2],std::map<int,std::vector<MElement * >> elements[10],std::map<int,std::map<int,std::string>> physicals[4],std::map<int,int> newElemTags[4],std::map<int,int> newPhysTags[4],std::map<int,int> borderElemTags[2],std::map<int,int> borderPhysTags[2],DI_Point::Container & cp,std::vector<DI_Line * > & lines,std::vector<DI_Triangle * > & triangles,std::vector<DI_Quad * > & quads,std::vector<DI_Tetra * > & tetras,std::vector<DI_Hexa * > & hexas)833 static void elementCutMesh(
834   MElement *e, std::vector<gLevelset *> &RPN, fullMatrix<double> &verticesLs,
835   GEntity *ge, GModel *GM, std::size_t &numEle,
836   std::map<int, MVertex *> &vertexMap, newVerticesContainer &newVertices,
837   std::map<MElement *, MElement *> &newParents,
838   std::map<MElement *, MElement *> &newDomains,
839   std::multimap<MElement *, MElement *> borders[2],
840   std::map<int, std::vector<MElement *> > elements[10],
841   std::map<int, std::map<int, std::string> > physicals[4],
842   std::map<int, int> newElemTags[4], std::map<int, int> newPhysTags[4],
843   std::map<int, int> borderElemTags[2], std::map<int, int> borderPhysTags[2],
844   DI_Point::Container &cp, std::vector<DI_Line *> &lines,
845   std::vector<DI_Triangle *> &triangles, std::vector<DI_Quad *> &quads,
846   std::vector<DI_Tetra *> &tetras, std::vector<DI_Hexa *> &hexas)
847 {
848   int elementary = ge->tag();
849   int eType = e->getType();
850   int recur = e->getPolynomialOrder() - 1;
851   int ePart = e->getPartition();
852   MElement *eParent = e->getParent();
853   std::vector<int> gePhysicals = ge->physicals;
854 
855   int integOrder = 1;
856   std::vector<DI_IntegrationPoint *> ipV;
857   std::vector<DI_IntegrationPoint *> ipS;
858   bool isCut = false;
859   std::size_t nbL = lines.size();
860   std::size_t nbTr = triangles.size();
861   std::size_t nbQ = quads.size();
862   std::size_t nbTe = tetras.size();
863   std::size_t nbH = hexas.size();
864 
865   MElement *copy = e->copy(vertexMap, newParents, newDomains);
866   MElement *parent = eParent ? copy->getParent() : copy;
867 
868   double **nodeLs = new double *[e->getNumVertices()];
869 
870   switch(eType) {
871   case TYPE_TET:
872   case TYPE_HEX:
873   case TYPE_PYR:
874   case TYPE_PRI:
875   case TYPE_POLYH: {
876     if(eType == TYPE_TET) {
877       DI_Tetra T(
878         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
879         e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
880         e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
881         e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
882       T.setPolynomialOrder(recur + 1);
883       for(std::size_t i = 0; i < e->getNumVertices(); i++)
884         nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
885       isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
886                     tetras, quads, triangles, recur, nodeLs);
887     }
888     else if(eType == TYPE_HEX) {
889       DI_Hexa H(
890         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
891         e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
892         e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
893         e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z(),
894         e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z(),
895         e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z(),
896         e->getVertex(6)->x(), e->getVertex(6)->y(), e->getVertex(6)->z(),
897         e->getVertex(7)->x(), e->getVertex(7)->y(), e->getVertex(7)->z());
898       H.setPolynomialOrder(recur + 1);
899       for(std::size_t i = 0; i < e->getNumVertices(); i++)
900         nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
901       isCut =
902         H.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
903               hexas, tetras, quads, triangles, lines, recur, nodeLs);
904     }
905     else if(eType == TYPE_PRI) {
906       DI_Tetra T1(
907         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
908         e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
909         e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
910         e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
911       // T1.setPolynomialOrder(recur+1);
912       for(int i = 0; i < 3; i++)
913         nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
914       nodeLs[3] = &verticesLs(0, e->getVertex(5)->getIndex());
915       bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
916                         tetras, quads, triangles, recur, nodeLs);
917       DI_Tetra T2(
918         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
919         e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z(),
920         e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
921         e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
922       // T2.setPolynomialOrder(recur+1);
923       nodeLs[0] = &verticesLs(0, e->getVertex(0)->getIndex());
924       nodeLs[1] = &verticesLs(0, e->getVertex(4)->getIndex());
925       nodeLs[2] = &verticesLs(0, e->getVertex(1)->getIndex());
926       nodeLs[3] = &verticesLs(0, e->getVertex(5)->getIndex());
927       bool iC2 = T2.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
928                         tetras, quads, triangles, recur, nodeLs);
929       DI_Tetra T3(
930         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
931         e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z(),
932         e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z(),
933         e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
934       // T3.setPolynomialOrder(recur+1);
935       for(int i = 1; i < 4; i++)
936         nodeLs[i] = &verticesLs(0, e->getVertex(i + 2)->getIndex());
937       bool iC3 = T3.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
938                         tetras, quads, triangles, recur, nodeLs);
939       isCut = iC1 || iC2 || iC3;
940     }
941     else if(eType == TYPE_PYR) {
942       DI_Tetra T1(
943         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
944         e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
945         e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
946         e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
947       // T1.setPolynomialOrder(recur+1);
948       for(int i = 0; i < 3; i++)
949         nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
950       nodeLs[3] = &verticesLs(0, e->getVertex(4)->getIndex());
951       bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
952                         tetras, quads, triangles, recur, nodeLs);
953       DI_Tetra T2(
954         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
955         e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
956         e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z(),
957         e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
958       // T2.setPolynomialOrder(recur+1);
959       nodeLs[0] = &verticesLs(0, e->getVertex(0)->getIndex());
960       for(int i = 1; i < 4; i++)
961         nodeLs[i] = &verticesLs(0, e->getVertex(i + 1)->getIndex());
962       bool iC2 = T2.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
963                         tetras, quads, triangles, recur, nodeLs);
964       isCut = iC1 || iC2;
965     }
966     else if(eType == TYPE_POLYH) {
967       for(int i = 0; i < e->getNumChildren(); i++) {
968         MTetrahedron *t = (MTetrahedron *)e->getChild(i);
969         DI_Tetra Tet(
970           t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
971           t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
972           t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
973           t->getVertex(3)->x(), t->getVertex(3)->y(), t->getVertex(3)->z());
974         Tet.setPolynomialOrder(recur + 1);
975         for(std::size_t i = 0; i < t->getNumVertices(); i++)
976           nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
977         bool iC = Tet.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
978                           tetras, quads, triangles, recur, nodeLs);
979         isCut = (isCut || iC);
980       }
981     }
982     MPolyhedron *p1 = nullptr, *p2 = nullptr;
983     if(isCut) {
984       std::vector<MTetrahedron *> poly[2];
985 
986       for(std::size_t i = nbTe; i < tetras.size(); i++) {
987         MVertex *mv[4] = {nullptr, nullptr, nullptr, nullptr};
988         for(int j = 0; j < 4; j++) {
989           int numV = getElementVertexNum(tetras[i]->pt(j), e);
990           if(numV == -1) {
991             MVertex *newv =
992               new MVertex(tetras[i]->x(j), tetras[i]->y(j), tetras[i]->z(j));
993             auto it = newVertices.insert(newv);
994             mv[j] = *(it.first);
995             if(!it.second) newv->deleteLast();
996           }
997           else {
998             auto it = vertexMap.find(numV);
999             if(it == vertexMap.end()) {
1000               mv[j] = new MVertex(tetras[i]->x(j), tetras[i]->y(j),
1001                                   tetras[i]->z(j), nullptr, numV);
1002               vertexMap[numV] = mv[j];
1003             }
1004             else
1005               mv[j] = it->second;
1006           }
1007         }
1008         MTetrahedron *mt = new MTetrahedron(mv[0], mv[1], mv[2], mv[3], 0, 0);
1009         if(tetras[i]->lsTag() < 0)
1010           poly[0].push_back(mt);
1011         else
1012           poly[1].push_back(mt);
1013       }
1014       bool own = (eParent && !e->ownsParent()) ? false : true;
1015       if(poly[0].size()) {
1016         int n = (e->getParent()) ? e->getNum() : ++numEle;
1017         p1 = new MPolyhedron(poly[0], n, ePart, own, parent);
1018         own = false;
1019         int reg = getElementaryTag(-1, elementary, newElemTags[3]);
1020         getPhysicalTag(-1, gePhysicals, newPhysTags[3]);
1021         elements[9][reg].push_back(p1);
1022         assignPhysicals(GM, gePhysicals, reg, 3, physicals, newPhysTags[3], -1);
1023       }
1024       if(poly[1].size()) {
1025         int n =
1026           (e->getParent() && poly[0].size() == 0) ? e->getNum() : ++numEle;
1027         p2 = new MPolyhedron(poly[1], n, ePart, own, parent);
1028         getPhysicalTag(1, gePhysicals, newPhysTags[3]);
1029         elements[9][elementary].push_back(p2);
1030         assignPhysicals(GM, gePhysicals, elementary, 3, physicals,
1031                         newPhysTags[3], 1);
1032       }
1033       // check for border surfaces cut earlier along the polyhedra
1034       auto itr = borders[1].equal_range(copy);
1035       std::vector<std::pair<MElement *, MElement *> > bords;
1036       for(auto it = itr.first; it != itr.second; it++) {
1037         MElement *tb = it->second;
1038         int match = 0;
1039         for(std::size_t i = 0; i < p1->getNumPrimaryVertices(); i++) {
1040           if(tb->getVertex(0) == p1->getVertex(i) ||
1041              tb->getVertex(1) == p1->getVertex(i) ||
1042              tb->getVertex(2) == p1->getVertex(i))
1043             match++;
1044           if(match == 3) break;
1045         }
1046         MElement *dom = (match == 3) ? p1 : p2;
1047         tb->setDomain(dom, (tb->getDomain(0) == copy) ? 0 : 1);
1048         bords.push_back(std::make_pair(dom, tb));
1049       }
1050       borders[1].erase(itr.first, itr.second);
1051       for(std::size_t i = 0; i < bords.size(); i++) borders[1].insert(bords[i]);
1052       if(eParent) {
1053         copy->setParent(nullptr, false);
1054         delete copy;
1055       }
1056     }
1057     else { // no cut
1058       int lsTag;
1059       if(eType == TYPE_HEX)
1060         lsTag = hexas[nbH]->lsTag();
1061       else
1062         lsTag = tetras[nbTe]->lsTag();
1063       int reg = getElementaryTag(lsTag, elementary, newElemTags[3]);
1064       getPhysicalTag(lsTag, gePhysicals, newPhysTags[3]);
1065       if(eType == TYPE_TET)
1066         elements[4][reg].push_back(copy);
1067       else if(eType == TYPE_HEX)
1068         elements[5][reg].push_back(copy);
1069       else if(eType == TYPE_PRI)
1070         elements[6][reg].push_back(copy);
1071       else if(eType == TYPE_PYR)
1072         elements[7][reg].push_back(copy);
1073       else if(eType == TYPE_POLYH)
1074         elements[9][reg].push_back(copy);
1075       assignPhysicals(GM, gePhysicals, reg, 3, physicals, newPhysTags[3],
1076                       lsTag);
1077     }
1078 
1079     for(std::size_t i = nbTr; i < triangles.size(); i++) {
1080       MVertex *mv[3] = {nullptr, nullptr, nullptr};
1081       for(int j = 0; j < 3; j++) {
1082         int numV = getElementVertexNum(triangles[i]->pt(j), e);
1083         if(numV == -1) {
1084           MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
1085                                       triangles[i]->z(j));
1086           auto it = newVertices.insert(newv);
1087           mv[j] = *(it.first);
1088           if(!it.second) newv->deleteLast();
1089         }
1090         else {
1091           auto it = vertexMap.find(numV);
1092           if(it == vertexMap.end()) {
1093             mv[j] = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
1094                                 triangles[i]->z(j), nullptr, numV);
1095             vertexMap[numV] = mv[j];
1096           }
1097           else
1098             mv[j] = it->second;
1099         }
1100       }
1101       MTriangle *tri;
1102       if(p1 || p2) {
1103         if(!p1)
1104           tri =
1105             new MTriangleBorder(mv[0], mv[1], mv[2], ++numEle, ePart, p2, p1);
1106         else
1107           tri =
1108             new MTriangleBorder(mv[0], mv[1], mv[2], ++numEle, ePart, p1, p2);
1109       }
1110       else
1111         tri = new MTriangle(mv[0], mv[1], mv[2], ++numEle, ePart);
1112       int lsTag = triangles[i]->lsTag();
1113       int cR = elements[2].count(lsTag) + elements[3].count(lsTag) +
1114                elements[8].count(lsTag);
1115       int cP = 0;
1116       for(auto it = physicals[2].begin(); it != physicals[2].end(); it++)
1117         for(auto it2 = it->second.begin(); it2 != it->second.end(); it2++)
1118           if(it2->first == lsTag) {
1119             cP = 1;
1120             break;
1121           }
1122       // the surfaces are cut before the volumes!
1123       int reg = getBorderTag(lsTag, cR, newElemTags[2][0], borderElemTags[1]);
1124       int physTag =
1125         (!gePhysicals.size()) ?
1126           0 :
1127           getBorderTag(lsTag, cP, newPhysTags[2][0], borderPhysTags[1]);
1128       elements[2][reg].push_back(tri);
1129       if(physTag) assignLsPhysical(GM, reg, 2, physicals, physTag, lsTag);
1130       for(int i = 0; i < 2; i++)
1131         if(tri->getDomain(i))
1132           borders[1].insert(std::make_pair(tri->getDomain(i), tri));
1133     }
1134   } break;
1135   case TYPE_TRI:
1136   case TYPE_QUA:
1137   case TYPE_POLYG: {
1138     if(eType == TYPE_TRI) {
1139       DI_Triangle T(
1140         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
1141         e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
1142         e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z());
1143       T.setPolynomialOrder(recur + 1);
1144       for(std::size_t i = 0; i < e->getNumVertices(); i++)
1145         nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
1146       isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
1147                     quads, triangles, lines, recur, nodeLs);
1148     }
1149     else if(eType == TYPE_QUA) {
1150       DI_Quad Q(
1151         e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
1152         e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
1153         e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
1154         e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
1155       Q.setPolynomialOrder(recur + 1);
1156       for(std::size_t i = 0; i < e->getNumVertices(); i++)
1157         nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
1158       isCut = Q.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
1159                     quads, triangles, lines, recur, nodeLs);
1160     }
1161     else if(eType == TYPE_POLYG) {
1162       for(int i = 0; i < e->getNumChildren(); i++) {
1163         MElement *t = e->getChild(i);
1164         DI_Triangle Tri(
1165           t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
1166           t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
1167           t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z());
1168         Tri.setPolynomialOrder(recur + 1);
1169         for(std::size_t i = 0; i < t->getNumVertices(); i++)
1170           nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
1171         bool iC = Tri.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
1172                           quads, triangles, lines, recur, nodeLs);
1173         isCut = (isCut || iC);
1174       }
1175     }
1176 
1177     MPolygon *p1 = nullptr, *p2 = nullptr;
1178     if(isCut) {
1179       std::vector<MTriangle *> poly[2];
1180 
1181       for(std::size_t i = nbTr; i < triangles.size(); i++) {
1182         MVertex *mv[3] = {nullptr, nullptr, nullptr};
1183         for(int j = 0; j < 3; j++) {
1184           int numV = getElementVertexNum(triangles[i]->pt(j), e);
1185           if(numV == -1) {
1186             MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
1187                                         triangles[i]->z(j));
1188             auto it = newVertices.insert(newv);
1189             mv[j] = *(it.first);
1190             if(!it.second) newv->deleteLast();
1191           }
1192           else {
1193             auto it = vertexMap.find(numV);
1194             if(it == vertexMap.end()) {
1195               mv[j] = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
1196                                   triangles[i]->z(j), nullptr, numV);
1197               vertexMap[numV] = mv[j];
1198             }
1199             else
1200               mv[j] = it->second;
1201           }
1202         }
1203         MTriangle *mt = new MTriangle(mv[0], mv[1], mv[2], 0, 0);
1204         if(triangles[i]->lsTag() < 0)
1205           poly[0].push_back(mt);
1206         else
1207           poly[1].push_back(mt);
1208       }
1209       // if quads
1210       for(std::size_t i = nbQ; i < quads.size(); i++) {
1211         MVertex *mv[4] = {nullptr, nullptr, nullptr, nullptr};
1212         for(int j = 0; j < 4; j++) {
1213           int numV = getElementVertexNum(quads[i]->pt(j), e);
1214           if(numV == -1) {
1215             MVertex *newv =
1216               new MVertex(quads[i]->x(j), quads[i]->y(j), quads[i]->z(j));
1217             auto it = newVertices.insert(newv);
1218             mv[j] = *(it.first);
1219             if(!it.second) newv->deleteLast();
1220           }
1221           else {
1222             auto it = vertexMap.find(numV);
1223             if(it == vertexMap.end()) {
1224               mv[j] = new MVertex(quads[i]->x(j), quads[i]->y(j),
1225                                   quads[i]->z(j), nullptr, numV);
1226               vertexMap[numV] = mv[j];
1227             }
1228             else
1229               mv[j] = it->second;
1230           }
1231         }
1232         MTriangle *mt0 = new MTriangle(mv[0], mv[1], mv[2], 0, 0);
1233         MTriangle *mt1 = new MTriangle(mv[0], mv[2], mv[3], 0, 0);
1234         if(quads[i]->lsTag() < 0) {
1235           poly[0].push_back(mt0);
1236           poly[0].push_back(mt1);
1237         }
1238         else {
1239           poly[1].push_back(mt0);
1240           poly[1].push_back(mt1);
1241         }
1242       }
1243 
1244       bool own = (eParent && !e->ownsParent()) ? false : true;
1245       if(poly[0].size()) {
1246         int n = (e->getParent()) ? e->getNum() : ++numEle;
1247         if(eType == MSH_TRI_B || eType == MSH_POLYG_B)
1248           p1 = new MPolygonBorder(poly[0], n, ePart, own, parent,
1249                                   copy->getDomain(0), copy->getDomain(1));
1250         else
1251           p1 = new MPolygon(poly[0], n, ePart, own, parent);
1252         own = false;
1253         int reg = getElementaryTag(-1, elementary, newElemTags[2]);
1254         getPhysicalTag(-1, gePhysicals, newPhysTags[2]);
1255         elements[8][reg].push_back(p1);
1256         assignPhysicals(GM, gePhysicals, reg, 2, physicals, newPhysTags[2], -1);
1257         for(int i = 0; i < 2; i++)
1258           if(p1->getDomain(i))
1259             borders[1].insert(std::make_pair(p1->getDomain(i), p1));
1260       }
1261       if(poly[1].size()) {
1262         int n =
1263           (e->getParent() && poly[0].size() == 0) ? e->getNum() : ++numEle;
1264         if(eType == MSH_TRI_B || eType == MSH_POLYG_B)
1265           p2 = new MPolygonBorder(poly[1], n, ePart, own, parent,
1266                                   copy->getDomain(0), copy->getDomain(1));
1267         else
1268           p2 = new MPolygon(poly[1], n, ePart, own, parent);
1269         getPhysicalTag(1, gePhysicals, newPhysTags[2]);
1270         elements[8][elementary].push_back(p2);
1271         assignPhysicals(GM, gePhysicals, elementary, 2, physicals,
1272                         newPhysTags[2], 1);
1273         for(int i = 0; i < 2; i++)
1274           if(p2->getDomain(i))
1275             borders[1].insert(std::make_pair(p2->getDomain(i), p2));
1276       }
1277       // check for border lines cut earlier along the polygons
1278       auto itr = borders[0].equal_range(copy);
1279       std::vector<std::pair<MElement *, MElement *> > bords;
1280       for(auto it = itr.first; it != itr.second; ++it) {
1281         MElement *lb = it->second;
1282         int match = 0;
1283         for(std::size_t i = 0; i < p1->getNumPrimaryVertices(); i++) {
1284           if(lb->getVertex(0) == p1->getVertex(i) ||
1285              lb->getVertex(1) == p1->getVertex(i))
1286             match++;
1287           if(match == 2) break;
1288         }
1289         MElement *dom = (match == 2) ? p1 : p2;
1290         lb->setDomain(dom, (lb->getDomain(0) == copy) ? 0 : 1);
1291         bords.push_back(std::make_pair(dom, lb));
1292       }
1293       borders[0].erase(itr.first, itr.second);
1294       for(std::size_t i = 0; i < bords.size(); i++) borders[0].insert(bords[i]);
1295       if(eParent) {
1296         copy->setParent(nullptr, false);
1297         delete copy;
1298       }
1299     }
1300     else { // no cut
1301       int lsTag;
1302       if(eType == TYPE_QUA)
1303         lsTag = quads[nbQ]->lsTag();
1304       else
1305         lsTag = triangles[nbTr]->lsTag();
1306       int reg = getElementaryTag(lsTag, elementary, newElemTags[2]);
1307       getPhysicalTag(lsTag, gePhysicals, newPhysTags[2]);
1308       if(eType == TYPE_TRI)
1309         elements[2][reg].push_back(copy);
1310       else if(eType == TYPE_QUA)
1311         elements[3][reg].push_back(copy);
1312       else if(eType == TYPE_POLYG)
1313         elements[8][reg].push_back(copy);
1314       assignPhysicals(GM, gePhysicals, reg, 2, physicals, newPhysTags[2],
1315                       lsTag);
1316       for(int i = 0; i < 2; i++)
1317         if(copy->getDomain(i))
1318           borders[1].insert(std::make_pair(copy->getDomain(i), copy));
1319     }
1320 
1321     for(std::size_t i = nbL; i < lines.size(); i++) {
1322       MVertex *mv[2] = {nullptr, nullptr};
1323       for(int j = 0; j < 2; j++) {
1324         int numV = getElementVertexNum(lines[i]->pt(j), e);
1325         if(numV == -1) {
1326           MVertex *newv =
1327             new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j));
1328           auto it = newVertices.insert(newv);
1329           mv[j] = *(it.first);
1330           if(!it.second) newv->deleteLast();
1331         }
1332         else {
1333           auto it = vertexMap.find(numV);
1334           if(it == vertexMap.end()) {
1335             mv[j] = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j),
1336                                 nullptr, numV);
1337             vertexMap[numV] = mv[j];
1338           }
1339           else
1340             mv[j] = it->second;
1341         }
1342       }
1343       MLine *lin;
1344       if(p1 || p2) {
1345         if(!p1)
1346           lin = new MLineBorder(mv[0], mv[1], ++numEle, ePart, p2, p1);
1347         else
1348           lin = new MLineBorder(mv[0], mv[1], ++numEle, ePart, p1, p2);
1349       }
1350       else
1351         lin = new MLine(mv[0], mv[1], ++numEle, ePart);
1352       int lsTag = lines[i]->lsTag();
1353       int cR = elements[1].count(lsTag);
1354       int cP = 0;
1355       for(auto it = physicals[1].begin(); it != physicals[1].end(); it++)
1356         for(auto it2 = it->second.begin(); it2 != it->second.end(); it2++)
1357           if(it2->first == lsTag) {
1358             cP = 1;
1359             break;
1360           }
1361       // the lines are cut before the surfaces!
1362       int reg = getBorderTag(lsTag, cR, newElemTags[1][0], borderElemTags[0]);
1363       int physTag =
1364         (!gePhysicals.size()) ?
1365           0 :
1366           getBorderTag(lsTag, cP, newPhysTags[1][0], borderPhysTags[0]);
1367       elements[1][reg].push_back(lin);
1368       if(physTag) assignLsPhysical(GM, reg, 1, physicals, physTag, lsTag);
1369       for(int i = 0; i < 2; i++)
1370         if(lin->getDomain(i))
1371           borders[0].insert(std::make_pair(lin->getDomain(i), lin));
1372     }
1373   } break;
1374   case TYPE_LIN: {
1375     DI_Line L(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
1376               e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z());
1377     L.setPolynomialOrder(recur + 1);
1378     for(std::size_t i = 0; i < e->getNumVertices(); i++)
1379       nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
1380     isCut = L.cut(RPN, ipV, cp, integOrder, lines, recur, nodeLs);
1381 
1382     if(isCut) {
1383       bool own = (eParent && !e->ownsParent()) ? false : true;
1384       for(std::size_t i = nbL; i < lines.size(); i++) {
1385         MVertex *mv[2] = {nullptr, nullptr};
1386         for(int j = 0; j < 2; j++) {
1387           int numV = getElementVertexNum(lines[i]->pt(j), e);
1388           if(numV == -1) {
1389             MVertex *newv =
1390               new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j));
1391             auto it = newVertices.insert(newv);
1392             mv[j] = *(it.first);
1393             if(!it.second) newv->deleteLast();
1394           }
1395           else {
1396             auto it = vertexMap.find(numV);
1397             if(it == vertexMap.end()) {
1398               mv[j] = new MVertex(lines[i]->x(j), lines[i]->y(j),
1399                                   lines[i]->z(j), nullptr, numV);
1400               vertexMap[numV] = mv[j];
1401             }
1402             else
1403               mv[j] = it->second;
1404           }
1405         }
1406         MLine *ml;
1407         int n = (e->getParent() && i == nbL) ? e->getNum() : ++numEle;
1408         if(eType != MSH_LIN_B)
1409           ml = new MLineChild(mv[0], mv[1], n, ePart, own, parent);
1410         else
1411           ml = new MLineBorder(mv[0], mv[1], n, ePart, copy->getDomain(0),
1412                                copy->getDomain(1));
1413         own = false;
1414         int reg =
1415           getElementaryTag(lines[i]->lsTag(), elementary, newElemTags[1]);
1416         int lsTag = lines[i]->lsTag();
1417         getPhysicalTag(lsTag, gePhysicals, newPhysTags[1]);
1418         elements[1][reg].push_back(ml);
1419         assignPhysicals(GM, gePhysicals, reg, 1, physicals, newPhysTags[1],
1420                         lsTag);
1421         for(int i = 0; i < 2; i++)
1422           if(ml->getDomain(i))
1423             borders[0].insert(std::make_pair(ml->getDomain(i), ml));
1424       }
1425       if(eParent) {
1426         copy->setParent(nullptr, false);
1427         delete copy;
1428       }
1429     }
1430     else { // no cut
1431       int lsTag = lines[nbL]->lsTag();
1432       int reg = getElementaryTag(lsTag, elementary, newElemTags[1]);
1433       getPhysicalTag(lsTag, gePhysicals, newPhysTags[1]);
1434       elements[1][reg].push_back(copy);
1435       assignPhysicals(GM, gePhysicals, reg, 1, physicals, newPhysTags[1],
1436                       lsTag);
1437       for(int i = 0; i < 2; i++)
1438         if(copy->getDomain(i))
1439           borders[0].insert(std::make_pair(copy->getDomain(i), copy));
1440     }
1441   } break;
1442   case TYPE_PNT: {
1443     DI_Point P(e->getVertex(0)->x(), e->getVertex(0)->y(),
1444                e->getVertex(0)->z());
1445     P.computeLs(RPN.back());
1446     int lsTag = P.lsTag();
1447     int reg = getElementaryTag(lsTag, elementary, newElemTags[0]);
1448     getPhysicalTag(lsTag, gePhysicals, newPhysTags[0]);
1449     elements[0][reg].push_back(copy);
1450     assignPhysicals(GM, gePhysicals, reg, 0, physicals, newPhysTags[0], lsTag);
1451   } break;
1452   default: Msg::Error("This type of element cannot be cut %d.", eType); break;
1453   }
1454 
1455   for(std::size_t i = 0; i < ipS.size(); i++) delete ipS[i];
1456   for(std::size_t i = 0; i < ipV.size(); i++) delete ipV[i];
1457   delete[] nodeLs;
1458 }
1459 
1460 #endif // HAVE_DINTEGRATION
1461 
buildCutMesh(GModel * gm,gLevelset * ls,std::map<int,std::vector<MElement * >> elements[10],std::map<int,MVertex * > & vertexMap,std::map<int,std::map<int,std::string>> physicals[4],bool cutElem)1462 GModel *buildCutMesh(GModel *gm, gLevelset *ls,
1463                      std::map<int, std::vector<MElement *> > elements[10],
1464                      std::map<int, MVertex *> &vertexMap,
1465                      std::map<int, std::map<int, std::string> > physicals[4],
1466                      bool cutElem)
1467 {
1468 #if defined(HAVE_DINTEGRATION)
1469 
1470   GModel *cutGM = new GModel(gm->getName() + "_cut");
1471   cutGM->setFileName(cutGM->getName());
1472 
1473   std::vector<GEntity *> gmEntities;
1474   gm->getEntities(gmEntities);
1475 
1476   std::vector<gLevelset *> primitives;
1477   ls->getPrimitivesPO(primitives);
1478   int primS = primitives.size();
1479   std::vector<gLevelset *> RPN;
1480   ls->getRPN(RPN);
1481   int numVert = gm->indexMeshVertices(true);
1482   int nbLs = (primS > 1) ? primS + 1 : 1;
1483   fullMatrix<double> verticesLs(nbLs, numVert + 1);
1484 
1485   // compute all at once for ls POINTS (type = 11)
1486   bool lsPoints = false;
1487   for(int i = 0; i < primS; i++)
1488     if(primitives[i]->type() == LSPOINTS) {
1489       lsPoints = true;
1490       break;
1491     }
1492   if(lsPoints) {
1493     std::vector<MVertex *> vert;
1494     for(std::size_t i = 0; i < gmEntities.size(); i++) {
1495       for(std::size_t j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
1496         vert.push_back(gmEntities[i]->getMeshVertex(j));
1497       }
1498     }
1499     for(int k = 0; k < primS; k++) {
1500       if(primitives[k]->type() == LSPOINTS) {
1501         ((gLevelsetPoints *)primitives[k])->computeLS(vert);
1502       }
1503     }
1504   }
1505 
1506   // compute and store levelset values + create new nodes
1507   for(std::size_t i = 0; i < gmEntities.size(); i++) {
1508     for(std::size_t j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
1509       MVertex *vi = gmEntities[i]->getMeshVertex(j);
1510       if(vi->getIndex() < 0) continue;
1511       int k = 0;
1512       for(; k < primS; k++)
1513         verticesLs(k, vi->getIndex()) =
1514           (*primitives[k])(vi->x(), vi->y(), vi->z());
1515       if(primS > 1)
1516         verticesLs(k, vi->getIndex()) = (*ls)(vi->x(), vi->y(), vi->z());
1517 
1518       MVertex *vn =
1519         new MVertex(vi->x(), vi->y(), vi->z(), nullptr, vi->getNum());
1520       vertexMap[vi->getNum()] = vn;
1521     }
1522   }
1523 
1524   // element number increment
1525   std::size_t numEle =
1526     gm->getNumMeshElements() + gm->getNumMeshParentElements();
1527   for(std::size_t i = 0; i < gmEntities.size(); i++) {
1528     for(std::size_t j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
1529       MElement *e = gmEntities[i]->getMeshElement(j);
1530       if(e->getNum() > numEle) numEle = e->getNum();
1531       if(e->getParent())
1532         if(e->getParent()->getNum() > numEle) numEle = e->getParent()->getNum();
1533     }
1534   }
1535 
1536   std::map<int, int> newElemTags[4]; // map<oldElementary,newElementary>[dim]
1537   std::map<int, int> newPhysTags[4]; // map<oldPhysical,newPhysical>[dim]
1538   for(int d = 0; d < 4; d++) {
1539     newElemTags[d][0] = gm->getMaxElementaryNumber(d); // max value at [dim][0]
1540     newPhysTags[d][0] = gm->getMaxPhysicalNumber(d); // max value at [dim][0]
1541   }
1542 
1543   std::map<MElement *, MElement *> newParents; // map<oldParent, newParent>
1544   std::map<MElement *, MElement *> newDomains; // map<oldDomain, newDomain>
1545   std::map<int, int>
1546     borderElemTags[2]; // map<lsTag,elementary>[line=0,surface=1]
1547   std::map<int, int> borderPhysTags[2]; // map<lstag,physical>[line=0,surface=1]
1548 
1549   // SplitMesh
1550   if(!cutElem) {
1551     for(std::size_t i = 0; i < gmEntities.size(); i++) {
1552       for(std::size_t j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
1553         MElement *e = gmEntities[i]->getMeshElement(j);
1554         e->setVolumePositive();
1555         elementSplitMesh(e, RPN, verticesLs, gmEntities[i], gm, numEle,
1556                          vertexMap, newParents, newDomains, elements, physicals,
1557                          newElemTags, newPhysTags, borderElemTags,
1558                          borderPhysTags);
1559         if(e->getPartition() > static_cast<int>(cutGM->getNumPartitions())) {
1560           cutGM->setNumPartitions(e->getPartition());
1561         }
1562       }
1563     }
1564     return cutGM;
1565   }
1566 
1567   // CutMesh
1568   newVerticesContainer newVertices;
1569   // multimap<domain,border>[polyg=0,polyh=1]
1570   std::multimap<MElement *, MElement *> borders[2];
1571   DI_Point::Container cp;
1572   std::vector<DI_Line *> lines;
1573   std::vector<DI_Triangle *> triangles;
1574   std::vector<DI_Quad *> quads;
1575   std::vector<DI_Tetra *> tetras;
1576   std::vector<DI_Hexa *> hexas;
1577   std::vector<int> lsLineRegs;
1578   for(std::size_t i = 0; i < gmEntities.size(); i++) {
1579     std::vector<int> oldLineRegs;
1580     for(auto it = elements[1].begin(); it != elements[1].end(); it++)
1581       oldLineRegs.push_back(it->first);
1582     int nbBorders = borders[0].size();
1583     for(std::size_t j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
1584       MElement *e = gmEntities[i]->getMeshElement(j);
1585       e->setVolumePositive();
1586       elementCutMesh(e, RPN, verticesLs, gmEntities[i], gm, numEle, vertexMap,
1587                      newVertices, newParents, newDomains, borders, elements,
1588                      physicals, newElemTags, newPhysTags, borderElemTags,
1589                      borderPhysTags, cp, lines, triangles, quads, tetras,
1590                      hexas);
1591       if(e->getPartition() > static_cast<int>(cutGM->getNumPartitions())) {
1592         cutGM->setNumPartitions(e->getPartition());
1593       }
1594     }
1595 
1596     // Create elementary and physical for non connected border lines
1597     if((int)borders[0].size() > nbBorders && gmEntities[i]->dim() == 2 &&
1598        i == gm->getNumVertices() + gm->getNumEdges() + gm->getNumFaces() - 1) {
1599       int k = 0;
1600       for(auto it = elements[1].begin(); it != elements[1].end(); it++) {
1601         if(oldLineRegs.size() && it->first == oldLineRegs[k])
1602           k++;
1603         else
1604           lsLineRegs.push_back(it->first);
1605       }
1606       for(std::size_t j = 0; j < lsLineRegs.size(); j++) {
1607         int nLR = lsLineRegs[j];
1608         bool havePhys = physicals[1][nLR].size();
1609         while(1) {
1610           std::vector<MElement *> conLines;
1611           conLines.push_back(elements[1][nLR][0]);
1612           elements[1][nLR].erase(elements[1][nLR].begin());
1613           MVertex *v1 = conLines[0]->getVertex(0);
1614           MVertex *v2 = conLines[0]->getVertex(1);
1615           for(std::size_t k = 0; k < elements[1][nLR].size();) {
1616             MVertex *va = elements[1][nLR][k]->getVertex(0);
1617             MVertex *vb = elements[1][nLR][k]->getVertex(1);
1618             if(va == v1 || vb == v1 || va == v2 || vb == v2) {
1619               conLines.push_back(elements[1][nLR][k]);
1620               elements[1][nLR].erase(elements[1][nLR].begin() + k);
1621               if(v1 == va)
1622                 v1 = vb;
1623               else if(v1 == vb)
1624                 v1 = va;
1625               else if(v2 == va)
1626                 v2 = vb;
1627               else if(v2 == vb)
1628                 v2 = va;
1629               k = 0;
1630             }
1631             else
1632               k++;
1633           }
1634           if(!elements[1][nLR].empty()) {
1635             int newReg = ++newElemTags[1][0];
1636             int newPhys = (!havePhys) ? 0 : ++newPhysTags[1][0];
1637             if(newPhys)
1638               assignLsPhysical(gm, newReg, 1, physicals, newPhys,
1639                                lines[lines.size() - 1]->lsTag());
1640             for(std::size_t k = 0; k < conLines.size(); k++)
1641               elements[1][newReg].push_back(conLines[k]);
1642           }
1643           else {
1644             for(std::size_t k = 0; k < conLines.size(); k++)
1645               elements[1][nLR].push_back(conLines[k]);
1646             break;
1647           }
1648         }
1649       }
1650     }
1651 
1652     for(auto it = cp.begin(); it != cp.end(); it++) delete *it;
1653     for(std::size_t k = 0; k < lines.size(); k++) delete lines[k];
1654     for(std::size_t k = 0; k < triangles.size(); k++) delete triangles[k];
1655     for(std::size_t k = 0; k < quads.size(); k++) delete quads[k];
1656     for(std::size_t k = 0; k < tetras.size(); k++) delete tetras[k];
1657     for(std::size_t k = 0; k < hexas.size(); k++) delete hexas[k];
1658     cp.clear();
1659     lines.clear();
1660     triangles.clear();
1661     quads.clear();
1662     tetras.clear();
1663     hexas.clear();
1664   }
1665 
1666 #if 0
1667   int numElements = 0;
1668   for(int i = 0; i < 10; i++) {
1669     printf(" - element type : %d\n", i);
1670     for(auto it = elements[i].begin(); it != elements[i].end(); it++){
1671       printf(" elementary : %d\n",it->first);
1672       for(std::size_t j = 0; j < it->second.size(); j++){
1673         MElement *e = it->second[j];
1674         printf("element %d",e->getNum());
1675         if(e->getParent()) printf(" par=%d (%d)",e->getParent()->getNum(),e->ownsParent());
1676         if(e->getDomain(0)) printf(" d0=%d",e->getDomain(0)->getNum());
1677         if(e->getDomain(1)) printf(" d1=%d",e->getDomain(1)->getNum());
1678         printf("\n"); numElements++;
1679       }
1680     }
1681   }
1682   printf("PHYS\n");
1683   for(int i = 0; i < 4; i++)
1684     for(auto it = physicals[i].begin(); it != physicals[i].end(); it++)
1685       for(auto it2 = it->second.begin(); it2 != it->second.end(); it2++)
1686         printf(" dim=%d reg=%d phys=%d \"%s\"\n", i, it->first, it2->first,
1687                it2->second.c_str());
1688   printf("new Model : %d elements %d nodes\n\n", numElements, vertexMap.size());
1689 #endif
1690 
1691   for(auto it = newVertices.begin(); it != newVertices.end(); ++it) {
1692     vertexMap[(*it)->getNum()] = *it;
1693   }
1694 
1695   return cutGM;
1696 #else
1697   Msg::Error("Gmsh need to be compiled with levelset support to cut mesh");
1698   return 0;
1699 #endif
1700 }
1701