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