diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 92b5f3d4b68d4fdc531879ab24b4df0acc063f29..b9d91aa18cd1d1a92936c7fb6215700c63eb931a 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -124,8 +124,9 @@
#define MSH_LIN_B 67
#define MSH_TRI_B 68
#define MSH_POLYG_B 69
+#define MSH_LIN_C 70
-#define MSH_NUM_TYPE 35
+#define MSH_NUM_TYPE 70
// Geometric entities
#define ENT_NONE 0
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 2023e9272e3b2375054b1ca7d416add7e86150e1..ec453a6abb4e2f859b2c031d24a3c2096e0d3783 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -621,16 +621,20 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
// get the number of elements we need to save
int numElements = getNumElementsMSH(this, saveAll, saveSinglePartition);
- // get parent elements for polygons and polyhedra
- std::map<MElement*, GEntity*> parents[2];
+ // get parent elements for lines, polygons and polyhedra
+ std::map<MElement*, GEntity*> parents[3];
+ for(eiter it = firstEdge(); it != lastEdge(); ++it)
+ for(unsigned int i = 0; i < (*it)->lines.size(); i++)
+ if((*it)->lines[i]->ownsParent())
+ parents[0][(*it)->lines[i]->getParent()] = (*it);
for(fiter it = firstFace(); it != lastFace(); ++it)
for(unsigned int i = 0; i < (*it)->polygons.size(); i++)
if((*it)->polygons[i]->ownsParent())
- parents[0][(*it)->polygons[i]->getParent()] = (*it);
+ parents[1][(*it)->polygons[i]->getParent()] = (*it);
for(riter it = firstRegion(); it != lastRegion(); ++it)
for(unsigned int i = 0; i < (*it)->polyhedra.size(); i++)
if((*it)->polyhedra[i]->ownsParent())
- parents[1][(*it)->polyhedra[i]->getParent()] = (*it);
+ parents[2][(*it)->polyhedra[i]->getParent()] = (*it);
if(version >= 2.0){
fprintf(fp, "$MeshFormat\n");
@@ -692,13 +696,21 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version, binary, num, (*it)->tag(), (*it)->physicals);
// lines
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin();
+ it != parents[0].end(); it++)
+ if(it->first->getType() == TYPE_LIN) {
+ writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+ it->second->tag(), it->second->physicals, 0);
+ parentsNum[it->first] = num;
+ }
for(eiter it = firstEdge(); it != lastEdge(); ++it)
writeElementsMSH(fp, this, (*it)->lines, saveAll, saveSinglePartition,
- version, binary, num, (*it)->tag(), (*it)->physicals);
+ version, binary, num, (*it)->tag(), (*it)->physicals,
+ &parentsNum);
// triangles
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin();
- it != parents[0].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+ it != parents[1].end(); it++)
if(it->first->getType() == TYPE_TRI) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
@@ -709,8 +721,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version, binary, num, (*it)->tag(), (*it)->physicals);
// quads
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin();
- it != parents[0].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+ it != parents[1].end(); it++)
if(it->first->getType() == TYPE_QUA) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
@@ -721,8 +733,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version, binary, num, (*it)->tag(), (*it)->physicals);
// polygons
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin();
- it != parents[0].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+ it != parents[1].end(); it++)
if(it->first->getType() == TYPE_POLYG) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
@@ -734,8 +746,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
&parentsNum);
// tets
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
- it != parents[1].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
+ it != parents[2].end(); it++)
if(it->first->getType() == TYPE_TET) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
@@ -746,8 +758,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version, binary, num, (*it)->tag(), (*it)->physicals);
// hexas
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
- it != parents[1].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
+ it != parents[2].end(); it++)
if(it->first->getType() == TYPE_HEX) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
@@ -758,8 +770,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version, binary, num, (*it)->tag(), (*it)->physicals);
// prisms
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
- it != parents[1].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
+ it != parents[2].end(); it++)
if(it->first->getType() == TYPE_PRI) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
@@ -770,8 +782,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version, binary, num, (*it)->tag(), (*it)->physicals);
// pyramids
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
- it != parents[1].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
+ it != parents[2].end(); it++)
if(it->first->getType() == TYPE_PYR) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
@@ -782,8 +794,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version, binary, num, (*it)->tag(), (*it)->physicals);
// polyhedra
- for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
- it != parents[1].end(); it++)
+ for(std::map<MElement*, GEntity*>::const_iterator it = parents[2].begin();
+ it != parents[2].end(); it++)
if(it->first->getType() == TYPE_POLYH) {
writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
it->second->tag(), it->second->physicals, 0);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 8d055784ee716ec86d862cd884d44c2bbc88e8bc..e1d8f5d555034c1908b018edf5440dd1848976bb 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -776,6 +776,7 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_LIN_10 : if(name) *name = "Line 10"; return 2 + 8;
case MSH_LIN_11 : if(name) *name = "Line 11"; return 2 + 9;
case MSH_LIN_B : if(name) *name = "Line Border"; return 2;
+ case MSH_LIN_C : if(name) *name = "Line Child"; return 2;
case MSH_TRI_3 : if(name) *name = "Triangle 3"; return 3;
case MSH_TRI_6 : if(name) *name = "Triangle 6"; return 3 + 3;
case MSH_TRI_9 : if(name) *name = "Triangle 9"; return 3 + 6;
@@ -926,6 +927,7 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_LIN_10: return new MLineN(v, num, part);
case MSH_LIN_11: return new MLineN(v, num, part);
case MSH_LIN_B: return new MLineBorder(v, num, part);
+ case MSH_LIN_C: return new MLineChild(v, num, part);
case MSH_TRI_3: return new MTriangle(v, num, part);
case MSH_TRI_6: return new MTriangle6(v, num, part);
case MSH_TRI_9: return new MTriangleN(v, 3, num, part);
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 388746427fbae8bfe57b40453ea96d59362fd1a2..d87e34eaaff2682abf29cb905c5e7323dafc6233 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -86,11 +86,22 @@ void MPolyhedron::_init()
bool MPolyhedron::isInside(double u, double v, double w)
{
- double ksi[3] = {u, v, w};
+ double uvw[3] = {u, v, w};
for(unsigned int i = 0; i < _parts.size(); i++) {
- double uvw[3];
- _parts[i]->xyz2uvw(ksi,uvw);
- if(_parts[i]->isInside(uvw[0],uvw[1],uvw[2]))
+ double verts[4][3];
+ for(int j = 0; j < 4; j++) {
+ MVertex *vij = _parts[i]->getVertex(j);
+ double v_xyz[3] = {vij->x(), vij->y(), vij->z()};
+ _orig->xyz2uvw(v_xyz, verts[j]);
+ }
+ MVertex v0(verts[0][0], verts[0][1], verts[0][2]);
+ MVertex v1(verts[1][0], verts[1][1], verts[1][2]);
+ MVertex v2(verts[2][0], verts[2][1], verts[2][2]);
+ MVertex v3(verts[3][0], verts[3][1], verts[3][2]);
+ MTetrahedron t(&v0, &v1, &v2, &v3);
+ double ksi[3];
+ t.xyz2uvw(uvw, ksi);
+ if(t.isInside(ksi[0], ksi[1], ksi[2]))
return true;
}
return false;
@@ -176,11 +187,21 @@ void MPolygon::_initVertices()
}
bool MPolygon::isInside(double u, double v, double w)
{
- double ksi[3] = {u, v, w};
+ double uvw[3] = {u, v, w};
for(unsigned int i = 0; i < _parts.size(); i++) {
- double uvw[3];
- _parts[i]->xyz2uvw(ksi,uvw);
- if(_parts[i]->isInside(uvw[0],uvw[1],uvw[2]))
+ double v_uvw[3][3];
+ for(int j = 0; j < 3; j++) {
+ MVertex *vij = _parts[i]->getVertex(j);
+ double v_xyz[3] = {vij->x(), vij->y(), vij->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[j]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
+ MTriangle t(&v0, &v1, &v2);
+ double ksi[3];
+ t.xyz2uvw(uvw, ksi);
+ if(t.isInside(ksi[0], ksi[1], ksi[2]))
return true;
}
return false;
@@ -223,6 +244,61 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*pts = _intpt;
}
+//----------------------------------- MLineChild ------------------------------
+
+bool MLineChild::isInside(double u, double v, double w)
+{
+ double uvw[3] = {u, v, w};
+ double v_uvw[2][3];
+ for(int i = 0; i < 2; i++) {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MLine l(&v0, &v1);
+ double ksi[3];
+ l.xyz2uvw(uvw, ksi);
+ if(l.isInside(ksi[0], ksi[1], ksi[2]))
+ return true;
+ return false;
+}
+
+void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+ *npts = 0;
+ double jac[3][3];
+ if(_intpt) delete [] _intpt;
+ _intpt = new IntPt[getNGQLPts(pOrder)];
+ int nptsi;
+ IntPt *ptsi;
+ double v_uvw[2][3];
+ for(int i = 0; i < 2; i++) {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MLine l(&v0, &v1);
+ l.getIntegrationPoints(pOrder, &nptsi, &ptsi);
+ for(int ip = 0; ip < nptsi; ip++){
+ const double u = ptsi[ip].pt[0];
+ const double v = ptsi[ip].pt[1];
+ const double w = ptsi[ip].pt[2];
+ const double weight = ptsi[ip].weight;
+ const double detJ = l.getJacobian(u, v, w, jac);
+ SPoint3 p; l.pnt(u, v, w, p);
+ _intpt[*npts + ip].pt[0] = p.x();
+ _intpt[*npts + ip].pt[1] = p.y();
+ _intpt[*npts + ip].pt[2] = p.z();
+ _intpt[*npts + ip].weight = detJ * weight;
+ }
+ *npts = nptsi;
+ *pts = _intpt;
+}
+
//----------------------------------- MTriangleBorder ------------------------------
void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
@@ -636,9 +712,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
else
poly[1].push_back(mt);
}
+ bool own = (eParent && !e->ownsParent()) ? false : true;
if(poly[0].size()) {
- bool own = (eParent && !e->ownsParent()) ? false : true;
p1 = new MPolyhedron(poly[0], ++numEle, ePart, own, parent);
+ own = false;
int reg = getElementaryTag(-1, elementary, newElemTags[3]);
std::vector<int> phys;
getPhysicalTag(-1, gePhysicals, phys, newPhysTags[3]);
@@ -646,7 +723,6 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
assignPhysicals(GM, phys, reg, 3, physicals);
}
if(poly[1].size()) {
- bool own = poly[0].size() ? false : (eParent && !e->ownsParent()) ? false : true;
p2 = new MPolyhedron(poly[1], ++numEle, ePart, own, parent);
elements[9][elementary].push_back(p2);
assignPhysicals(GM, gePhysicals, elementary, 3, physicals);
@@ -794,9 +870,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
else
poly[1].push_back(mt);
}
+ bool own = (eParent && !e->ownsParent()) ? false : true;
if(poly[0].size()) {
- bool own = (eParent && !e->ownsParent()) ? false : true;
p1 = new MPolygon(poly[0], ++numEle, ePart, own, parent);
+ own = false;
int reg = getElementaryTag(-1, elementary, newElemTags[2]);
std::vector<int> phys;
getPhysicalTag(-1, gePhysicals, phys, newPhysTags[2]);
@@ -807,7 +884,6 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
borders[1].insert(std::pair<MElement*, MElement*>(p1->getDomain(i), p1));
}
if(poly[1].size()) {
- bool own = poly[0].size() ? false : (eParent && !e->ownsParent()) ? false : true;
p2 = new MPolygon(poly[1], ++numEle, ePart, own, parent);
elements[8][elementary].push_back(p2);
assignPhysicals(GM, gePhysicals, elementary, 2, physicals);
@@ -895,7 +971,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
isCut = L.cut(RPN, ipV, cp, integOrder, lines, 0, nodeLs);
if(isCut) {
- if(eType == MSH_LIN_2) numEle--; // delete copy
+ bool own = (eParent && !e->ownsParent()) ? false : true;
for (unsigned int i = nbL; i < lines.size(); i++){
MVertex *mv[2] = {NULL, NULL};
for(int j = 0; j < 2; j++){
@@ -916,9 +992,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
}
MLine *ml;
- if(eType != MSH_LIN_B) ml = new MLine(mv[0], mv[1], ++numEle, ePart);
+ if(eType != MSH_LIN_B) ml = new MLineChild(mv[0], mv[1], ++numEle, ePart, own, parent);
else ml = new MLineBorder(mv[0], mv[1], ++numEle, ePart,
copy->getDomain(0), copy->getDomain(1));
+ own = false;
int reg = getElementaryTag(lines[i]->lsTag(), elementary, newElemTags[1]);
std::vector<int> phys;
getPhysicalTag(lines[i]->lsTag(), gePhysicals, phys, newPhysTags[1]);
@@ -928,7 +1005,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
if(ml->getDomain(i))
borders[0].insert(std::pair<MElement*, MElement*>(ml->getDomain(i), ml));
}
- delete copy;
+ if(eParent) {copy->setParent(NULL, false); delete copy;}
}
else { // no cut
int reg = getElementaryTag(lines[nbL]->lsTag(), elementary, newElemTags[1]);
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 65f7a9d47e12e9db113caf362cd76109315c5dc3..6af8126073f432b25c90c1007bf02a71400a86a1 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -126,6 +126,8 @@ class MPolyhedron : public MElement {
{
_orig->getGradShapeFunctions(u, v, w, s, o);
}
+ // the parametric coordinates of the polyhedron are
+ // the coordinates in the local parent element.
virtual void xyz2uvw(double xyz[3], double uvw[3])
{
_orig->xyz2uvw(xyz,uvw);
@@ -245,6 +247,8 @@ class MPolygon : public MElement {
{
_orig->getGradShapeFunctions(u, v, w, s, o);
}
+ // the parametric coordinates of the polygon are
+ // the coordinates in the local parent element.
virtual void xyz2uvw(double xyz[3], double uvw[3])
{
_orig->xyz2uvw(xyz,uvw);
@@ -268,6 +272,51 @@ class MPolygon : public MElement {
}
};
+class MLineChild : public MLine {
+ protected:
+ bool _owner;
+ MElement* _orig;
+ IntPt *_intpt;
+ public:
+ MLineChild(MVertex *v0, MVertex *v1, int num = 0, int part = 0,
+ bool owner = false, MElement* orig = NULL)
+ : MLine(v0, v1, num, part), _owner(owner), _orig(orig) {}
+ MLineChild(std::vector<MVertex*> v, int num = 0, int part = 0,
+ bool owner = false, MElement* orig = NULL)
+ : MLine(v, num, part), _owner(owner), _orig(orig) {}
+ ~MLineChild()
+ {
+ if(_owner)
+ delete _orig;
+ }
+ virtual int getTypeForMSH() const { return MSH_LIN_C; }
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const
+ {
+ return _orig->getFunctionSpace(order);
+ }
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+ {
+ _orig->getShapeFunctions(u, v, w, s, o);
+ }
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+ {
+ _orig->getGradShapeFunctions(u, v, w, s, o);
+ }
+ // the parametric coordinates of the LineChildren are
+ // the coordinates in the local parent element.
+ virtual void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ _orig->xyz2uvw(xyz,uvw);
+ }
+ virtual bool isInside(double u, double v, double w);
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+ virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
+ virtual bool ownsParent() const { return _owner; }
+};
+
+// -------------------- Border classes
+
class MTriangleBorder : public MTriangle {
protected:
MElement* _domains[2];
@@ -401,6 +450,16 @@ class MLineBorder : public MLine {
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
+// Build a new GModel with elements on each side of the levelset ls.
+// New physical and elementary entities are created.
+// The physical and elementary numbers of the elements with ls < 0 are
+// the physical and elementary number of the elements cut.
+// The physical and elementary numbers of the elements with ls > 0 are
+// the maximum physical and elementary numbers existing in their dimension + 1.
+// The physical and elementary numbers of the elements on the border (ls=0) are
+// the levelset tag, unless an entity of the same dimension has already this number,
+// knowing that the elements are cut in ascending dimension order (points, lines,
+// surfaces and then volumes).
GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::map<int, std::vector<MElement*> > elements[10],
std::map<int, MVertex*> &vertexMap,