Select Git revision
GeomMeshMatcher.cpp
GeomMeshMatcher.cpp 31.26 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
//
// Contributor(s):
// Bastien Gorissen, Koen Hillewaert
//
#include <cstdarg>
#include <algorithm>
#include <list>
#include <vector>
#include "GeomMeshMatcher.h"
#include "Pair.h"
#include "discreteVertex.h"
#include "GmshMessage.h"
#include "SOrientedBoundingBox.h"
#include "MElement.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MVertex.h"
#include "MLine.h"
#include "MPyramid.h"
#include "MTrihedron.h"
#include "MPrism.h"
#include "MPoint.h"
#include "MHexahedron.h"
#include "MTetrahedron.h"
#include "OS.h"
#include "Context.h"
#include "closestVertex.h"
GeomMeshMatcher *GeomMeshMatcher::_gmm_instance = 0;
template <class T, class container>
void getIntersection(std::vector<T> &res, std::vector<container> &lists)
{
res.clear();
container const &first_list = lists[0];
bool allsame = true;
for(typename container::const_iterator item = first_list.begin();
item != first_list.end(); item++) {
bool found = true;
for(typename std::vector<container>::iterator list_iter = lists.begin();
list_iter != lists.end(); list_iter++) {
if(*(list_iter) != first_list) {
allsame = false;
if(std::find((*list_iter).begin(), (*list_iter).end(), (*item)) ==
(*list_iter).end()) {
found = false;
}
else {
found = true;
break;
}
}
}
if(found || allsame) { res.push_back(*item); }
}
}
template <class T> T findMatching(std::vector<Pair<T, T> > &matching, T &entity)
{
for(typename std::vector<Pair<T, T> >::iterator pair = matching.begin();
pair != matching.end(); pair++) {
if((*pair).left() == entity) return ((*pair).right()); }
return (0);
}
// 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 *> >;
int num_matched_vertices = 0;
int num_total_vertices = m2->getNumVertices();
std::vector<GVertex *> vertices;
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 = CTX::instance()->geom.matchMeshTolerance;
discreteVertex *choice = 0;
double best_score = DBL_MAX;
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) {
choice = v2;
best_score = score;
}
}
if(choice && best_score != DBL_MAX) {
choice->physicals = v1->physicals;
coresp_v->push_back(Pair<GVertex *, GVertex *>(v1, choice));
num_matched_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,
std::vector<Pair<GVertex *, GVertex *> > *coresp_v,
bool &ok)
{
int num_matched_edges = 0;
int num_total_edges = m2->getNumEdges();
// Vector that will be returned.
std::vector<Pair<GEdge *, GEdge *> > *coresp_e =
new std::vector<Pair<GEdge *, GEdge *> >;
std::vector<GEdge *> closed_curves;
for(GModel::eiter eit = m1->firstEdge(); eit != m1->lastEdge(); eit++) {
GEdge *e1 = (GEdge *)*eit;
GVertex *v1 = e1->getBeginVertex();
GVertex *v2 = e1->getEndVertex();
std::vector<GEdge *> common_edges;
std::vector<std::vector<GEdge *> > 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 {
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;
if(common_edges.size() == 1) { choice = common_edges[0]; }
else {
// More than one edge between the two points ? No worries, let
// 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 = e1->getOBB();
double best_score = DBL_MAX;
// Next, let's iterate over the mesh entities.
for(std::vector<GEdge *>::iterator candidate = common_edges.begin();
candidate != common_edges.end(); candidate++) {
SOrientedBoundingBox mesh_obb = (*candidate)->getOBB();
double score = SOrientedBoundingBox::compare(geo_obb, mesh_obb);
if(score < best_score) {
best_score = score;
choice = (*candidate);
}
}
}
coresp_e->push_back(Pair<GEdge *, GEdge *>(e1, choice));
// copy topological information
if(choice){
// choice->setTag(e1->tag());
choice->physicals = e1->physicals;
}
num_matched_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,
std::vector<Pair<GEdge *, GEdge *> > *coresp_e,
bool &ok)
{
int num_matched_faces = 0;
int num_total_faces = m2->getNumFaces();
std::vector<Pair<GFace *, GFace *> > *coresp_f =
new std::vector<Pair<GFace *, GFace *> >;
for(GModel::fiter fit = m1->firstFace(); fit != m1->lastFace(); fit++) {
GFace *f1 = (GFace *)*fit;
std::vector<std::vector<GFace *> > lists;
std::vector<GEdge *> boundary_edges = f1->edges();
for(std::vector<GEdge *>::iterator boundary_edge = boundary_edges.begin();
boundary_edge != boundary_edges.end(); boundary_edge++) {
// if (boundary_edge->getBeginVertex() ==
// boundary_edge->getEndVertex() &&
if(!(*boundary_edge)->isSeam(f1)) {
GEdge *ge = findMatching<GEdge *>(*coresp_e, *boundary_edge);
if(!ge) {
Msg::Error(
"Could not find matching edge %i in face %i during matching",
(*boundary_edge)->tag(), f1->tag());
}
lists.push_back(ge->faces());
}
}
std::vector<GFace *> common_faces;
getIntersection<GFace *>(common_faces, lists);
GFace *choice = 0;
if(common_faces.size() == 0) {
Msg::Debug("Could not match face %i (geom).", f1->tag());
continue;
}
if(common_faces.size() == 1) { choice = common_faces[0]; }
else {
// Then, compute the minimal bounding box
SOrientedBoundingBox geo_obb = f1->getOBB();
double best_score = DBL_MAX;
// Next, let's iterate over the mesh entities.
for(std::vector<GFace *>::iterator candidate = common_faces.begin();
candidate != common_faces.end(); candidate++) {
SOrientedBoundingBox mesh_obb = (*candidate)->getOBB();
Msg::Info("Comparing score : %f",
SOrientedBoundingBox::compare(geo_obb, mesh_obb));
double score = SOrientedBoundingBox::compare(geo_obb, mesh_obb);
if(score < best_score) {
best_score = score;
choice = (*candidate);
}
} }
if(choice) {
Msg::Debug("Faces %i (geom) and %i (mesh) match.", f1->tag(),
choice->tag());
coresp_f->push_back(Pair<GFace *, GFace *>(f1, choice));
// copy topological information
choice->setTag(f1->tag());
f1->physicals = choice->physicals;
num_matched_faces++;
}
}
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)
{
int num_matched_regions = 0;
int num_total_regions = 0;
std::vector<Pair<GRegion *, GRegion *> > *coresp_r =
new std::vector<Pair<GRegion *, GRegion *> >;
std::vector<GEntity *> m1_entities;
m1->getEntities(m1_entities, 3);
std::vector<GEntity *> m2_entities;
m2->getEntities(m2_entities, 3);
if(m1_entities.empty() || m2_entities.empty()) {
Msg::Info(
"No regions could be matched since one of the models doesn't have any");
return coresp_r;
}
for(std::vector<GEntity *>::iterator entity1 = m1_entities.begin();
entity1 != m1_entities.end(); entity1++) {
// if ((*entity1)->dim() != 3) continue;
num_total_regions++;
// std::vector<list<GRegion*> > lists;
std::vector<GFace *> boundary_faces = ((GFace *)(*entity1))->faces();
std::vector<GFace *> coresp_bound_faces;
std::vector<GRegion *> common_regions;
for(std::vector<GFace *>::iterator boundary_face = boundary_faces.begin();
boundary_face != boundary_faces.end(); boundary_face++) {
coresp_bound_faces.push_back(
findMatching<GFace *>(*coresp_f, *boundary_face));
}
for(std::vector<GEntity *>::iterator entity2 = m2_entities.begin();
entity2 != m2_entities.end(); entity2++) {
if((*entity2)->dim() != 3) continue;
std::vector<std::vector<GFace *> > lists;
lists.push_back(coresp_bound_faces);
lists.push_back(((GRegion *)*entity2)->faces());
std::vector<GFace *> common_faces;
getIntersection<GFace *>(common_faces, lists);
if(common_faces.size() == coresp_bound_faces.size()) {
common_regions.push_back((GRegion *)*entity2);
}
}
if(common_regions.size() == 1) {
coresp_r->push_back(
Pair<GRegion *, GRegion *>((GRegion *)*entity1, common_regions[0]));
common_regions[0]->setTag(((GRegion *)*entity1)->tag());
(*entity1)->physicals = common_regions[0]->physicals;
num_matched_regions++;
}
else if(common_regions.size() > 1) {
// So, first step is to build an array of points taken on the geo entity
/*
This is made in a backward fashion compared to the other entities...
*/
std::vector<GEdge *> boundaries = ((GRegion *)*entity1)->edges();
// Then, compute the minimal bounding box
SOrientedBoundingBox geo_obb = ((GRegion *)*entity1)->getOBB();
GRegion *choice = 0;
double best_score = DBL_MAX;
// Next, let's iterate over the mesh entities.
for(std::vector<GRegion *>::iterator candidate = common_regions.begin();
candidate != common_regions.end(); candidate++) {
// Again, build an array with the vertices.
SOrientedBoundingBox mesh_obb = (*candidate)->getOBB();
Msg::Info("Comparing score : %f",
SOrientedBoundingBox::compare(geo_obb, mesh_obb));
double score = SOrientedBoundingBox::compare(geo_obb, mesh_obb);
if(score < best_score) {
best_score = score;
choice = (*candidate);
}
}
coresp_r->push_back(
Pair<GRegion *, GRegion *>((GRegion *)*entity1, choice));
if(choice){
choice->setTag(((GRegion *)*entity1)->tag());
(*entity1)->physicals = choice->physicals;
}
num_matched_regions++;
}
}
Msg::Info("Regions matched : %i / %i", num_matched_regions,
num_total_regions);
if(num_matched_regions != num_total_regions) ok = false;
return coresp_r;
}
// Public
GeomMeshMatcher *GeomMeshMatcher::instance()
{
if(!GeomMeshMatcher::_gmm_instance) {
GeomMeshMatcher::_gmm_instance = new GeomMeshMatcher();
}
return (GeomMeshMatcher::_gmm_instance);
}
void GeomMeshMatcher::destroy()
{
if(GeomMeshMatcher::_gmm_instance) delete GeomMeshMatcher::_gmm_instance;
}
static GVertex *getGVertex(MVertex *v1, GModel *gm, const double TOL)
{
GVertex *best = 0;
double bestScore = TOL;
for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) { {
GVertex *v2 = (*it)->getBeginVertex();
double score = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
(v1->y() - v2->y()) * (v1->y() - v2->y()) +
(v1->z() - v2->z()) * (v1->z() - v2->z()));
if(score < bestScore) {
bestScore = score;
best = v2;
}
}
{
GVertex *v2 = (*it)->getEndVertex();
double score = sqrt((v1->x() - v2->x()) * (v1->x() - v2->x()) +
(v1->y() - v2->y()) * (v1->y() - v2->y()) +
(v1->z() - v2->z()) * (v1->z() - v2->z()));
if(score < bestScore) {
bestScore = score;
best = v2;
}
}
}
return best;
}
static GPoint getGEdge(MVertex *v1, GModel *gm, const double TOL)
{
GPoint gpBest;
double bestScore = TOL;
for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
GEdge *e = *it;
double pp;
GPoint gp = e->closestPoint(SPoint3(v1->x(), v1->y(), v1->z()), pp);
double score = sqrt((v1->x() - gp.x()) * (v1->x() - gp.x()) +
(v1->y() - gp.y()) * (v1->y() - gp.y()) +
(v1->z() - gp.z()) * (v1->z() - gp.z()));
if(score < bestScore) {
bestScore = score;
gpBest = gp;
}
}
return gpBest;
}
static GPoint getGFace(MVertex *v1, GModel *gm, const double TOL)
{
GPoint gpBest;
double bestScore = TOL;
for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
GFace *gf = *it;
SPoint2 pp;
double guess[2] = {0, 0};
GPoint gp = gf->closestPoint(SPoint3(v1->x(), v1->y(), v1->z()), guess);
double score = sqrt((v1->x() - gp.x()) * (v1->x() - gp.x()) +
(v1->y() - gp.y()) * (v1->y() - gp.y()) +
(v1->z() - gp.z()) * (v1->z() - gp.z()));
if(score < bestScore) {
bestScore = score;
gpBest = gp;
}
}
return gpBest;
}
int GeomMeshMatcher::forceTomatch(GModel *geom, GModel *mesh, const double TOL)
{
// assume that the geometry is the right one
std::vector<GEntity *> entities; mesh->getEntities(entities);
for(std::size_t i = 0; i < entities.size(); i++) {
for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++) {
MVertex *v = entities[i]->mesh_vertices[j];
GVertex *gv = getGVertex(v, geom, TOL);
bool found = 0;
if(gv) {
found = 1;
MVertex *vvv = new MVertex(v->x(), v->y(), v->z(), gv, v->getNum());
gv->mesh_vertices.push_back(vvv);
gv->points.push_back(new MPoint(vvv, v->getNum()));
}
else if(v->onWhat()->dim() == 1) {
GPoint gp = getGEdge(v, geom, 1.e22);
if(gp.g()) {
GEntity *gg = (GEntity *)gp.g();
found = 1;
gg->mesh_vertices.push_back(new MEdgeVertex(
gp.x(), gp.y(), gp.z(), gg, gp.u(), v->getNum()));
}
}
if(!found && v->onWhat()->dim() <= 2) {
GPoint gp = getGFace(v, geom, TOL);
if(gp.g()) {
GEntity *gg = (GEntity *)gp.g();
found = 1;
gg->mesh_vertices.push_back(new MFaceVertex(
gp.x(), gp.y(), gp.z(), gg, gp.u(), gp.v(), v->getNum()));
}
}
if(!found)
Msg::Error("vertex %d classified on %d %d not matched", v->getNum(),
v->onWhat()->dim(), v->onWhat()->tag());
}
}
for(GModel::eiter it = mesh->firstEdge(); it != mesh->lastEdge(); ++it) {
for(std::size_t i = 0; i < (*it)->lines.size(); i++) {
MVertex *v1 =
geom->getMeshVertexByTag((*it)->lines[i]->getVertex(0)->getNum());
MVertex *v2 =
geom->getMeshVertexByTag((*it)->lines[i]->getVertex(1)->getNum());
if(v1 && v2) {
GEdge *ge = 0;
if(v1->onWhat()->dim() == 1) ge = (GEdge *)v1->onWhat();
if(v2->onWhat()->dim() == 1) ge = (GEdge *)v2->onWhat();
if(ge) {
double u1, u2;
reparamMeshVertexOnEdge(v1, ge, u1);
reparamMeshVertexOnEdge(v2, ge, u2);
if(u1 < u2)
ge->lines.push_back(new MLine(v1, v2));
else
ge->lines.push_back(new MLine(v2, v1));
}
else
Msg::Error("Could not find curve");
}
else {
if(!v1)
Msg::Error("Could not find node %lu",
(*it)->lines[i]->getVertex(0)->getNum());
if(!v2)
Msg::Error("Could not find node %lu",
(*it)->lines[i]->getVertex(1)->getNum());
}
}
}
for(GModel::fiter it = mesh->firstFace(); it != mesh->lastFace(); ++it) {
for(std::size_t i = 0; i < (*it)->triangles.size(); i++) {
MVertex *v1 = geom->getMeshVertexByTag((*it)->triangles[i]->getVertex(0)->getNum());
MVertex *v2 =
geom->getMeshVertexByTag((*it)->triangles[i]->getVertex(1)->getNum());
MVertex *v3 =
geom->getMeshVertexByTag((*it)->triangles[i]->getVertex(2)->getNum());
if(v1->onWhat()->dim() == 2)
((GFace *)v1->onWhat())->triangles.push_back(new MTriangle(v1, v2, v3));
else if(v2->onWhat()->dim() == 2)
((GFace *)v2->onWhat())->triangles.push_back(new MTriangle(v1, v2, v3));
else if(v3->onWhat()->dim() == 2)
((GFace *)v3->onWhat())->triangles.push_back(new MTriangle(v1, v2, v3));
}
for(std::size_t i = 0; i < (*it)->quadrangles.size(); i++) {
MVertex *v1 =
geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(0)->getNum());
MVertex *v2 =
geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(1)->getNum());
MVertex *v3 =
geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(2)->getNum());
MVertex *v4 =
geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(3)->getNum());
if(v1->onWhat()->dim() == 2)
((GFace *)v1->onWhat())
->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
else if(v2->onWhat()->dim() == 2)
((GFace *)v2->onWhat())
->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
else if(v3->onWhat()->dim() == 2)
((GFace *)v3->onWhat())
->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
else if(v4->onWhat()->dim() == 2)
((GFace *)v4->onWhat())
->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
}
}
geom->writeMSH("hopla.msh", 2.2, false, false, true);
return 0;
}
template <class GEType>
static void copy_periodicity(std::vector<Pair<GEType *, GEType *> > &eCor,
std::map<MVertex *, MVertex *> &mesh_to_geom)
{
typename std::multimap<GEType *, GEType *> eMap; // (eCor.begin(),eCor.end());
typename std::vector<Pair<GEType *, GEType *> >::iterator eIter =
eCor.begin();
for(; eIter != eCor.end(); ++eIter) {
eMap.insert(std::make_pair(eIter->second(), eIter->first()));
}
typename std::multimap<GEType *, GEType *>::iterator srcIter = eMap.begin();
for(; srcIter != eMap.end(); ++srcIter) {
GEType *oldTgt = srcIter->first;
GEType *oldSrc = dynamic_cast<GEType *>(oldTgt->getMeshMaster());
if(oldSrc != NULL && oldSrc != oldTgt) {
GEType *newTgt = srcIter->second;
typename std::multimap<GEType *, GEType *>::iterator tgtIter =
eMap.find(oldSrc);
if(tgtIter == eMap.end()) {
Msg::Error("Could not find matched entity for %d",
"which has a matched periodic counterpart %d", oldSrc->tag(),
oldTgt->tag());
}
GEType *newSrc = tgtIter->second;
newTgt->setMeshMaster(newSrc, oldTgt->affineTransform);
std::map<MVertex *, MVertex *> &oldV2v = oldTgt->correspondingVertices; std::map<MVertex *, MVertex *> &newV2v = newTgt->correspondingVertices;
std::map<MVertex *, MVertex *>::iterator vIter = oldV2v.begin();
for(; vIter != oldV2v.end(); ++vIter) {
MVertex *oldTgtV = vIter->first;
MVertex *oldSrcV = vIter->second;
std::map<MVertex *, MVertex *>::iterator newTvIter =
mesh_to_geom.find(oldTgtV);
std::map<MVertex *, MVertex *>::iterator newSvIter =
mesh_to_geom.find(oldSrcV);
if(newTvIter == mesh_to_geom.end()) {
Msg::Error(
"Could not find copy of target vertex %d in entity %d of dim",
oldTgtV->getIndex(), oldTgt->tag(), oldTgt->dim());
}
if(newSvIter == mesh_to_geom.end()) {
Msg::Error(
"Could not find copy of source vertex %d in entity %d of dim",
oldSrcV->getIndex(), oldSrc->tag(), oldSrc->dim());
}
newV2v[newTvIter->second] = newSvIter->second;
}
}
}
}
template <class GEType>
static bool apply_periodicity(std::vector<Pair<GEType *, GEType *> > &eCor)
{
typename std::multimap<GEType *, GEType *> eMap; // (eCor.begin(),eCor.end());
typename std::vector<Pair<GEType *, GEType *> >::iterator eIter =
eCor.begin();
for(; eIter != eCor.end(); ++eIter) {
eMap.insert(std::make_pair(eIter->second(), eIter->first()));
}
typename std::multimap<GEType *, GEType *>::iterator srcIter = eMap.begin();
int dim = -1;
for(; srcIter != eMap.end(); ++srcIter) {
GEType *newTgt = srcIter->second;
newTgt->updateCorrespondingVertices();
newTgt->alignElementsWithMaster();
if(dim == -1) dim = newTgt->dim();
}
if(dim < 2) { // required for multiple periodic directions
for(srcIter = eMap.begin(); srcIter != eMap.end(); ++srcIter) {
GEType *newTgt = srcIter->second;
newTgt->copyMasterCoordinates();
newTgt->alignElementsWithMaster();
}
}
return false;
}
static void copy_vertices(GVertex *to, GVertex *from,
std::map<MVertex *, MVertex *> &_mesh_to_geom)
{
to->deleteMesh();
if(from) {
for(std::size_t i = 0; i < 1; i++) {
MVertex *v_from = from->mesh_vertices[i];
MVertex *v_to = new MVertex(v_from->x(), v_from->y(), v_from->z(), to);
to->mesh_vertices.push_back(v_to); _mesh_to_geom[v_from] = v_to;
}
}
}
static void copy_vertices(GRegion *to, GRegion *from,
std::map<MVertex *, MVertex *> &_mesh_to_geom)
{
to->deleteMesh();
if(from) {
for(std::size_t i = 0; i < from->mesh_vertices.size(); i++) {
MVertex *v_from = from->mesh_vertices[i];
MVertex *v_to = new MVertex(v_from->x(), v_from->y(), v_from->z(), to);
to->mesh_vertices.push_back(v_to);
_mesh_to_geom[v_from] = v_to;
}
}
}
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;
}
for(std::size_t i = 0; i < from->mesh_vertices.size(); i++) {
MVertex *v_from = from->mesh_vertices[i];
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;
}
}
static void copy_vertices(GFace *to, GFace *from,
std::map<MVertex *, MVertex *> &_mesh_to_geom)
{
for(std::size_t i = 0; i < from->mesh_vertices.size(); i++) {
MVertex *v_from = from->mesh_vertices[i];
double uv[2];
GPoint gp =
to->closestPoint(SPoint3(v_from->x(), v_from->y(), v_from->z()), uv);
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-1)
Msg::Error("Impossible to match 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(), to->tag(), gp.x(),
gp.y(), gp.z());
MFaceVertex *v_to = new MFaceVertex(v_from->x(), v_from->y(), v_from->z(),
to, gp.u(), gp.v()); to->mesh_vertices.push_back(v_to);
_mesh_to_geom[v_from] = v_to;
}
}
template <class ELEMENT>
static void copy_elements(std::vector<ELEMENT *> &to,
std::vector<ELEMENT *> &from,
std::map<MVertex *, MVertex *> &_mesh_to_geom)
{
MElementFactory toto;
to.clear();
for(std::size_t i = 0; i < from.size(); i++) {
ELEMENT *e = from[i];
std::vector<MVertex *> nodes;
for(std::size_t j = 0; j < e->getNumVertices(); j++) {
std::map<MVertex *, MVertex *>::iterator vIter =
_mesh_to_geom.find(e->getVertex(j));
if(vIter == _mesh_to_geom.end()) {
Msg::Error("Could not find match for vertex %i during element copy "
"while matching discrete to actual CAD",
e->getVertex(j)->getNum());
}
else
nodes.push_back(vIter->second);
}
to.push_back((ELEMENT *)(toto.create(e->getTypeForMSH(), nodes)));
}
}
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,
std::vector<Pair<GRegion *, GRegion *> > *coresp_r)
{
// copy all elements
for(std::size_t i = 0; i < coresp_v->size(); ++i)
copy_vertices((*coresp_v)[i].first(), (*coresp_v)[i].second(),
_mesh_to_geom);
for(std::size_t i = 0; i < coresp_e->size(); ++i)
copy_vertices((*coresp_e)[i].first(), (*coresp_e)[i].second(),
_mesh_to_geom);
for(std::size_t i = 0; i < coresp_f->size(); ++i)
copy_vertices((*coresp_f)[i].first(), (*coresp_f)[i].second(),
_mesh_to_geom);
for(std::size_t i = 0; i < coresp_r->size(); ++i)
copy_vertices((*coresp_r)[i].first(), (*coresp_r)[i].second(),
_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,
std::vector<Pair<GRegion *, GRegion *> > *coresp_r)
{
// copy all elements
for(std::size_t i = 0; i < coresp_v->size(); ++i) {
GVertex *dest = (*coresp_v)[i].first();
GVertex *orig = (*coresp_v)[i].second();
copy_elements<MPoint>(dest->points, orig->points, _mesh_to_geom);
}
for(std::size_t i = 0; i < coresp_e->size(); ++i) {
GEdge *dest = (*coresp_e)[i].first();
GEdge *orig = (*coresp_e)[i].second(); copy_elements<MLine>(dest->lines, orig->lines, _mesh_to_geom);
}
for(std::size_t i = 0; i < coresp_f->size(); ++i) {
GFace *dest = (*coresp_f)[i].first();
GFace *orig = (*coresp_f)[i].second();
copy_elements<MTriangle>(dest->triangles, orig->triangles, _mesh_to_geom);
copy_elements<MQuadrangle>(dest->quadrangles, orig->quadrangles,
_mesh_to_geom);
}
for(std::size_t i = 0; i < coresp_r->size(); ++i) {
GRegion *dest = (*coresp_r)[i].first();
GRegion *orig = (*coresp_r)[i].second();
copy_elements<MTetrahedron>(dest->tetrahedra, orig->tetrahedra,
_mesh_to_geom);
copy_elements<MHexahedron>(dest->hexahedra, orig->hexahedra, _mesh_to_geom);
copy_elements<MPrism>(dest->prisms, orig->prisms, _mesh_to_geom);
copy_elements<MPyramid>(dest->pyramids, orig->pyramids, _mesh_to_geom);
copy_elements<MTrihedron>(dest->trihedra, orig->trihedra, _mesh_to_geom);
}
}
int GeomMeshMatcher::match(GModel *geom, GModel *mesh)
{
Msg::StatusBar(true, "Matching discrete mesh to actual CAD ...");
double t1 = Cpu();
GModel::setCurrent(geom);
geom->setPhysicalNames(mesh->getPhysicalNames());
bool ok = true;
std::vector<Pair<GVertex *, GVertex *> > *coresp_v(NULL);
std::vector<Pair<GEdge *, GEdge *> > *coresp_e(NULL);
std::vector<Pair<GFace *, GFace *> > *coresp_f(NULL);
std::vector<Pair<GRegion *, GRegion *> > *coresp_r(NULL);
coresp_v = matchVertices(geom, mesh, ok);
if(ok) {
coresp_e = matchEdges(geom, mesh, coresp_v, ok);
if(ok) {
coresp_f = matchFaces(geom, mesh, coresp_e, ok);
if(ok) { coresp_r = matchRegions(geom, mesh, coresp_f, ok); }
else
Msg::Error("Could only match %i faces out of %i ... stopping match",
coresp_f->size(), mesh->getNumFaces());
}
else
Msg::Error("Could only match %i edges out of %i ... stopping match",
coresp_e->size(), mesh->getNumEdges());
}
else
Msg::Error("Could only match %i nodes out of %i ... "
"check mesh/CAD or increase Geom.MatchMeshTolerance (now %g)",
coresp_v->size(), mesh->getNumVertices(),
CTX::instance()->geom.matchMeshTolerance);
std::map<MVertex *, MVertex *> _mesh_to_geom;
if(ok) {
copy_vertices(geom, mesh, _mesh_to_geom, coresp_v, coresp_e, coresp_f,
coresp_r);
double t00 = Cpu();
copy_elements(geom, mesh, _mesh_to_geom, coresp_v, coresp_e, coresp_f,
coresp_r);
Msg::Info("Copying mesh elements to CAD model entities took %g s",
Cpu() - t00);
t00 = Cpu();
if(!apply_periodicity(*coresp_v))
copy_periodicity(*coresp_v, _mesh_to_geom);
if(!apply_periodicity(*coresp_e))
copy_periodicity(*coresp_e, _mesh_to_geom);
if(!apply_periodicity(*coresp_f))
copy_periodicity(*coresp_f, _mesh_to_geom);
Msg::Info("Applying periodicity to CAD model entities took %g s",
Cpu() - t00);
}
if(coresp_v) delete coresp_v;
if(coresp_e) delete coresp_e;
if(coresp_f) delete coresp_f;
if(coresp_r) delete coresp_r;
if(ok)
Msg::StatusBar(true, "Matched successfully mesh to CAD in %g s",
Cpu() - t1);
else
Msg::Error("Failed to match mesh to CAD, please check");
return ok ? 1 : 0;
}