From a38005448fdcf47639358cd126b35a9f5a310dc0 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 10 Mar 2010 14:20:28 +0000
Subject: [PATCH]

---
 Common/LuaBindings.h         |   2 +-
 Geo/GFaceCompound.cpp        |  31 +++++------
 Geo/GModel.cpp               | 102 +++++++++++++++++------------------
 Geo/discreteEdge.cpp         |  29 +++++++---
 Geo/discreteEdge.h           |   2 +-
 Mesh/multiscalePartition.cpp |   1 +
 6 files changed, 92 insertions(+), 75 deletions(-)

diff --git a/Common/LuaBindings.h b/Common/LuaBindings.h
index 78f8913036..36cbe04951 100644
--- a/Common/LuaBindings.h
+++ b/Common/LuaBindings.h
@@ -499,7 +499,7 @@ static int luaCall(lua_State *L, void (*_f)(t0)) {
 };
 
 template <>
-#if (__GNUC__< 4 || (_GNUC__ == 4 && __GNUC_MINOR__ < 3))
+#if (__GNUC__< 4 || (__GNUC__ == 4 && __GNUC_MINOR__ < 3))
 static 
 #endif
 int luaCall<void>(lua_State *L,void (*_f)()) {
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index fb88fe0c13..300bba8c1d 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -485,7 +485,7 @@ bool GFaceCompound::parametrize() const
   //-----------------
   if (_mapping == HARMONIC){
     Msg::Debug("Parametrizing surface %d with 'harmonic map'", tag()); 
-    //fillNeumannBCS();
+    fillNeumannBCS();
     parametrize(ITERU,HARMONIC); 
     parametrize(ITERV,HARMONIC);
   }
@@ -504,8 +504,8 @@ bool GFaceCompound::parametrize() const
   else if (_mapping == CONFORMAL){
     Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
     fillNeumannBCS();
-    //bool withoutFolding = parametrize_conformal_spectral() ;
-    bool withoutFolding = parametrize_conformal();
+    bool withoutFolding = parametrize_conformal_spectral() ;
+    //bool withoutFolding = parametrize_conformal();
     if ( withoutFolding == false ){
       //printStuff(); exit(1);
       Msg::Warning("$$$ Parametrization switched to harmonic map");
@@ -533,7 +533,7 @@ bool GFaceCompound::parametrize() const
 
   if (checkAspectRatio() > AR_MAX){
     Msg::Warning("Geometrical aspect ratio too high");
-    //exit(1);
+    exit(1);
     paramOK = false;
   }
 
@@ -950,7 +950,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
 	vR = ge->mesh_vertices[j];
 	vR->getParameter(0,tR);
 	if(!vR->getParameter(0,tR)) {
-	  Msg::Error("vertex vr %p not medgevertex \n", vR);
+	  Msg::Error("vertex vr %p not MedgeVertex \n", vR);
 	  Msg::Exit(1);
 	}
 	if(tLoc >= tL && tLoc <= tR){
@@ -1807,7 +1807,8 @@ bool GFaceCompound::checkTopology() const
   double D = H; 
   if (_interior_loops.size() > 0)    D =  getSizeBB(_U0); 
   int AR = (int) checkAspectRatio();
-  //int AR = (int) ceil(H/D);
+  //int AR2 = (int) ceil(H/D);
+  //int AR = std::min(AR, AR2);
 
   if (G != 0 || Nb < 1){
     correctTopo = false;
@@ -1847,8 +1848,8 @@ double GFaceCompound::checkAspectRatio() const
 
   if(allNodes.empty()) buildAllNodes();
   
-  double limit =  1.e-15;
-  double areaMin = 1.e15;
+  double limit =  1.e-17;
+  double areaMin = 1.e20;
   double area3D = 0.0;
   int nb = 0;
   std::list<GFace*>::const_iterator it = _compound.begin();
@@ -1916,34 +1917,34 @@ void GFaceCompound::coherenceNormals()
     MTriangle *t = triangles[i];
     for (int j = 0; j < 3; j++){
       MEdge me = t->getEdge(j);
-      std::map<MEdge, std::set<MTriangle*>, Less_Edge >::iterator it = edge2tris.find(me);
+      std::map<MEdge, std::set<MTriangle*, std::less<MTriangle*> >, Less_Edge >::iterator it = edge2tris.find(me);
       if (it == edge2tris.end()) {
-	std::set<MTriangle*> mySet;
+	std::set<MTriangle*, std::less<MTriangle*> > mySet;
 	mySet.insert(t);
 	edge2tris.insert(std::make_pair(me, mySet));
       }
       else{
-	std::set<MTriangle*> mySet = it->second;
+	std::set<MTriangle*, std::less<MTriangle*> > mySet = it->second;
 	mySet.insert(t);
 	it->second = mySet;
       }
     }
   }
   
-  std::set<MTriangle* > touched;
+  std::set<MTriangle* , std::less<MTriangle*> > touched;
   int iE, si, iE2, si2;
   touched.insert(triangles[0]);
   while(touched.size() != triangles.size()){
     for(unsigned int i = 0; i < triangles.size(); i++){
       MTriangle *t = triangles[i];
-      std::set<MTriangle*>::iterator it2 = touched.find(t);
+      std::set<MTriangle*, std::less<MTriangle*> >::iterator it2 = touched.find(t);
       if(it2 != touched.end()){
 	for (int j = 0; j < 3; j++){
 	  MEdge me = t->getEdge(j);
 	  t->getEdgeInfo(me, iE,si);
 	  std::map<MEdge, std::set<MTriangle*>, Less_Edge >::iterator it = edge2tris.find(me);
-	  std::set<MTriangle*> mySet = it->second;
-	  for(std::set<MTriangle*>::iterator itt = mySet.begin(); itt != mySet.end(); itt++){
+	  std::set<MTriangle*, std::less<MTriangle*> > mySet = it->second;
+	  for(std::set<MTriangle*, std::less<MTriangle*> >::iterator itt = mySet.begin(); itt != mySet.end(); itt++){
 	    if (*itt != t){
 	      (*itt)->getEdgeInfo(me,iE2,si2);
 	      if(si == si2)  (*itt)->revert();
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 66f2a059ca..b995197777 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1322,31 +1322,31 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
   }
 
   // for each discreteEdge, create Topology
-  std::map<MVertex*, MVertex*, std::less<MVertex*> > old2new;
+  std::map<GFace*, std::map<MVertex*, MVertex*, std::less<MVertex*> > > face2Vert;
+  face2Vert.clear();
   for (std::vector<discreteEdge*>::iterator it = discEdges.begin(); it != discEdges.end(); it++){
-    double t1 = Cpu();
     (*it)->createTopo();
-    double t2 = Cpu();
-    (*it)->parametrize(old2new);
-    double t3 = Cpu();
+    (*it)->parametrize(face2Vert);//region2Vert
   }
 
-  //we need to recreate lines, triangles and tets
+ //we need to recreate lines, triangles and tets
   //that contain those new MEdgeVertices
-  double t5 = Cpu();  
-  for (std::vector<discreteFace*>::iterator iFace = discFaces.begin();  iFace != discFaces.end(); iFace++){
+  //for (std::list<GFace*>::iterator iFace = lFaces.begin();  iFace != lFaces.end(); iFace++){
+  for (std::map<GFace*, std::map<MVertex*, MVertex*, std::less<MVertex*> > >::iterator iFace= face2Vert.begin(); iFace != face2Vert.end(); iFace++){
+    std::map<MVertex*, MVertex*, std::less<MVertex*> > old2new = iFace->second;
+    GFace *gf = iFace->first;
     std::vector<MTriangle*> newTriangles;
     std::vector<MQuadrangle*> newQuadrangles;
-    for (unsigned int i = 0; i < (*iFace)->getNumMeshElements(); ++i){
-      MElement *e = (*iFace)->getMeshElement(i);
+    for (unsigned int i = 0; i < gf->getNumMeshElements(); ++i){
+      MElement *e = gf->getMeshElement(i);
       int N = e->getNumVertices();
       std::vector<MVertex *> v(N);
       for(int j = 0; j < N; j++) v[j] = e->getVertex(j);
       for (int j = 0; j < N; j++){
         std::map<MVertex*, MVertex*, std::less<MVertex*> >::iterator itmap = old2new.find(v[j]);
+	MVertex *vNEW;
         if (itmap != old2new.end())  {
-	  MVertex *vNEW;
-          vNEW = itmap->second;
+	  vNEW = itmap->second;
           v[j]=vNEW;
         }
       }
@@ -1354,47 +1354,47 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
       else if ( N == 4)  newQuadrangles.push_back(new  MQuadrangle(v[0], v[1], v[2], v[3]));
       
     }
-    (*iFace)->deleteVertexArrays();
-    (*iFace)->triangles.clear();
-    (*iFace)->triangles = newTriangles;
-    (*iFace)->quadrangles.clear();
-    (*iFace)->quadrangles = newQuadrangles;
+    gf->deleteVertexArrays();
+    gf->triangles.clear();
+    gf->triangles = newTriangles;
+    gf->quadrangles.clear();
+    gf->quadrangles = newQuadrangles;
   }
 
-  for(GModel::riter iRegion = firstRegion();  iRegion != lastRegion(); iRegion++){
-     std::vector<MTetrahedron*> newTetrahedra;
-     std::vector<MHexahedron*> newHexahedra;
-     std::vector<MPrism*> newPrisms;
-     std::vector<MPyramid*> newPyramids;
-     for (unsigned int i = 0; i < (*iRegion)->getNumMeshElements(); ++i){
-       MElement *e = (*iRegion)->getMeshElement(i);
-       int N = e->getNumVertices();
-       std::vector<MVertex *> v(N);
-       for(int j = 0; j < N; j++) v[j] = e->getVertex(j);
-       for (int j = 0; j < N; j++){
-         std::map<MVertex*, MVertex*, std::less<MVertex*> >::iterator itmap = old2new.find(v[j]);
-         MVertex *vNEW;
-         if (itmap != old2new.end())  {
-           vNEW = itmap->second;
-           v[j]=vNEW;
-         }
-       }
-       if (N == 4) newTetrahedra.push_back(new MTetrahedron(v[0], v[1], v[2], v[3]));
-       else if (N == 5) newPyramids.push_back(new MPyramid(v[0], v[1], v[2], v[3], v[4]));
-       else if (N == 6) newPrisms.push_back(new MPrism(v[0], v[1], v[2], v[3], v[4], v[5]));
-       else if (N == 8) newHexahedra.push_back(new MHexahedron(v[0], v[1], v[2], v[3],
-                                                               v[4], v[5], v[6], v[7]));
-     }
-     (*iRegion)->deleteVertexArrays();
-     (*iRegion)->tetrahedra.clear();
-     (*iRegion)->tetrahedra = newTetrahedra;
-     (*iRegion)->pyramids.clear();
-     (*iRegion)->pyramids = newPyramids;
-     (*iRegion)->prisms.clear();
-     (*iRegion)->prisms = newPrisms;
-     (*iRegion)->hexahedra.clear();
-     (*iRegion)->hexahedra = newHexahedra;
-   }
+//   for(GModel::riter iRegion = firstRegion();  iRegion != lastRegion(); iRegion++){
+//      std::vector<MTetrahedron*> newTetrahedra;
+//      std::vector<MHexahedron*> newHexahedra;
+//      std::vector<MPrism*> newPrisms;
+//      std::vector<MPyramid*> newPyramids;
+//      for (unsigned int i = 0; i < (*iRegion)->getNumMeshElements(); ++i){
+//        MElement *e = (*iRegion)->getMeshElement(i);
+//        int N = e->getNumVertices();
+//        std::vector<MVertex *> v(N);
+//        for(int j = 0; j < N; j++) v[j] = e->getVertex(j);
+//        for (int j = 0; j < N; j++){
+//          std::map<MVertex*, MVertex*, std::less<MVertex*> >::iterator itmap = old2new.find(v[j]);
+//          MVertex *vNEW;
+//          if (itmap != old2new.end())  {
+//            vNEW = itmap->second;
+//            v[j]=vNEW;
+//          }
+//        }
+//        if (N == 4) newTetrahedra.push_back(new MTetrahedron(v[0], v[1], v[2], v[3]));
+//        else if (N == 5) newPyramids.push_back(new MPyramid(v[0], v[1], v[2], v[3], v[4]));
+//        else if (N == 6) newPrisms.push_back(new MPrism(v[0], v[1], v[2], v[3], v[4], v[5]));
+//        else if (N == 8) newHexahedra.push_back(new MHexahedron(v[0], v[1], v[2], v[3],
+//                                                                v[4], v[5], v[6], v[7]));
+//      }
+//      (*iRegion)->deleteVertexArrays();
+//      (*iRegion)->tetrahedra.clear();
+//      (*iRegion)->tetrahedra = newTetrahedra;
+//      (*iRegion)->pyramids.clear();
+//      (*iRegion)->pyramids = newPyramids;
+//      (*iRegion)->prisms.clear();
+//      (*iRegion)->prisms = newPrisms;
+//      (*iRegion)->hexahedra.clear();
+//      (*iRegion)->hexahedra = newHexahedra;
+//    }
 
 }
 
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 877cea98d9..bda5879bec 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -243,19 +243,18 @@ void discreteEdge::setBoundVertices()
     +---------------+--------------+----------- ... ----------+
     _pars[0]=0   _pars[1]=1    _pars[2]=2             _pars[N+1]=N+1
 */
-void discreteEdge::parametrize(std::map<MVertex*, MVertex*, std::less<MVertex*> > &old2new)
+void discreteEdge::parametrize( std::map<GFace*,  std::map<MVertex*, MVertex*, std::less<MVertex*> > > &face2Vert)
 { 
   for (unsigned int i = 0; i < lines.size() + 1; i++){
     _pars.push_back(i);
   }
 
   //Replace MVertex by MedgeVertex
-  
+  std::map<MVertex*, MVertex*, std::less<MVertex*> > old2new;	
+  old2new.clear();
   std::vector<MVertex* > newVertices;
   std::vector<MLine*> newLines;
   
-  double t1 = Cpu();
-
   MVertex *vL = getBeginVertex()->mesh_vertices[0];
   int i = 0;
   for(i = 0; i < (int)lines.size() - 1; i++){
@@ -274,15 +273,31 @@ void discreteEdge::parametrize(std::map<MVertex*, MVertex*, std::less<MVertex*>
   newLines.push_back(new MLine(vL, vR));
   _orientation[i] = 1;
 
-  mesh_vertices = newVertices;
-
   deleteVertexArrays();
   model()->destroyMeshCaches();
 
+  mesh_vertices = newVertices;
   lines.clear();
   lines = newLines;
+  
+  //  we need to recreate lines, triangles and tets
+  //  that contain those new MEdgeVertices
+  
+   for(std::list<GFace*>::iterator iFace = l_faces.begin(); iFace != l_faces.end(); ++iFace){
+     std::map<GFace*,  std::map<MVertex*, MVertex*, std::less<MVertex*> > >::iterator itmap = face2Vert.find(*iFace);
+     if (itmap == face2Vert.end()) {
+       face2Vert.insert(std::make_pair(*iFace, old2new));
+     }
+     else{
+       std::map<MVertex*, MVertex*, std::less<MVertex*> > mapVert = itmap->second;
+       for (std::map<MVertex*, MVertex*, std::less<MVertex*> >::iterator it = old2new.begin(); it != old2new.end(); it++)
+	 mapVert.insert(*it);
+       itmap->second = mapVert;
+     }
+   }
+
+  //computeNormals();
 
-  computeNormals();
  }
 
 void discreteEdge::computeNormals () const
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 63d8de1d9a..37e1a62e1a 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -25,7 +25,7 @@ class discreteEdge : public GEdge {
   virtual GPoint point(double p) const;
   virtual SVector3 firstDer(double par) const;
   virtual Range<double> parBounds(int) const;
-  void parametrize(std::map<MVertex*, MVertex*, std::less<MVertex*> > &old2new);
+  void parametrize( std::map<GFace*, std::map<MVertex*, MVertex*, std::less<MVertex*> > > &face2Verts);
   void orderMLines();
   void setBoundVertices();
   void createTopo();
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index dad3b0777e..4767606e9e 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -363,6 +363,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
       int nbParts = 2;
       Msg::Info("Mesh partition: level (%d-%d)  is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
  		nextLevel->recur,nextLevel->region, AR, NB, nbParts);  
+      //partition(*nextLevel, nbParts, MULTILEVEL);
       partition(*nextLevel, nbParts, LAPLACIAN);
     }
     else {
-- 
GitLab