gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
MElementCut.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // 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 
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 
84 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 
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 
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 
257 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 
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 
321 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 
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 
376 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 
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 
435 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 
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
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
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 
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 }
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 
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 
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 
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 
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 
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 
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 }
MPolygonBorder
Definition: MElementCut.h:461
MElement::getNum
virtual std::size_t getNum() const
Definition: MElement.h:68
MTriangleBorder
Definition: MElementCut.h:425
MLineBorder::getParent
virtual MElement * getParent() const
Definition: MElementCut.h:518
MElement::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MElement.cpp:1128
MTetrahedron::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MTetrahedron.cpp:121
MPolygon::_edges
std::vector< MEdge > _edges
Definition: MElementCut.h:206
MElement::setVolumePositive
virtual bool setVolumePositive()
Definition: MElement.cpp:591
MLineChild::_intpt
IntPt * _intpt
Definition: MElementCut.h:366
MEdge
Definition: MEdge.h:14
GModel::getNumMeshParentElements
std::size_t getNumMeshParentElements() const
Definition: GModel.cpp:1551
dot
double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
Definition: STensor3.h:71
MLineBorder::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MElementCut.cpp:454
MElement::getDomain
virtual MElement * getDomain(int i) const
Definition: MElement.h:243
GModel::getMaxElementaryNumber
int getMaxElementaryNumber(int dim)
Definition: GModel.cpp:817
MTriangle::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3]) const
Definition: MTriangle.cpp:124
TYPE_LIN
#define TYPE_LIN
Definition: GmshDefines.h:65
MTetrahedron
Definition: MTetrahedron.h:34
MLine::getVertex
virtual MVertex * getVertex(int num)
Definition: MLine.h:45
Msg::Info
static void Info(const char *fmt,...)
Definition: GmshMessage.cpp:587
MPolygon::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MElementCut.cpp:279
MLineChild
Definition: MElementCut.h:362
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
MTriangle::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MTriangle.cpp:353
MVertex
Definition: MVertex.h:24
MElement::getDim
virtual int getDim() const =0
GModel::getName
std::string getName() const
Definition: GModel.h:329
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
buildCutMesh
GModel * 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)
Definition: MElementCut.cpp:1462
GEntity::physicals
std::vector< int > physicals
Definition: GEntity.h:65
MPolyhedron::getEdge
virtual MEdge getEdge(int num) const
Definition: MElementCut.h:78
SPoint3
Definition: SPoint3.h:14
MLineBorder::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MElementCut.cpp:435
MPolygon::_intpt
IntPt * _intpt
Definition: MElementCut.h:202
MElement::getParent
virtual MElement * getParent() const
Definition: MElement.h:231
MVertex::getNum
std::size_t getNum() const
Definition: MVertex.h:86
MFace::normal
SVector3 normal() const
Definition: MFace.cpp:204
TYPE_PNT
#define TYPE_PNT
Definition: GmshDefines.h:64
MPolygon::getVertex
virtual MVertex * getVertex(int num)
Definition: MElementCut.h:240
SVector3
Definition: SVector3.h:16
TYPE_TRI
#define TYPE_TRI
Definition: GmshDefines.h:66
MPolyhedron::getVolume
virtual double getVolume()
Definition: MElementCut.h:122
MTriangle::getVertex
virtual MVertex * getVertex(int num)
Definition: MTriangle.h:62
MLineBorder
Definition: MElementCut.h:495
MPolyhedron::getFace
virtual MFace getFace(int num) const
Definition: MElementCut.h:96
MTriangleBorder::_intpt
IntPt * _intpt
Definition: MElementCut.h:428
MPolygon::_parts
std::vector< MTriangle * > _parts
Definition: MElementCut.h:203
GModel::getNumMeshElements
std::size_t getNumMeshElements(int dim=-1) const
Definition: GModel.cpp:1540
MElement::copy
virtual MElement * copy(std::map< int, MVertex * > &vertexMap, std::map< MElement *, MElement * > &newParents, std::map< MElement *, MElement * > &newDomains)
Definition: MElement.cpp:2486
GEntity
Definition: GEntity.h:31
IntPt::pt
double pt[3]
Definition: GaussIntegration.h:13
TYPE_PRI
#define TYPE_PRI
Definition: GmshDefines.h:70
MTetrahedron::getNumVertices
virtual std::size_t getNumVertices() const
Definition: MTetrahedron.h:66
MPolyhedron::_init
void _init()
Definition: MElementCut.cpp:25
GModel::getNumPartitions
std::size_t getNumPartitions() const
Definition: GModel.h:602
MElement::getNumVertices
virtual std::size_t getNumVertices() const =0
GModel::setNumPartitions
void setNumPartitions(std::size_t npart)
Definition: GModel.h:603
getNGQLPts
int getNGQLPts(int order)
Definition: GaussQuadratureLin.cpp:33
MElement::ownsParent
virtual bool ownsParent() const
Definition: MElement.h:236
MLine::_v
MVertex * _v[2]
Definition: MLine.h:23
MLine
Definition: MLine.h:21
gLevelsetPoints
Definition: gmshLevelset.h:217
MLineBorder::_intpt
IntPt * _intpt
Definition: MElementCut.h:498
MElement::getVertex
virtual const MVertex * getVertex(int num) const =0
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
fullMatrix< double >
MElement::getEdge
virtual MEdge getEdge(int num) const =0
MFace
Definition: MFace.h:20
GModel::getPhysicalName
std::string getPhysicalName(int dim, int num) const
Definition: GModel.cpp:961
IntPt::weight
double weight
Definition: GaussIntegration.h:14
getNGQTPts
int getNGQTPts(int order, bool forceTensorRule=false)
Definition: GaussQuadratureTri.cpp:904
MEdge::getVertex
MVertex * getVertex(std::size_t i) const
Definition: MEdge.h:39
MTriangleBorder::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MElementCut.cpp:376
GModel::setFileName
void setFileName(const std::string &fileName)
Definition: GModel.cpp:123
MElement::getType
virtual int getType() const =0
gLevelset
Definition: gmshLevelset.h:64
MElementCut.h
MElement::getNumChildren
virtual int getNumChildren() const
Definition: MElement.h:234
MLine::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MLine.cpp:15
MElement::getNumFaces
virtual int getNumFaces()=0
GModel::getMaxPhysicalNumber
int getMaxPhysicalNumber(int dim)
Definition: GModel.cpp:917
MLineChild::_orig
MElement * _orig
Definition: MElementCut.h:365
MTetrahedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MTetrahedron.h:67
MTriangleBorder::getParent
virtual MElement * getParent() const
Definition: MElementCut.h:448
getNGQTetPts
int getNGQTetPts(int order, bool forceTensorRule=false)
Definition: GaussQuadratureTet.cpp:3362
GModel::getNumVertices
std::size_t getNumVertices() const
Definition: GModel.h:347
MPolygon::_innerVertices
std::vector< MVertex * > _innerVertices
Definition: MElementCut.h:205
MElement::getFace
virtual MFace getFace(int num) const =0
GModel
Definition: GModel.h:44
MSH_TRI_B
#define MSH_TRI_B
Definition: GmshDefines.h:147
LSPOINTS
#define LSPOINTS
Definition: gmshLevelset.h:48
MLine::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MLine.h:91
MTriangle::_v
MVertex * _v[3]
Definition: MTriangle.h:28
MFace::getVertex
MVertex * getVertex(std::size_t i) const
Definition: MFace.h:30
Tet
Definition: delaunay3d.cpp:179
MPolygon::getParent
virtual MElement * getParent() const
Definition: MElementCut.h:293
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
MLineChild::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MElementCut.cpp:340
MSH_POLYG_B
#define MSH_POLYG_B
Definition: GmshDefines.h:148
MElement
Definition: MElement.h:30
GEntity::tag
int tag() const
Definition: GEntity.h:280
MPolyhedron
Definition: MElementCut.h:21
GModel::getNumFaces
std::size_t getNumFaces() const
Definition: GModel.h:345
MSH_LIN_B
#define MSH_LIN_B
Definition: GmshDefines.h:146
MPolyhedron::_intpt
IntPt * _intpt
Definition: MElementCut.h:25
MTriangle::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MTriangle.h:155
GModel::getEntities
void getEntities(std::vector< GEntity * > &entities, int dim=-1) const
Definition: GModel.cpp:651
MVertex::deleteLast
void deleteLast()
Definition: MVertex.cpp:41
MLineChild::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MElementCut.cpp:321
MPolygon::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MElementCut.cpp:257
gLevelset::getPrimitivesPO
void getPrimitivesPO(std::vector< gLevelset * > &primitives)
Definition: gmshLevelset.cpp:291
MPolyhedron::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MElementCut.cpp:84
GModel::getNumEdges
std::size_t getNumEdges() const
Definition: GModel.h:346
MTriangle
Definition: MTriangle.h:26
TYPE_PYR
#define TYPE_PYR
Definition: GmshDefines.h:69
MElement::pnt
virtual void pnt(double u, double v, double w, SPoint3 &p) const
Definition: MElement.cpp:1072
MPolygon::getNumEdges
virtual int getNumEdges() const
Definition: MElementCut.h:252
TYPE_QUA
#define TYPE_QUA
Definition: GmshDefines.h:67
MPolyhedron::_orig
MElement * _orig
Definition: MElementCut.h:24
TYPE_POLYG
#define TYPE_POLYG
Definition: GmshDefines.h:72
MVertex::getIndex
long int getIndex() const
Definition: MVertex.h:93
IntPt
Definition: GaussIntegration.h:12
MTriangleBorder::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MElementCut.cpp:396
gLevelset::getRPN
void getRPN(std::vector< gLevelset * > &gLsRPN)
Definition: gmshLevelset.cpp:321
fullMatrix::size1
int size1() const
Definition: fullMatrix.h:274
MElement::getChild
virtual MElement * getChild(int i) const
Definition: MElement.h:235
z
const double z
Definition: GaussQuadratureQuad.cpp:56
MQuadrangle.h
MElement::getNumPrimaryVertices
std::size_t getNumPrimaryVertices() const
Definition: MElement.h:160
MElement.h
picojson::copy
void copy(const std::string &s, Iter oi)
Definition: picojson.h:510
MPolyhedron::_parts
std::vector< MTetrahedron * > _parts
Definition: MElementCut.h:26
TYPE_HEX
#define TYPE_HEX
Definition: GmshDefines.h:71
TYPE_TET
#define TYPE_TET
Definition: GmshDefines.h:68
MPolygon::_vertices
std::vector< MVertex * > _vertices
Definition: MElementCut.h:204
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
MPolyhedron::getVertex
virtual MVertex * getVertex(int num)
Definition: MElementCut.h:65
MElement::getJacobian
virtual double getJacobian(const fullMatrix< double > &gsf, double jac[3][3]) const
Definition: MElement.cpp:868
MPolygon
Definition: MElementCut.h:198
MPolyhedron::_vertices
std::vector< MVertex * > _vertices
Definition: MElementCut.h:27
GModel.h
MFace::getNumVertices
std::size_t getNumVertices() const
Definition: MFace.h:29
MVertex::y
double y() const
Definition: MVertex.h:61
MPolyhedron::_edges
std::vector< MEdge > _edges
Definition: MElementCut.h:29
MElement::getPolynomialOrder
virtual int getPolynomialOrder() const
Definition: MElement.h:78
MPolyhedron::_faces
std::vector< MFace > _faces
Definition: MElementCut.h:30
MElement::getPartition
virtual int getPartition() const
Definition: MElement.h:92
MTetrahedron::isInside
virtual bool isInside(double u, double v, double w) const
Definition: MTetrahedron.h:174
MElement::getNumEdges
virtual int getNumEdges() const =0
MPolyhedron::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MElementCut.cpp:107
TYPE_POLYH
#define TYPE_POLYH
Definition: GmshDefines.h:73
MQuadrangle
Definition: MQuadrangle.h:26
gmshLevelset.h
MPolygon::_orig
MElement * _orig
Definition: MElementCut.h:201
GModel::indexMeshVertices
std::size_t indexMeshVertices(bool all, int singlePartition=0, bool renumber=true)
Definition: GModel.cpp:2135
MPolyhedron::_innerVertices
std::vector< MVertex * > _innerVertices
Definition: MElementCut.h:28
MEdge::isInside
bool isInside(MVertex *v) const
Definition: MEdge.cpp:16
MVertex::x
double x() const
Definition: MVertex.h:60
MElement::setDomain
virtual void setDomain(MElement *e, int i)
Definition: MElement.h:244
MPolygon::_initVertices
void _initVertices()
Definition: MElementCut.cpp:152