diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index ff9005b3fcb2dca8c5ba3ac4ea62f1b354a188e2..71bb65fb6f38fb4ec35d786f786ccb75c667e5a3 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 63ccf471370d0d9c5ac9fae78e50994a3285ed7b..b12a9968cc5fd8e07855b39541ede513507dc4a6 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 d89c3d4f4a107ac5e1c613e083eabf6c870d47b0..ae0429b364f8db67ff0b727b8c10fd6570dbe2ec 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 8d2379699af0fd01516c99d2c824b3da0f366a91..302e9b2e787d69cbded6d27a33dd88b0e0fd4f54 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 4f8698d1a15b8e2372f4a4dbe75c6e9ed9389007..f1f96f1a89b587ed5e01eafd6a6b4da393ba2b75 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 5b7a814ab7e3d385bcc2a9b4a9dd45b7320d949f..e95e6cb2bf5e04ab9a2abce5f5a2a8b87ee88a8f 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";