better handling of embedded entities in extrusions, to avoid creating duplicate mesh nodes

CHANGELOG.txt

 4.5.5 (Work-in-progress): tooltips in GUI to help discovery of scripting
 options; fixed MED IO of high-order elements; fixed OCC attribute search by
 bounding box; fix parsing of mac-encoded scripts; new RecombineMesh command;
-small bug fixes.
+added support for extrusion of mixed-dimension entities with OCC; small bug
+fixes.
 
 4.5.4 (February 29, 2020): periodic mesh optimization now ensures that the
 master mesh is not modified; code cleanup; API tutorials; small bug fixes.

Geo/GFace.cpp

@@ -336,8 +336,8 @@ std::vector<MVertex *> GFace::getEmbeddedMeshVertices(bool force) const
       tmp.insert((*it)->getEndVertex()->mesh_vertices.begin(),
                  (*it)->getEndVertex()->mesh_vertices.end());
   }
-  for(std::set<GVertex *>::const_iterator it = embedded_vertices.begin();
-      it != embedded_vertices.end(); it++) {
+  for(std::set<GVertex *, GEntityPtrLessThan>::const_iterator it =
+        embedded_vertices.begin(); it != embedded_vertices.end(); it++) {
     tmp.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
   }
   return std::vector<MVertex *>(tmp.begin(), tmp.end());
@@ -416,8 +416,8 @@ std::string GFace::getAdditionalInfoString(bool multline)
 
   if(embedded_vertices.size()) {
     sstream << "Embedded points: ";
-    for(std::set<GVertex *>::iterator it = embedded_vertices.begin();
-        it != embedded_vertices.end(); ++it) {
+    for(std::set<GVertex *, GEntityPtrLessThan>::iterator it =
+          embedded_vertices.begin(); it != embedded_vertices.end(); ++it) {
       if(it != embedded_vertices.begin()) sstream << ", ";
       sstream << (*it)->tag();
     }
@@ -480,8 +480,8 @@ void GFace::writeGEO(FILE *fp)
       it != embedded_edges.end(); it++)
     fprintf(fp, "Line {%d} In Surface {%d};\n", (*it)->tag(), tag());
 
-  for(std::set<GVertex *>::iterator it = embedded_vertices.begin();
-      it != embedded_vertices.end(); it++)
+  for(std::set<GVertex *, GEntityPtrLessThan>::iterator it =
+        embedded_vertices.begin(); it != embedded_vertices.end(); it++)
     fprintf(fp, "Point {%d} In Surface {%d};\n", (*it)->tag(), tag());
 
   if(meshAttributes.method == MESH_TRANSFINITE) {

Geo/GRegion.cpp

@@ -529,6 +529,41 @@ double GRegion::computeSolidProperties(std::vector<double> cg,
   return volume;
 }
 
+std::vector<MVertex *> GRegion::getEmbeddedMeshVertices() const
+{
+  std::set<MVertex *> tmp;
+  for(std::vector<GFace *>::const_iterator it = embedded_faces.begin();
+      it != embedded_faces.end(); it++) {
+    tmp.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
+    std::vector<GEdge*> ed = (*it)->edges();
+    for(std::vector<GEdge *>::const_iterator it2 = ed.begin();
+        it2 != ed.end(); it2++) {
+      tmp.insert((*it2)->mesh_vertices.begin(), (*it2)->mesh_vertices.end());
+      if((*it2)->getBeginVertex())
+        tmp.insert((*it2)->getBeginVertex()->mesh_vertices.begin(),
+                   (*it2)->getBeginVertex()->mesh_vertices.end());
+      if((*it2)->getEndVertex())
+        tmp.insert((*it2)->getEndVertex()->mesh_vertices.begin(),
+                   (*it2)->getEndVertex()->mesh_vertices.end());
+    }
+  }
+  for(std::vector<GEdge *>::const_iterator it = embedded_edges.begin();
+      it != embedded_edges.end(); it++) {
+    tmp.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
+    if((*it)->getBeginVertex())
+      tmp.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
+                 (*it)->getBeginVertex()->mesh_vertices.end());
+    if((*it)->getEndVertex())
+      tmp.insert((*it)->getEndVertex()->mesh_vertices.begin(),
+                 (*it)->getEndVertex()->mesh_vertices.end());
+  }
+  for(std::vector<GVertex *>::const_iterator it = embedded_vertices.begin();
+      it != embedded_vertices.end(); it++) {
+    tmp.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
+  }
+  return std::vector<MVertex *>(tmp.begin(), tmp.end());
+}
+
 std::vector<GVertex *> GRegion::vertices() const
 {
   std::set<GVertex *> v;

Geo/GRegion.h

@@ -76,6 +76,7 @@ public:
   std::vector<GVertex *> &embeddedVertices() { return embedded_vertices; }
   std::vector<GEdge *> &embeddedEdges() { return embedded_edges; }
   std::vector<GFace *> &embeddedFaces() { return embedded_faces; }
+  std::vector<MVertex *> getEmbeddedMeshVertices() const;
 
   // edges that bound the region
   virtual std::vector<GEdge *> const &edges() const;

Mesh/meshGFaceExtruded.cpp

@@ -100,7 +100,8 @@ extrudeMesh(GEdge *from, GFace *to, MVertexRTree &pos,
           double x = v->x(), y = v->y(), z = v->z();
           ep->Extrude(j, k + 1, x, y, z);
           if(j != ep->mesh.NbLayer - 1 || k != ep->mesh.NbElmLayer[j] - 1) {
-            MVertex *newv = 0;
+            MVertex *newv = pos.find(x, y, z);
+            if(!newv) {
               if(to->geomType() != GEntity::DiscreteSurface &&
                  to->geomType() != GEntity::BoundaryLayerSurface) {
                 // This can be inefficient, and sometimes useless. We could add an
@@ -113,6 +114,7 @@ extrudeMesh(GEdge *from, GFace *to, MVertexRTree &pos,
                 newv = new MVertex(x, y, z, to);
               }
               to->mesh_vertices.push_back(newv);
+            }
             pos.insert(newv);
             extruded_vertices.push_back(newv);
           }
@@ -190,7 +192,8 @@ static void copyMesh(GFace *from, GFace *to, MVertexRTree &pos)
     double x = v->x(), y = v->y(), z = v->z();
     ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
                 x, y, z);
-    MVertex *newv = 0;
+    MVertex *newv = pos.find(x, y, z);
+    if(!newv) {
       if(to->geomType() != GEntity::DiscreteSurface &&
          to->geomType() != GEntity::BoundaryLayerSurface) {
         SPoint3 xyz(x, y, z);
@@ -203,6 +206,7 @@ static void copyMesh(GFace *from, GFace *to, MVertexRTree &pos)
       to->mesh_vertices.push_back(newv);
       pos.insert(newv);
     }
+  }
 
 #if defined(HAVE_QUADTRI)
   // if performing QuadToTri mesh, cannot simply copy the mesh from the source.
@@ -298,24 +302,20 @@ int MeshExtrudedSurface(
 
   Msg::Info("Meshing surface %d (Extruded)", gf->tag());
 
-  // build an rtree with all the vertices on the boundary of the face gf
+  // build an rtree with all the vertices on the face and its boundary
   MVertexRTree pos(CTX::instance()->geom.tolerance * CTX::instance()->lc);
+  pos.insert(gf->mesh_vertices);
   std::vector<GEdge *> const &edges = gf->edges();
-  std::vector<GEdge *>::const_iterator it = edges.begin();
-  while(it != edges.end()) {
+  for(std::vector<GEdge *>::const_iterator it = edges.begin();
+      it != edges.end(); it++) {
     pos.insert((*it)->mesh_vertices);
     if((*it)->getBeginVertex())
       pos.insert((*it)->getBeginVertex()->mesh_vertices);
     if((*it)->getEndVertex())
       pos.insert((*it)->getEndVertex()->mesh_vertices);
-    ++it;
-  }
-
-  // if the edges of the mesh are constrained, the vertices already
-  // exist on the face--so we add them to the set
-  if(constrainedEdges) {
-    pos.insert(gf->mesh_vertices);
   }
+  std::vector<MVertex*> embedded = gf->getEmbeddedMeshVertices();
+  pos.insert(embedded);
 
   if(ep->geo.Mode == EXTRUDED_ENTITY) {
     // surface is extruded from a curve

Mesh/meshGRegionExtruded.cpp

@@ -174,6 +174,7 @@ static void extrudeMesh(GFace *from, GRegion *to, MVertexRTree &pos)
         double x = v->x(), y = v->y(), z = v->z();
         ep->Extrude(j, k + 1, x, y, z);
         if(j != ep->mesh.NbLayer - 1 || k != ep->mesh.NbElmLayer[j] - 1) {
+          if(!pos.find(x, y, z)) {
             MVertex *newv = new MVertex(x, y, z, to);
             to->mesh_vertices.push_back(newv);
             pos.insert(newv);
@@ -181,6 +182,7 @@ static void extrudeMesh(GFace *from, GRegion *to, MVertexRTree &pos)
         }
       }
     }
+  }
 
 #if defined(HAVE_QUADTRI)
   if(ep && ep->mesh.ExtrudeMesh && ep->mesh.QuadToTri && ep->mesh.Recombine) {
@@ -223,23 +225,23 @@ static void extrudeMesh(GFace *from, GRegion *to, MVertexRTree &pos)
 static void insertAllVertices(GRegion *gr, MVertexRTree &pos)
 {
   pos.insert(gr->mesh_vertices);
+  std::vector<MVertex*> embedded = gr->getEmbeddedMeshVertices();
+  pos.insert(embedded);
   std::vector<GFace *> faces = gr->faces();
-  std::vector<GFace *>::iterator itf = faces.begin();
-  while(itf != faces.end()) {
+  for(std::vector<GFace *>::iterator itf = faces.begin(); itf != faces.end();
+      itf++) {
     pos.insert((*itf)->mesh_vertices);
     std::vector<MVertex *> embedded = (*itf)->getEmbeddedMeshVertices();
     pos.insert(embedded);
     std::vector<GEdge *> const &edges = (*itf)->edges();
-    std::vector<GEdge *>::const_iterator ite = edges.begin();
-    while(ite != edges.end()) {
+    for(std::vector<GEdge *>::const_iterator ite = edges.begin(); ite != edges.end();
+        ite++) {
       pos.insert((*ite)->mesh_vertices);
       if((*ite)->getBeginVertex())
         pos.insert((*ite)->getBeginVertex()->mesh_vertices);
       if((*ite)->getEndVertex())
         pos.insert((*ite)->getEndVertex()->mesh_vertices);
-      ++ite;
     }
-    ++itf;
   }
 }