diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 6028b116cf6b8fcb33c9485060547171de2cce90..8236a46c119f86f28d1a08f4394c49e68acbc747 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -5,6 +5,7 @@
#include <stdlib.h>
#include <sstream>
+#include <stack>
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "GModel.h"
@@ -1263,10 +1264,39 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
return num;
}
+
static void recurConnectMElementsByMFace(const MFace &f,
std::multimap<MFace, MElement*, Less_Face> &e2f,
std::set<MElement*> &group,
- std::set<MFace, Less_Face> &touched)
+ std::set<MFace, Less_Face> &touched,
+ int recur_level)
+{
+ // this is very slow...
+ std::stack<MFace> _stack;
+ _stack.push(f);
+
+ while(!_stack.empty()){
+ MFace ff = _stack.top();
+ _stack.pop();
+ if (touched.find(ff) == touched.end()){
+ touched.insert(ff);
+ for (std::multimap<MFace, MElement*, Less_Face>::iterator it = e2f.lower_bound(ff);
+ it != e2f.upper_bound(ff); ++it){
+ group.insert(it->second);
+ for (int i = 0; i < it->second->getNumFaces(); ++i){
+ _stack.push(it->second->getFace(i));
+ }
+ }
+ }
+ }
+ printf("group pf %d elements found\n",group.size());
+}
+
+static void recurConnectMElementsByMFaceOld(const MFace &f,
+ std::multimap<MFace, MElement*, Less_Face> &e2f,
+ std::set<MElement*> &group,
+ std::set<MFace, Less_Face> &touched,
+ int recur_level)
{
if (touched.find(f) != touched.end()) return;
touched.insert(f);
@@ -1274,7 +1304,7 @@ static void recurConnectMElementsByMFace(const MFace &f,
it != e2f.upper_bound(f); ++it){
group.insert(it->second);
for (int i = 0; i < it->second->getNumFaces(); ++i){
- recurConnectMElementsByMFace(it->second->getFace(i), e2f, group, touched);
+ recurConnectMElementsByMFace(it->second->getFace(i), e2f, group, touched, recur_level+1);
}
}
}
@@ -1291,7 +1321,7 @@ static int connectedVolumes(std::vector<MElement*> &elements,
while(!e2f.empty()){
std::set<MElement*> group;
std::set<MFace, Less_Face> touched;
- recurConnectMElementsByMFace(e2f.begin()->first, e2f, group, touched);
+ recurConnectMElementsByMFace(e2f.begin()->first, e2f, group, touched, 0);
std::vector<MElement*> temp;
temp.insert(temp.begin(), group.begin(), group.end());
regs.push_back(temp);
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index 76f9ec103000f13a3ff0d0f7d853e8cc675dda93..2a0d44dd6bc91ddd12122c3db7d705268c47e4e7 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/Geo/GeomMeshMatcher.cpp
@@ -618,6 +618,10 @@ static void copy_vertices (GRegion *to, GRegion *from, std::map<MVertex*,MVertex
}
static void copy_vertices (GEdge *to, GEdge *from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
to->deleteMesh();
+ if (!from){
+ Msg::Warning("Edge %d in the mesh do not match any edge of the model",to->tag());
+ return;
+ }
for (int i=0;i<from->mesh_vertices.size();i++){
MVertex *v_from = from->mesh_vertices[i];
double t;
@@ -639,8 +643,10 @@ static void copy_vertices (GFace *geom, GFace *mesh, std::map<MVertex*,MVertex*>
double DDD = ( v_from->x() - gp.x()) * ( v_from->x() - gp.x()) +
( v_from->y() - gp.y()) * ( v_from->y() - gp.y()) +
( v_from->z() - gp.z()) * ( v_from->z() - gp.z()) ;
- if (sqrt(DDD) > 1.e-3)Msg::Error("Impossible to match one point Original point %f %f %f New point %f %f %f",
+ if (sqrt(DDD) > 1.e-1)Msg::Error("Impossible to match one point Original point %f %f %f New point %f %f %f",
v_from->x(), v_from->y(), v_from->z(),gp.x(), gp.y(), gp.z());
+ else if (sqrt(DDD) > 1.e-3)Msg::Warning("One mesh vertex %f %f %f of GFace %d \n is difficult to match : closest point %f %f %f",
+ v_from->x(), v_from->y(), v_from->z(),geom->tag(),gp.x(), gp.y(), gp.z());
MFaceVertex *v_to = new MFaceVertex (v_from->x(),v_from->y(),v_from->z(), geom, gp.u(),gp.v() );
geom->mesh_vertices.push_back(v_to);
_mesh_to_geom[v_from] = v_to;
@@ -649,7 +655,9 @@ static void copy_vertices (GFace *geom, GFace *mesh, std::map<MVertex*,MVertex*>
}
template <class ELEMENT>
-static void copy_elements (std::vector<ELEMENT*> &to, std::vector<ELEMENT*> &from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+static void copy_elements (std::vector<ELEMENT*> &to,
+ std::vector<ELEMENT*> &from,
+ std::map<MVertex*,MVertex*> &_mesh_to_geom){
MElementFactory toto;
to.clear();
for (int i=0;i < from.size();i++){
@@ -663,26 +671,35 @@ static void copy_elements (std::vector<ELEMENT*> &to, std::vector<ELEMENT*> &fro
}
-void copy_vertices (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+void copy_vertices (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_mesh_to_geom,
+ std::vector<Pair<GVertex*, GVertex*> > *coresp_v,
+ std::vector<Pair<GEdge*, GEdge*> > *coresp_e,
+ std::vector<Pair<GFace*, GFace*> > *coresp_f){
+
// copy all elements
- for (GModel::viter vit = geom->firstVertex() ; vit != geom->lastVertex(); ++vit)
- copy_vertices(*vit,mesh->getVertexByTag((*vit)->tag()),_mesh_to_geom);
- for (GModel::eiter eit = geom->firstEdge() ; eit != geom->lastEdge(); ++eit)
- copy_vertices(*eit,mesh->getEdgeByTag((*eit)->tag()),_mesh_to_geom);
- for (GModel::fiter fit = geom->firstFace() ; fit != geom->lastFace(); ++fit)
- copy_vertices(*fit,mesh->getFaceByTag((*fit)->tag()),_mesh_to_geom);
+ for (int i=0;i<coresp_v->size();++i)
+ copy_vertices((*coresp_v)[i].first(),(*coresp_v)[i].second(),_mesh_to_geom);
+ for (int i=0;i<coresp_e->size();++i)
+ copy_vertices((*coresp_e)[i].first(),(*coresp_e)[i].second(),_mesh_to_geom);
+ for (int i=0;i<coresp_f->size();++i)
+ copy_vertices((*coresp_f)[i].first(),(*coresp_f)[i].second(),_mesh_to_geom);
for (GModel::riter rit = geom->firstRegion() ; rit != geom->lastRegion(); ++rit)
copy_vertices(*rit,mesh->getRegionByTag((*rit)->tag()),_mesh_to_geom);
}
-void copy_elements (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+void copy_elements (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_mesh_to_geom,
+ std::vector<Pair<GVertex*, GVertex*> > *coresp_v,
+ std::vector<Pair<GEdge*, GEdge*> > *coresp_e,
+ std::vector<Pair<GFace*, GFace*> > *coresp_f){
// copy all elements
- for (GModel::viter vit = geom->firstVertex() ; vit != geom->lastVertex(); ++vit)
- if (mesh->getVertexByTag((*vit)->tag())) copy_elements<MPoint>((*vit)->points,mesh->getVertexByTag((*vit)->tag())->points,_mesh_to_geom);
- for (GModel::eiter eit = geom->firstEdge() ; eit != geom->lastEdge(); ++eit)
- copy_elements<MLine>((*eit)->lines,mesh->getEdgeByTag((*eit)->tag())->lines,_mesh_to_geom);
- for (GModel::fiter fit = geom->firstFace() ; fit != geom->lastFace(); ++fit) {
- copy_elements<MTriangle>((*fit)->triangles,mesh->getFaceByTag((*fit)->tag())->triangles,_mesh_to_geom);
- copy_elements<MQuadrangle>((*fit)->quadrangles,mesh->getFaceByTag((*fit)->tag())->quadrangles,_mesh_to_geom);
+ for (int i=0;i<coresp_v->size();++i)
+ copy_elements<MPoint>((*coresp_v)[i].first()->points,(*coresp_v)[i].second()->points,_mesh_to_geom);
+
+ for (int i=0;i<coresp_e->size();++i)
+ copy_elements<MLine>((*coresp_e)[i].first()->lines,(*coresp_e)[i].second()->lines,_mesh_to_geom);
+
+ for (int i=0;i<coresp_f->size();++i){
+ copy_elements<MTriangle>((*coresp_f)[i].first()->triangles,(*coresp_f)[i].second()->triangles,_mesh_to_geom);
+ copy_elements<MQuadrangle>((*coresp_f)[i].first()->quadrangles,(*coresp_f)[i].second()->quadrangles,_mesh_to_geom);
}
for (GModel::riter rit = geom->firstRegion() ; rit != geom->lastRegion(); ++rit) {
copy_elements<MTetrahedron>((*rit)->tetrahedra,mesh->getRegionByTag((*rit)->tag())->tetrahedra,_mesh_to_geom);
@@ -713,8 +730,8 @@ int GeomMeshMatcher::match(GModel *geom, GModel *mesh)
matchRegions(geom, mesh, coresp_f,ok);
if (ok) {
std::map<MVertex*,MVertex*> _mesh_to_geom;
- copy_vertices(geom, mesh, _mesh_to_geom);
- copy_elements(geom, mesh, _mesh_to_geom);
+ copy_vertices(geom, mesh, _mesh_to_geom,coresp_v,coresp_e,coresp_f);
+ copy_elements(geom, mesh, _mesh_to_geom,coresp_v,coresp_e,coresp_f);
return 1;
} else {
Msg::Error("Could not match every region !");
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 5849d5c0f57a7d394f36f58199c360db41a406ec..f42fdde41d5a012329d3cdd5f123e773e9c2a92f 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -998,6 +998,7 @@ void optimalPointFrontal (GFace *gf,
newPoint[1] = midpoint[1] + d * dir[1];
}
+///
void optimalPointFrontalB (GFace *gf,
MTri3* worst,