Merge branch 'enhance-crack-plugin' into 'master'

enhance crack plugin

See merge request !468

CHANGELOG.txt

 4.10.2 (Work-in-progress): fixed regression introduced in 4.9 for recombined
 meshes with boundary layers; new HealShapes command in .geo files; simplified
-calculation of OCC STL bounding boxes; small bug fixes.
+calculation of OCC STL bounding boxes; generalized Crack plugin; small bug fixes.
 
 4.10.1 (May 1, 2022): small bug fixes.
 

examples/api/crack.py

@@ -31,8 +31,14 @@ gmsh.model.addPhysicalGroup(1, new_lines, 101)
 gmsh.model.mesh.generate(2)
 
 gmsh.plugin.setNumber("Crack", "PhysicalGroup", 101)
+# gmsh.plugin.setNumber("Crack", "DebugView", 1)
 gmsh.plugin.run("Crack")
 
+# save all the elements in the mesh (even those that do not belong to any
+# physical group):
+gmsh.option.setNumber("Mesh.SaveAll", 1)
+gmsh.write("crack.msh")
+
 if '-nopopup' not in sys.argv:
     gmsh.fltk.run()
 

examples/api/crack3d.py

@@ -28,6 +28,7 @@ gmsh.model.mesh.generate(3)
 # "crack" the mesh by duplicating the elements and nodes on the small surface
 gmsh.plugin.setNumber("Crack", "Dimension", 2)
 gmsh.plugin.setNumber("Crack", "PhysicalGroup", phys)
+#gmsh.plugin.setNumber("Crack", "DebugView", 1)
 gmsh.plugin.run("Crack")
 
 # save all the elements in the mesh (even those that do not belong to any

src/plugin/Crack.cpp

@@ -21,7 +21,8 @@ StringXNumber CrackOptions_Number[] = {
   {GMSH_FULLRC, "NormalX", nullptr, 0.},
   {GMSH_FULLRC, "NormalY", nullptr, 0.},
   {GMSH_FULLRC, "NormalZ", nullptr, 1.},
-  {GMSH_FULLRC, "NewPhysicalGroup", nullptr, 0}
+  {GMSH_FULLRC, "NewPhysicalGroup", nullptr, 0},
+  {GMSH_FULLRC, "DebugView", nullptr, 0}
 };
 
 extern "C" {
@@ -86,6 +87,7 @@ PView *GMSH_CrackPlugin::execute(PView *view)
   SVector3 normal1d(CrackOptions_Number[3].def, CrackOptions_Number[4].def,
                     CrackOptions_Number[5].def);
   int newPhysical = (int)CrackOptions_Number[6].def;
+  int debug = (int)CrackOptions_Number[7].def;
 
   if(dim != 1 && dim != 2) {
     Msg::Error("Crack dimension should be 1 or 2");
@@ -122,13 +124,19 @@ PView *GMSH_CrackPlugin::execute(PView *view)
     for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++)
       crackElements.push_back(entities[i]->getMeshElement(j));
 
-  // get internal crack nodes and boundary nodes
-  std::set<MVertex *> crackVertices, bndVertices;
+  // get internal crack nodes (and list of connected crack elements), as well as
+  // and boundary nodes
+  std::map<MVertex *, std::vector<MElement *> > crackVertices;
+  std::set<MVertex *> bndVertices;
   if(dim == 1) {
     for(std::size_t i = 0; i < crackElements.size(); i++) {
       for(std::size_t j = 0; j < crackElements[i]->getNumVertices(); j++) {
         MVertex *v = crackElements[i]->getVertex(j);
-        crackVertices.insert(v);
+        auto it = crackVertices.find(v);
+        if(it == crackVertices.end())
+          crackVertices[v] = {crackElements[i]};
+        else
+          it->second.push_back(crackElements[i]);
       }
       for(std::size_t j = 0; j < crackElements[i]->getNumPrimaryVertices();
           j++) {
@@ -145,7 +153,11 @@ PView *GMSH_CrackPlugin::execute(PView *view)
     for(std::size_t i = 0; i < crackElements.size(); i++) {
       for(std::size_t j = 0; j < crackElements[i]->getNumVertices(); j++) {
         MVertex *v = crackElements[i]->getVertex(j);
-        crackVertices.insert(v);
+        auto it = crackVertices.find(v);
+        if(it == crackVertices.end())
+          crackVertices[v] = {crackElements[i]};
+        else
+          it->second.push_back(crackElements[i]);
       }
       for(int j = 0; j < crackElements[i]->getNumEdges(); j++) {
         EdgeData ed(crackElements[i]->getEdge(j));
@@ -199,8 +211,11 @@ PView *GMSH_CrackPlugin::execute(PView *view)
   for(auto it = bndVertices.begin(); it != bndVertices.end(); it++)
     crackVertices.erase(*it);
 
-  // compute elements on one side of the crack
-  std::set<MElement *> oneside;
+  // compute elements on the positive side of the crack, and keep track of each
+  // node in the element that leads to categorizing the element on this side
+  // (this allows to handle the degenerate case where the same element touches
+  // the same crack on both sides, with different nodes - cf. #1750)
+  std::map<MElement *, std::vector<std::size_t> > oneside;
   std::vector<GEntity *> allentities;
   m->getEntities(allentities);
   for(std::size_t ent = 0; ent < allentities.size(); ent++) {
@@ -208,40 +223,43 @@ PView *GMSH_CrackPlugin::execute(PView *view)
     for(std::size_t i = 0; i < allentities[ent]->getNumMeshElements(); i++) {
       MElement *e = allentities[ent]->getMeshElement(i);
       for(std::size_t j = 0; j < e->getNumVertices(); j++) {
-        if(crackVertices.find(e->getVertex(j)) != crackVertices.end()) {
-          // element touches the crack: find the closest crack element
-          SPoint3 b = e->barycenter();
-          double d = 1e200;
-          MElement *ce = nullptr;
-          for(std::size_t k = 0; k < crackElements.size(); k++) {
-            double d2 = b.distance(crackElements[k]->barycenter());
-            if(d2 < d) {
-              d = d2;
-              ce = crackElements[k];
-            }
-          }
-          SVector3 dv = SVector3(e->barycenter(), ce->barycenter());
+        auto it = crackVertices.find(e->getVertex(j));
+        if(it == crackVertices.end()) continue;
+        // the element touches the crack by at least one node: determine if the
+        // vector joining its barycenter to the barycenter of one of the
+        // connected crack elements is not in the same direction as the normal
+        // to the crack element; if the condition is fulfilled for one of the
+        // connected crack elements, we consider the element lies on the
+        // "positive side" of the crack
+        SPoint3 b = e->barycenter();
+        for(auto ce : it->second) {
+          SVector3 dv = SVector3(b, ce->barycenter());
           SVector3 n;
           if(dim == 1)
             n = crossprod(normal1d, ce->getEdge(0).tangent());
           else
             n = ce->getFace(0).normal();
-          if(dot(n, dv) < 0) { oneside.insert(e); }
+          if(dot(n, dv) < 0) {
+            auto it2 = oneside.find(e);
+            if(it2 == oneside.end())
+              oneside[e] = {j};
+            else {
+              if(std::find(it2->second.begin(), it2->second.end(), j) ==
+                 it2->second.end())
+                it2->second.push_back(j);
+            }
+          }
         }
       }
     }
   }
 
-  /*
-  FILE *fp = fopen("debug.pos", "w");
-  if(fp){
-    fprintf(fp, "View \"Ele < 0\" {\n");
-    for(auto it = oneside.begin(); it != oneside.end(); it++)
-      (*it)->writePOS(fp, false, true, false, false, false, false);
-    fprintf(fp, "};\n");
-    fclose(fp);
+  if(debug) {
+    std::map<int, std::vector<double> > d;
+    for(auto e : oneside) d[e.first->getNum()] = {(double)e.first->getNum()};
+    view = new PView("Elements on positive side of crack", "ElementData",
+                     GModel::current(), d);
   }
-  */
 
   // create new crack entity
 
@@ -278,7 +296,7 @@ PView *GMSH_CrackPlugin::execute(PView *view)
   // duplicate internal crack nodes
   std::map<MVertex *, MVertex *> vxv;
   for(auto it = crackVertices.begin(); it != crackVertices.end(); it++) {
-    MVertex *v = *it;
+    MVertex *v = it->first;
     MVertex *newv = new MVertex(v->x(), v->y(), v->z(), crackEntity);
     crackEntity->mesh_vertices.push_back(newv);
     vxv[v] = newv;
@@ -304,14 +322,16 @@ PView *GMSH_CrackPlugin::execute(PView *view)
 
   // replace vertices in elements on one side of the crack
   for(auto it = oneside.begin(); it != oneside.end(); it++) {
-    MElement *e = *it;
-    for(std::size_t i = 0; i < e->getNumVertices(); i++) {
-      if(vxv.count(e->getVertex(i))) e->setVertex(i, vxv[e->getVertex(i)]);
+    MElement *e = it->first;
+    for(auto i : it->second) {
+      MVertex *v = e->getVertex(i);
+      if(vxv.count(v))
+        e->setVertex(i, vxv[v]);
+      else
+        Msg::Warning("Mesh node %lu not found in cracked nodes", v->getNum());
     }
   }
 
-  // m->pruneMeshVertexAssociations();
-
   CTX::instance()->mesh.changed = ENT_ALL;
 
   return view;