diff --git a/Plugin/Crack.cpp b/Plugin/Crack.cpp index 403087080ec3599962e76e66d59f15f4ef9dd018..2387e1eee8b44653f76c5a5f69f4cea8179086ec 100644 --- a/Plugin/Crack.cpp +++ b/Plugin/Crack.cpp @@ -169,64 +169,56 @@ PView *GMSH_CrackPlugin::execute(PView *view) } } } - for(std::set<MVertex*>::iterator it = bndVertices.begin(); it != bndVertices.end(); it++) crackVertices.erase(*it); - // compute smoothed normals on crack vertices - std::map<MVertex*, std::vector<MElement*> > vxe; - for(unsigned int i = 0; i < crackElements.size(); i++){ - MElement *e = crackElements[i]; - for(int k = 0; k < e->getNumVertices(); k++) - vxe[e->getVertex(k)].push_back(e); - } - std::map<MVertex*, SVector3> vxn; - for(std::map<MVertex*, std::vector<MElement*> >::iterator it = vxe.begin(); - it != vxe.end(); it++){ - SVector3 n; - for(unsigned int i = 0; i < it->second.size(); i++){ - if(dim == 1) - n += crossprod(normal1d, it->second[i]->getEdge(0).tangent()); - else - n += it->second[i]->getFace(0).normal(); - } - n.normalize(); - vxn[it->first] = n; - } - // compute elements on one side of the crack - vxe.clear(); + std::set<MElement*> oneside; std::vector<GEntity*> allentities; m->getEntities(allentities); for(unsigned int ent = 0; ent < allentities.size(); ent++){ - //if(allentities[ent]->dim() != dim + 1) continue; if(crackEntities.find(allentities[ent]) != crackEntities.end()) continue; for(unsigned int i = 0; i < allentities[ent]->getNumMeshElements(); i++){ MElement *e = allentities[ent]->getMeshElement(i); for(int j = 0; j < e->getNumVertices(); j++){ - MVertex *v = e->getVertex(j); - if(crackVertices.find(v) != crackVertices.end()){ - MVertex *vclose = 0; - double d = 1e22; - SVector3 dv; - for(std::set<MVertex*>::iterator it = crackVertices.begin(); - it != crackVertices.end(); it++){ - MVertex *v = *it; - double d2 = v->point().distance(e->barycenter()); + 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 = 0; + for(unsigned int k = 0; k < crackElements.size(); k++){ + double d2 = b.distance(crackElements[k]->barycenter()); if(d2 < d){ d = d2; - vclose = v; - dv = SVector3(e->barycenter(), vclose->point()); + ce = crackElements[k]; } } - if(dot(vxn[vclose], dv) < 0) - vxe[v].push_back(e); + SVector3 dv = SVector3(e->barycenter(), 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); + } } } } } + /* + FILE *fp = fopen("debug.pos", "w"); + if(fp){ + fprintf(fp, "View \"Ele < 0\" {\n"); + for(std::set<MElement*>::iterator 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 GEdge *crackEdge = 0; GFace *crackFace = 0; @@ -271,13 +263,7 @@ PView *GMSH_CrackPlugin::execute(PView *view) } // replace vertices in elements on one side of the crack - std::set<MElement*> eles; - for(std::map<MVertex*, std::vector<MElement*> >::iterator it = vxe.begin(); - it != vxe.end(); it++){ - for(unsigned int i = 0; i < it->second.size(); i++) - eles.insert(it->second[i]); - } - for(std::set<MElement*>::iterator it = eles.begin(); it != eles.end(); it++){ + for(std::set<MElement*>::iterator it = oneside.begin(); it != oneside.end(); it++){ MElement *e = *it; for(int i = 0; i < e->getNumVertices(); i++){ if(vxv.count(e->getVertex(i)))