From 7cd86dbabba7fdd20211554f5c0cdc336fd65b90 Mon Sep 17 00:00:00 2001
From: Bastien Gorissen <bastien.gorissen@cenaero.be>
Date: Wed, 8 Feb 2012 16:33:29 +0000
Subject: [PATCH] Modifications to the matcher to make it work across seams and
with closed curves.
---
Geo/GeomMeshMatcher.cpp | 294 +++++++++++++++++++---------------------
1 file changed, 136 insertions(+), 158 deletions(-)
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index 2a0d44dd6b..6ea652e30c 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/Geo/GeomMeshMatcher.cpp
@@ -67,98 +67,62 @@ template <class T> T findMatching(std::vector<Pair<T,T> >& matching, T& entity)
// Private
+
+// ------------------------------------------------------------[ Matching vertices ]
+
std::vector<Pair<GVertex*,GVertex*> >*
GeomMeshMatcher::matchVertices(GModel* m1, GModel *m2, bool& ok)
{
// Vector that will be returned.
std::vector<Pair<GVertex*,GVertex*> >* coresp_v = new std::vector<Pair<GVertex*,
GVertex*> >;
-
- // Building two vectors with the models' entities.
- std::vector<GEntity*> m1_entities;
- m1->getEntities(m1_entities);
- std::vector<GEntity*> m2_entities;
- m2->getEntities(m2_entities);
-
int num_matched_vertices = 0;
- int num_total_vertices = 0;
- int num_mesh_vertices = 0;
-
- int counter1 = 0;
+ int num_total_vertices = m1->getNumVertices();
std::vector<GVertex*> vertices;
- for (std::vector<GEntity*>::iterator entity1 = m1_entities.begin();
- entity1 != m1_entities.end();
- entity1++)
- {
- // If the entity isn't a vertex, don't bother.
- if ((*entity1)->dim() != 0) continue;
- num_total_vertices++;
- GVertex* v1 = (GVertex*) *entity1;
+ for(GModel::viter vit = m1->firstVertex(); vit != m1->lastVertex(); vit++) {
+
+ GVertex* v1 = (GVertex*) *vit;
// FIXME: need a *much* better way to fix the tolerance...
double tol = 10e-8;
- int counter2 = 0;
discreteVertex* best_candidate;
GEntity* best_candidate_ge = 0;
double best_score = DBL_MAX;
- for (std::vector<GEntity*>::iterator entity2 = m2_entities.begin();
- entity2 != m2_entities.end();
- entity2++)
- {
- if ((*entity2)->dim() != 0) continue;
- num_mesh_vertices++;
- for (unsigned int ed = 0;
- ed < ((discreteVertex*) *entity2)->getNumMeshElements();
- ed++)
- {
- discreteVertex* v2 = (discreteVertex*) *entity2;
-
- // We match the vertices if their coordinates are the same under the
- // specified tolerance.
- double score = std::max(fabs(v1->x() - v2->x()),
- std::max(fabs(v1->y() - v2->y()),
- fabs(v1->z() - v2->z())));
- if (score < tol && score < best_score) {
- best_candidate = v2;
- best_candidate_ge = (*entity2);
- best_score = score;
- }
+ for(GModel::viter vit2 = m2->firstVertex(); vit2 != m2->lastVertex(); vit2++) {
+
+ discreteVertex* v2 = (discreteVertex*) *vit2;
+
+ // We match the vertices if their coordinates are the same under the
+ // specified tolerance.
+ double score = std::max(fabs(v1->x() - v2->x()),
+ std::max(fabs(v1->y() - v2->y()),
+ fabs(v1->z() - v2->z())));
+ if (score < tol && score < best_score) {
+ best_candidate = v2;
+ best_candidate_ge = (*vit2);
+ best_score = score;
}
- counter2++;
}
if (best_score != DBL_MAX) {
- Msg::Info("Model Vertex %i (geom) and %i (mesh) match",
- (*entity1)->tag(),
+ Msg::Debug("Model Vertex %i (geom) and %i (mesh) match",
+ v1->tag(),
best_candidate_ge->tag());
- coresp_v->push_back(Pair<GVertex*,GVertex*>((GVertex*) *entity1,
- (GVertex*) best_candidate_ge));
- ((GVertex*) best_candidate_ge)->setTag(((GVertex*) *entity1)->tag());
- //m2->remove((GVertex*) best_candidate_ge);
- //vertices.push_back((GVertex*) best_candidate_ge);
- for (unsigned int v = 0; v < ((GVertex*) best_candidate_ge)->getNumMeshVertices(); v++) {
- ((GVertex*) best_candidate_ge)->getMeshVertex(v)->setEntity((GVertex*) *entity1);
- }
+ coresp_v->push_back(Pair<GVertex*,GVertex*>(v1, best_candidate));
num_matched_vertices++;
}
- counter1++;
}
- //for (std::vector<GVertex*>::iterator vert = vertices.begin();
- //vert != vertices.end();
- //vert++)
- //m2->add(*vert);
-
- Msg::Info("Vertices matched : %i / %i (%i)", num_matched_vertices, num_total_vertices, num_mesh_vertices);
- //if(num_matched_vertices != num_total_vertices) ok = false;
+ Msg::Info("Matched %i vertices out of %i.", num_matched_vertices, num_total_vertices);
return (coresp_v);
}
+// ------------------------------------------------------------[ Matching edges ]
std::vector<Pair<GEdge*,GEdge*> >*
GeomMeshMatcher::matchEdges(GModel* m1, GModel* m2,
@@ -171,27 +135,49 @@ GeomMeshMatcher::matchEdges(GModel* m1, GModel* m2,
// Vector that will be returned.
std::vector<Pair<GEdge*,GEdge*> >* coresp_e = new std::vector<Pair<GEdge*,GEdge*> >;
- // Building two vectors with the models' entities.
- std::vector<GEntity*> m1_entities;
- m1->getEntities(m1_entities);
- std::vector<GEntity*> m2_entities;
- m2->getEntities(m2_entities);
+ std::vector<GEdge*> closed_curves;
+
+ for(GModel::eiter eit = m1->firstEdge(); eit != m1->lastEdge(); eit++) {
- for (std::vector<GEntity*>::iterator entity1 = m1_entities.begin();
- entity1 != m1_entities.end();
- entity1++)
- {
- if ((*entity1)->dim() != 1) continue;
num_total_edges++;
- GVertex* v1 = ((GEdge*)(*entity1))->getBeginVertex();
- GVertex* v2 = ((GEdge*)(*entity1))->getEndVertex();
+ GEdge* e1 = (GEdge*) *eit;
+
+ GVertex* v1 = e1->getBeginVertex();
+ GVertex* v2 = e1->getEndVertex();
std::vector<GEdge*> common_edges;
std::vector<std::list<GEdge*> > lists;
- lists.push_back((findMatching<GVertex*>(*coresp_v,v1))->edges());
- lists.push_back((findMatching<GVertex*>(*coresp_v,v2))->edges());
- getIntersection<GEdge*>(common_edges, lists);
+
+ if (v1 == v2) {
+ Msg::Debug("Found a closed curve");
+ closed_curves.push_back(e1);
+
+ for (GModel::eiter eit2 = m2->firstEdge(); eit2 != m2->lastEdge(); eit2++) {
+ GEdge* e2 = (GEdge*)*eit2;
+ GVertex* v3 = e2->getBeginVertex();
+ GVertex* v4 = e2->getEndVertex();
+ if (v3 == v4) {
+ Msg::Debug("Found a loop (%i) in the mesh %i %i", e2->tag(), v3->tag(), v3->tag());
+ common_edges.push_back(e2);
+ }
+
+ }
+ } else {
+ //if (coresp_v->count(vfindMatching<GVertex*>(*coresp_v,v1)1) > 0 && coresp_v->count(v2) > 0) {
+ bool ok1 = false;
+ bool ok2 = false;
+ if (findMatching<GVertex*>(*coresp_v,v1) != 0) {
+ ok1 = true;
+ lists.push_back((findMatching<GVertex*>(*coresp_v,v1))->edges());
+ }
+ if (findMatching<GVertex*>(*coresp_v,v2) != 0) {
+ ok2 = true;
+ lists.push_back((findMatching<GVertex*>(*coresp_v,v2))->edges());
+ }
+ if (ok1 && ok2)
+ getIntersection<GEdge*>(common_edges, lists);
+ }
GEdge* choice = 0;
if (common_edges.size() == 0) continue;
@@ -202,7 +188,7 @@ GeomMeshMatcher::matchEdges(GModel* m1, GModel* m2,
// us use those bounding boxes !
// So, first step is to build an array of points taken on the geo entity
// Then, compute the minimal bounding box
- SOrientedBoundingBox geo_obb = ((GEdge*)(*entity1))->getOBB();
+ SOrientedBoundingBox geo_obb = e1->getOBB();
double best_score = DBL_MAX;
// Next, let's iterate over the mesh entities.
@@ -217,34 +203,22 @@ GeomMeshMatcher::matchEdges(GModel* m1, GModel* m2,
}
}
}
- Msg::Info("Edges %i (geom) and %i (mesh) match.",
- ((GEdge*)*entity1)->tag(),
+ Msg::Debug("Edges %i (geom) and %i (mesh) match.",
+ e1->tag(),
choice->tag());
- coresp_e->push_back(Pair<GEdge*,GEdge*>((GEdge*) *entity1 ,
- choice));
- choice->setTag(((GEdge*) *entity1)->tag());
+ coresp_e->push_back(Pair<GEdge*,GEdge*>(e1, choice));
+ //choice->setTag(e1->tag());
- // This reverses the edge if it's not parametrized in the right direction
-/*
- if (choice->getBeginVertex() == findMatching<GVertex*>(*coresp_v,v2) &&
- choice->getEndVertex() == findMatching<GVertex*>(*coresp_v,v1)) {
- Msg::Info("Wrong parametrization direction, reversing.");
- choice->reverse();
- }
-*/
- /*for (unsigned int v = 0; v < ((GEdge*) choice)->getNumMeshVertices(); v++) {
- if (((GEdge*) choice)->getMeshVertex(v)->onWhat()->dim() > 0)
- ((GEdge*) choice)->getMeshVertex(v)->setEntity((GEdge*) *entity1);
- }*/
num_matched_edges++;
}
- Msg::Info("Edges matched : %i / %i", num_matched_edges, num_total_edges);
+ Msg::Info("Matched %i edges out of %i.", num_matched_edges, num_total_edges);
if(num_matched_edges != num_total_edges) ok = false;
return (coresp_e);
}
-//-----------------------------------------------------------------------------
+// ------------------------------------------------------------[ Matching faces ]
+
std::vector<Pair<GFace*,GFace*> >*
GeomMeshMatcher:: matchFaces(GModel* m1, GModel* m2,
@@ -255,23 +229,19 @@ GeomMeshMatcher:: matchFaces(GModel* m1, GModel* m2,
std::vector<Pair<GFace*,GFace*> >* coresp_f = new std::vector<Pair<GFace*,GFace*> >;
- std::vector<GEntity*> m1_entities;
- m1->getEntities(m1_entities);
- std::vector<GEntity*> m2_entities;
- m2->getEntities(m2_entities);
+ for(GModel::fiter fit = m1->firstFace(); fit != m1->lastFace(); fit++) {
- for (std::vector<GEntity*>::iterator entity1 = m1_entities.begin();
- entity1 != m1_entities.end();
- entity1++)
- {
- if ((*entity1)->dim() != 2) continue;
+ GFace* f1 = (GFace*) *fit;
num_total_faces++;
-
+
std::vector<std::list<GFace*> > lists;
- std::list<GEdge*> boundary_edges = ((GEdge*)(*entity1))->edges();
+
+ std::list<GEdge*> boundary_edges = f1->edges();
for (std::list<GEdge*>::iterator boundary_edge = boundary_edges.begin();
boundary_edge != boundary_edges.end(); boundary_edge++) {
- lists.push_back(findMatching<GEdge*>(*coresp_e,*boundary_edge)->faces());
+ // if (boundary_edge->getBeginVertex() == boundary_edge->getEndVertex() &&
+ if (!(*boundary_edge)->isSeam(f1))
+ lists.push_back(findMatching<GEdge*>(*coresp_e,*boundary_edge)->faces());
}
std::vector<GFace*> common_faces;
getIntersection<GFace*>(common_faces, lists);
@@ -281,9 +251,8 @@ GeomMeshMatcher:: matchFaces(GModel* m1, GModel* m2,
} else if (common_faces.size() > 1) {
// Then, compute the minimal bounding box
- SOrientedBoundingBox geo_obb = ((GFace*)(*entity1))->getOBB();
+ SOrientedBoundingBox geo_obb = f1->getOBB();
- //GFace* choice = 0;
double best_score = DBL_MAX;
// Next, let's iterate over the mesh entities.
for (std::vector<GFace*>::iterator candidate = common_faces.begin();
@@ -299,24 +268,24 @@ GeomMeshMatcher:: matchFaces(GModel* m1, GModel* m2,
}
}
- coresp_f->push_back(Pair<GFace*,GFace*>((GFace*) *entity1 ,
- choice));
- choice->setTag(((GFace*) *entity1)->tag());
- /*for (unsigned int v = 0; v < ((GFace*) choice)->getNumMeshVertices(); v++) {
- if(((GFace*) choice)->getMeshVertex(v)->onWhat()->dim() > 1)
- ((GFace*) choice)->getMeshVertex(v)->setEntity((GFace*) *entity1);
- }*/
+ Msg::Debug("Faces %i (geom) and %i (mesh) match.",
+ f1->tag(),
+ choice->tag());
+ coresp_f->push_back(Pair<GFace*,GFace*>(f1 ,
+ choice));
+ choice->setTag(f1->tag());
num_matched_faces++;
}
- Msg::Info("Faces matched : %i / %i", num_matched_faces, num_total_faces);
- if(num_matched_faces != num_total_faces) ok = false;
+ Msg::Info("Matched %i faces out of %i.", num_matched_faces, num_total_faces);
return coresp_f;
}
+// ------------------------------------------------------------[ Matching regions ]
+
std::vector<Pair<GRegion*,GRegion*> >*
GeomMeshMatcher::matchRegions(GModel* m1, GModel* m2,
std::vector<Pair<GFace*,GFace*> >* coresp_f, bool& ok)
@@ -415,6 +384,7 @@ GeomMeshMatcher::matchRegions(GModel* m1, GModel* m2,
}
+
// Public
GeomMeshMatcher* GeomMeshMatcher::instance()
{
@@ -501,7 +471,6 @@ static GPoint getGFace (MVertex *v1, GModel *gm, const double TOL){
return gpBest;
}
-
int GeomMeshMatcher::forceTomatch(GModel *geom, GModel *mesh, const double TOL)
{
// assume that the geometry is the right one
@@ -616,12 +585,26 @@ static void copy_vertices (GRegion *to, GRegion *from, std::map<MVertex*,MVertex
_mesh_to_geom[v_from] = v_to;
}
}
-static void copy_vertices (GEdge *to, GEdge *from, std::map<MVertex*,MVertex*> &_mesh_to_geom){
+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;
}
+ if (!to){
+ Msg::Warning("Edge %d in the geometry do not match any edge of the mesh",from->tag());
+ return;
+ }
+
+ if (from->getBeginVertex() == from->getEndVertex()) {
+ MVertex *v_from = from->getBeginVertex()->mesh_vertices[0];
+ double t;
+ GPoint gp = to->closestPoint(SPoint3(v_from->x(),v_from->y(),v_from->z()), t );
+ MEdgeVertex *v_to = new MEdgeVertex (gp.x(),gp.y(),gp.z(), to, gp.u() );
+ to->mesh_vertices.push_back(v_to);
+ _mesh_to_geom[v_from] = v_to;
+ }
+
for (int i=0;i<from->mesh_vertices.size();i++){
MVertex *v_from = from->mesh_vertices[i];
double t;
@@ -629,6 +612,9 @@ static void copy_vertices (GEdge *to, GEdge *from, std::map<MVertex*,MVertex*> &
MEdgeVertex *v_to = new MEdgeVertex (gp.x(),gp.y(),gp.z(), to, gp.u() );
to->mesh_vertices.push_back(v_to);
_mesh_to_geom[v_from] = v_to;
+ if (v_from->getNum() == 3646) {
+ printf("FOUND IT!!\n");
+ }
}
// printf("Ending Edge %d %d vertices to match\n",from->tag(),from->mesh_vertices.size());
}
@@ -649,6 +635,9 @@ static void copy_vertices (GFace *geom, GFace *mesh, std::map<MVertex*,MVertex*>
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);
+ //SPoint2 param;
+ //reparamMeshVertexOnFace(v_to, to, param, true);
+ //printf("PARAMS : %g %g\n", param.x(), param.y());
_mesh_to_geom[v_from] = v_to;
}
// printf("Ending Face %d %d vertices to match\n",geom->tag(),geom->mesh_vertices.size());
@@ -665,6 +654,9 @@ static void copy_elements (std::vector<ELEMENT*> &to,
std::vector<MVertex*> nodes;
for(int j=0;j<e->getNumVertices();j++) {
nodes.push_back(_mesh_to_geom[e->getVertex(j)]);
+ if (_mesh_to_geom[e->getVertex(j)] == 0) {
+ printf("Error vertex %i\n", e->getVertex(j)->getNum());
+ }
}
to.push_back( (ELEMENT*)(toto.create(e->getTypeForMSH(), nodes) ));
}
@@ -683,20 +675,21 @@ void copy_vertices (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_me
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);
+ 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,
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 (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);
@@ -706,48 +699,33 @@ void copy_elements (GModel *geom, GModel *mesh, std::map<MVertex*,MVertex*> &_me
copy_elements<MHexahedron>((*rit)->hexahedra,mesh->getRegionByTag((*rit)->tag())->hexahedra,_mesh_to_geom);
copy_elements<MPrism>((*rit)->prisms,mesh->getRegionByTag((*rit)->tag())->prisms,_mesh_to_geom);
copy_elements<MPyramid>((*rit)->pyramids,mesh->getRegionByTag((*rit)->tag())->pyramids,_mesh_to_geom);
- }
+ }
}
-
int GeomMeshMatcher::match(GModel *geom, GModel *mesh)
{
mesh->createTopologyFromMesh();
GModel::setCurrent(geom);
+
bool ok = true;
+
// This will match VERTICES
std::vector<Pair<GVertex*, GVertex*> > *coresp_v = matchVertices(geom, mesh,ok);
- ok = true;
- if (ok) {
- // This will match EDGES
- std::vector<Pair<GEdge*, GEdge*> > *coresp_e = matchEdges(geom, mesh, coresp_v,ok);
- if (ok) {
- // This will match SURFACES
- std::vector<Pair<GFace*, GFace*> > *coresp_f = matchFaces(geom, mesh, coresp_e,ok);
- if (ok) {
- // This will match REGIONS
- //std::vector<Pair<GRegion*, GRegion*> >* coresp_r =
- matchRegions(geom, mesh, coresp_f,ok);
- if (ok) {
- std::map<MVertex*,MVertex*> _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 !");
- return 0;
- }
- } else {
- Msg::Error("Could not match every surface !");
- return 0;
- }
- } else {
- Msg::Error("Could not match every edge !");
- return 0;
- }
- } else {
- Msg::Error("Could not match every vertex !");
- return 0;
- }
+
+ // This will match EDGES
+ std::vector<Pair<GEdge*, GEdge*> > *coresp_e = matchEdges(geom, mesh, coresp_v,ok);
+
+ // This will match SURFACES
+ std::vector<Pair<GFace*, GFace*> > *coresp_f = matchFaces(geom, mesh, coresp_e,ok);
+ matchRegions(geom, mesh, coresp_f,ok);
+
+ std::map<MVertex*,MVertex*> _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);
+
+ geom->writeMSH("testing.msh", 2.2, false, true);
+
+ return 1;
+
return 0;
}
--
GitLab