Skip to content
Snippets Groups Projects
Select Git revision
  • d477254996c764a77837db3555f79298fbec36c2
  • master default protected
  • hierarchical-basis
  • alphashapes
  • bl
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • fixedMeshIF
  • save_edges
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

Levelset.cpp

Blame
  • GeomMeshMatcher.cpp 31.37 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);
    }
    
    // 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*> >;
      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;
        GEntity* choice_ge = 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;
            choice_ge = (*vit2);
            best_score = score;
          }
        }
    
        if (best_score != DBL_MAX) {
          Msg::Debug("Model Vertex %i (geom) and %i (mesh) match",
                     v1->tag(),choice_ge->tag());
          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);
            }
          }
        }
        Msg::Debug("Edges %i (geom) and %i (mesh) match.",
                   e1->tag(),
                   choice->tag());
        coresp_e->push_back(Pair<GEdge*,GEdge*>(e1, choice));
    
        // --- copy topological information
        //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));
    
          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(unsigned int i = 0; i < entities.size(); i++){
        for(unsigned int 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){
            printf("vertex %d matches GVertex %d\n",v->getNum(),gv->tag());
            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(),-1.,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 (unsigned int 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 printf("argh !\n");
          }
          else{
            if (!v1)printf("Vertex %d has not been found\n", (*it)->lines[i]->getVertex(0)->getNum());
            if (!v2)printf("Vertex %d has not been found\n", (*it)->lines[i]->getVertex(1)->getNum());
          }
        }
      }
      for (GModel::fiter it = mesh->firstFace(); it != mesh->lastFace(); ++it){
        for (unsigned int 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 (unsigned int 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::map<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 (unsigned int 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 (unsigned int 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 (unsigned int 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 (unsigned int 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);
        //SPoint2 param;
        //reparamMeshVertexOnFace(v_to, to, param, true);
        //printf("PARAMS : %g %g\n", param.x(), param.y());
        _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 (unsigned int 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 (unsigned int i=0;i<coresp_v->size();++i)
        copy_vertices((*coresp_v)[i].first(),(*coresp_v)[i].second(),_mesh_to_geom);
      for (unsigned int i=0;i<coresp_e->size();++i)
        copy_vertices((*coresp_e)[i].first(),(*coresp_e)[i].second(),_mesh_to_geom);
      for (unsigned int i=0;i<coresp_f->size();++i)
        copy_vertices((*coresp_f)[i].first(),(*coresp_f)[i].second(),_mesh_to_geom);
      for (unsigned int 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 (unsigned int 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 (unsigned int 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 (unsigned int 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 (unsigned int 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;
    }