diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index bc965d74c91f46268aface68d889c048d131c769..27569c327bce349801f8ce1018e38acbe544850f 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,5 +1,5 @@
4.10.2 (Work-in-progress): fixed regression introduced in 4.9 for recombined
-meshes with boundary layers; small bug fixes.
+meshes with boundary layers; generalized Crack plugin; small bug fixes.
4.10.1 (May 1, 2022): small bug fixes.
diff --git a/examples/api/crack.py b/examples/api/crack.py
index 7a3cab9f4576ee689d183db2e70800a83fdd45ae..f5f72e43f172fb2da22848e4f915fbba368b1baa 100644
--- a/examples/api/crack.py
+++ b/examples/api/crack.py
@@ -33,6 +33,11 @@ gmsh.model.mesh.generate(2)
gmsh.plugin.setNumber("Crack", "PhysicalGroup", 101)
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()
diff --git a/src/plugin/Crack.cpp b/src/plugin/Crack.cpp
index 5b28d36defec2d61f1cd27e4b2f2e15e91c9ab5d..f78eb233b05b4ce004cb3b6e8549eb5cad917b65 100644
--- a/src/plugin/Crack.cpp
+++ b/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, "Debug", 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));
@@ -208,40 +220,46 @@ 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();
+ bool found = false;
+ 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) {
+ oneside.insert(e);
+ found = true;
+ break;
+ }
}
+ if(found) break; // we're done for the element
}
}
}
- /*
- 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) {
+ Msg::Info("Writing 'debug.pos' file with elements detected on one side "
+ "of the crack");
+ FILE *fp = fopen("debug.pos", "w");
+ if(fp){
+ fprintf(fp, "View \"Elements on one side\" {\n");
+ for(auto it = oneside.begin(); it != oneside.end(); it++)
+ (*it)->writePOS(fp, false, true, false, false, false, false);
+ fprintf(fp, "};\n");
+ fclose(fp);
+ }
}
- */
// 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;