From b3e15b075d8e30121e47d986e257fd84a597fe8e Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 22 Jun 2012 13:47:37 +0000
Subject: [PATCH]

---
 Geo/GFaceCompound.cpp                        |  7 +-
 Mesh/CenterlineField.cpp                     | 80 --------------------
 Mesh/CenterlineField.h                       |  4 -
 Mesh/Generator.cpp                           | 11 ---
 Mesh/meshGFace.cpp                           | 32 ++++++++
 benchmarks/centerlines/aorta_centerlines.geo |  4 +-
 6 files changed, 37 insertions(+), 101 deletions(-)

diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index ff9005b3fc..71bb65fb6f 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -902,10 +902,9 @@ bool GFaceCompound::parametrize() const
       Msg::Info("Parametrizing surface %d with 'FE conformal map'", tag());
       parametrize_conformal(0, NULL, NULL);
     }
-    printStuff(55);
-    oriented = false; //checkOrientation(0);
-    if (!oriented)  oriented = checkOrientation(0, true);
-    printStuff(77);
+    //printStuff(55);
+    oriented = checkOrientation(0, true);
+    //printStuff(77);
     if (_type==SPECTRAL &&  (!oriented  || checkOverlap(vert)) ){
       Msg::Warning("!!! parametrization switched to 'FE conformal' map");
       parametrize_conformal(0, NULL, NULL);
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 63ccf47137..b12a9968cc 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -648,82 +648,6 @@ void Centerline::createSplitCompounds()
   }
 }
 
-void Centerline::cleanMesh()
-{
-  return;  // does not work yet !
-  ////////////////////////////////
-
-  if (!is_cut || !is_closed) return;
-
-  for (int i=0; i < NF; i++){
-    std::ostringstream oss;
-    oss << "new" << NF+i+1 ;
-    std::string name = oss.str();
-    current->setPhysicalName(name, 2, i+1);
-  }
-  Msg::Info("Writing new splitted mesh mySPLITMESH.msh");
-  current->writeMSH("mySPLITMESH.msh", 2.2, false, false);
-
-  std::set<MVertex*> allNod;
-  discreteFace * mySplitMesh;
-  std::vector<std::set<MVertex*> > inOutNod;
-  std::vector<discreteFace* > inOutMesh;
-
-  mySplitMesh= new discreteFace(current, 2*NF+1);
-  mySplitMesh->addPhysicalEntity(2*NF+1);
-  current->setPhysicalName("surface mesh", 2, 2*NF+1);
-  current->add(mySplitMesh);
-  for (unsigned int i=0; i < discFaces.size(); i++){
-    GFace *gfc =  current->getFaceByTag(NF+i+1);
-    for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
-      MTriangle *t = gfc->triangles[j];
-      std::vector<MVertex *> v(3);
-      for(int k = 0; k < 3; k++){
-	v[k] = t->getVertex(k);
-	allNod.insert(v[k]);
-      }
-      mySplitMesh->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
-    }
-    for(unsigned int j = 0; j < gfc->quadrangles.size(); ++j){
-      MQuadrangle *q = gfc->quadrangles[j];
-      std::vector<MVertex *> v(4);
-      for(int k = 0; k < 4; k++){
-	v[k] = q->getVertex(k);
-	allNod.insert(v[k]);
-      }
-      mySplitMesh->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3]));
-    }
-  }
-
-  //Removing discrete Vertices - Edges - Faces
-  for (int i=0; i < NV; i++){
-    GVertex *gv = current->getVertexByTag(i+1);
-    current->remove(gv);
-  }
-  for (int i=0; i < NE; i++){
-    GEdge *ge = current->getEdgeByTag(i+1);
-    GEdge *gec = current->getEdgeByTag(NE+i+1);
-    current->remove(ge);
-    current->remove(gec);
-  }
-  for (int i=0; i < NF; i++){
-    GFace *gf  = current->getFaceByTag(i+1);
-    current->remove(gf);
-  }
-  for (unsigned int i=0; i < discFaces.size(); i++){
-    GFace *gfc = current->getFaceByTag(NF+i+1);
-    current->remove(gfc);
-  }
-
-  //Put new mesh in a new discreteFace
- for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it)
-   mySplitMesh->mesh_vertices.push_back(*it);
- mySplitMesh->meshStatistics.status = GFace::DONE;
-
- current->createTopologyFromMesh();
-
-}
-
 void Centerline::createFaces()
 {
   std::vector<std::vector<MTriangle*> > faces;
@@ -1262,10 +1186,6 @@ void  Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
 
    return;
 
-   // double lc_a_curv = sqrt(1./dot(dir_a,curvMetric,dir_a));
-   // if (ds < thickness) lc_a = std::min(lc_a, lc_a_curv);
-   // return;
-
 }
 
 SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dirMax,
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index d89c3d4f4a..ae0429b364 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -113,8 +113,6 @@ class Centerline : public Field{
 " using the following script:\n\n"
 "vmtk vmtkcenterlines -seedselector openprofiles -ifile mysurface.stl -ofile centerlines.vtp --pipe vmtksurfacewriter -ifile centerlines.vtp -ofile centerlines.vtk\n";
   }
-  
-  void cleanMesh();
 
   //isotropic operator for mesh size field function of distance to centerline
   double operator() (double x, double y, double z, GEntity *ge=0);
@@ -186,8 +184,6 @@ class Centerline : public Field{
 " using the following script:\n\n"
 "vmtk vmtkcenterlines -seedselector openprofiles -ifile mysurface.stl -ofile centerlines.vtp --pipe vmtksurfacewriter -ifile centerlines.vtp -ofile centerlines.vtk\n";
   }
-  
-  void cleanMesh();
 
   //isotropic operator for mesh size field function of distance to centerline
   double operator() (double x, double y, double z, GEntity *ge=0);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 8d2379699a..302e9b2e78 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -506,17 +506,6 @@ static void Mesh2D(GModel *m)
 
   // collapseSmallEdges(*m);
 
-#if defined(HAVE_ANN)
-  //For centerline field, clean the cut parts
-  Centerline *center = 0;
-  FieldManager *fields = GModel::current()->getFields();
-  if (fields->getBackgroundField() > 0 ){
-    Field *myField = fields->get(fields->getBackgroundField());
-    center = dynamic_cast<Centerline*> (myField);
-  }
-  if (center) center->cleanMesh();
-#endif
-
   double t2 = Cpu();
   CTX::instance()->meshTimer[1] = t2 - t1;
   Msg::StatusBar(2, true, "Done meshing 2D (%g s)", CTX::instance()->meshTimer[1]);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 4f8698d1a1..f1f96f1a89 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -24,6 +24,7 @@
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MQuadrangle.h"
+#include "CenterlineField.h"
 #include "Context.h"
 #include "GPoint.h"
 #include "GmshMessage.h"
@@ -1420,6 +1421,32 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
 
   return true;
 }
+static bool meshGeneratorElliptic(GFace *gf, bool debug = true)
+{
+
+#if defined(HAVE_ANN)
+  Centerline *center = 0;
+  FieldManager *fields = GModel::current()->getFields();
+  if (fields->getBackgroundField() > 0 ){
+    Field *myField = fields->get(fields->getBackgroundField());
+    center = dynamic_cast<Centerline*> (myField);
+  }
+  
+  bool recombine =  (CTX::instance()->mesh.recombineAll);
+  int nbBoundaries = gf->edges().size(); 
+  //printf(" nbBounds = %d  (face %d) \n", nbBoundaries, gf->tag());
+
+  if (center && recombine && nbBoundaries == 2) {
+    printf("need for elliptic grid generator \n");
+    //createRegularParamGrid();
+    //solveElliptic();
+    return true;
+  }
+  else return false;
+
+#endif
+
+}
 
 static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
 {
@@ -1857,6 +1884,11 @@ void meshGFace::operator() (GFace *gf, bool print)
   Msg::Debug("Computing edge loops");
 
   Msg::Debug("Generating the mesh");
+
+  if(meshGeneratorElliptic(gf)){
+    printf("elliptic grid generator for face %d  \n", gf->tag());
+  }
+  //else if
   if ((gf->getNativeType() != GEntity::AcisModel ||
        (!gf->periodic(0) && !gf->periodic(1))) &&
       (noSeam(gf) || gf->getNativeType() == GEntity::GmshModel ||
diff --git a/benchmarks/centerlines/aorta_centerlines.geo b/benchmarks/centerlines/aorta_centerlines.geo
index 5b7a814ab7..e95e6cb2bf 100644
--- a/benchmarks/centerlines/aorta_centerlines.geo
+++ b/benchmarks/centerlines/aorta_centerlines.geo
@@ -1,9 +1,9 @@
-Mesh.Algorithm = 7; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
 Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
 
 Mesh.LcIntegrationPrecision = 1.e-2;
 
-//Mesh.RecombineAll = 1;
+Mesh.RecombineAll = 1;
 //Mesh.Bunin = 120;
 
 Merge "aorta2.stl";
-- 
GitLab