From 9e32c794feecc83605b6c2710f70936e41675d88 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 12 Apr 2013 12:13:35 +0000
Subject: [PATCH] smooth first der in Gfacecompound
---
Geo/GFaceCompound.cpp | 157 +++++++++++-------
Geo/GFaceCompound.h | 1 +
Mesh/BackgroundMesh.cpp | 40 +----
Mesh/BackgroundMesh.h | 2 +-
Mesh/CenterlineField.cpp | 6 +-
Mesh/meshGFace.cpp | 2 +-
Mesh/meshGFaceLloyd.cpp | 2 +-
Mesh/meshMetric.cpp | 4 +
benchmarks/centerlines/aorta_centerlines.geo | 15 +-
.../centerlines/carotid_centerlines.geo | 9 +-
benchmarks/centerlines/cyl_centerlines.geo | 11 +-
11 files changed, 141 insertions(+), 108 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index d1cf5ac67e..a7b96fbf4f 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -913,10 +913,6 @@ bool GFaceCompound::parametrize() const
parametrize(ITERV,CONVEX);
if (_type==MEANPLANE){
checkOrientation(0, true);
- // _type = ALREADYFIXED;
- // parametrize(ITERU,CONVEX);
- // parametrize(ITERV,CONVEX);
- // checkOrientation(0, true);
}
}
// Laplace parametrization
@@ -1934,7 +1930,6 @@ SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p, bool onSurface) const
GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
{
- //printf("in remeshed oct for par =%g %g\n", par1,par2);
//if not meshed yet
if (meshStatistics.status != GFace::DONE || triangles.size()+quadrangles.size() == 0) {
@@ -2092,6 +2087,8 @@ GPoint GFaceCompound::point(double par1, double par2) const
GFaceCompoundTriangle *lt;
getTriangle(par1, par2, <, U,V);
if(!lt && _mapping != RBF){
+ printf("ARRG POINT NOT FOUND--> should improve octree search \n");
+ exit(1);
//printf("POINT no success %d tris %d quad \n", triangles.size(), quadrangles.size());
GPoint gp = pointInRemeshedOctree(par1,par2);
gp.setNoSuccess();
@@ -2171,43 +2168,74 @@ GPoint GFaceCompound::point(double par1, double par2) const
Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
{
- if(!oct) parametrize();
-
- if(trivial())
+
+ if(trivial()){
return (*(_compound.begin()))->firstDer(param);
+ }
- double U, V;
+ if(!oct) parametrize();
+
+ double U,V;
GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), <, U,V);
- if(!lt && _mapping != RBF)
+ if(!lt){
+ printf("ARRG FIRSTDER POINT NOT FOUND --> should improve octree search \n");
+ exit(1);
return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
- else if (!lt && _mapping == RBF){
- double x, y, z;
- SVector3 dXdu, dXdv ;
- _rbf->UVStoXYZ(param.x(), param.y(), x, y, z, dXdu, dXdv);
- return Pair<SVector3, SVector3>(dXdu, dXdv);
- }
-
- double mat[2][2] = {{lt->p2.x() - lt->p1.x(), lt->p3.x() - lt->p1.x()},
- {lt->p2.y() - lt->p1.y(), lt->p3.y() - lt->p1.y()}};
- double inv[2][2];
- double det = inv2x2(mat,inv);
- if (!det && _mapping == RBF){
- double x, y, z;
- SVector3 dXdu, dXdv ;
- _rbf->UVStoXYZ(param.x(), param.y(), x, y, z, dXdu, dXdv);
- return Pair<SVector3, SVector3>(dXdu, dXdv);
}
- SVector3 dXdxi(lt->v2 - lt->v1);
- SVector3 dXdeta(lt->v3 - lt->v1);
-
- SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]);
- SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]);
-
- return Pair<SVector3, SVector3>(dXdu, dXdv);
+ MTriangle *tri = lt->tri;
+ SVector3 dXdu1 = firstDerivatives[tri->getVertex(0)].first();
+ SVector3 dXdu2 = firstDerivatives[tri->getVertex(1)].first();
+ SVector3 dXdu3 = firstDerivatives[tri->getVertex(2)].first();
+ SVector3 dXdv1 = firstDerivatives[tri->getVertex(0)].second();
+ SVector3 dXdv2 = firstDerivatives[tri->getVertex(1)].second();
+ SVector3 dXdv3 = firstDerivatives[tri->getVertex(2)].second();
+ SVector3 dXdu = dXdu1*(1.-U-V) + dXdu2*U + dXdu3*V;
+ SVector3 dXdv = dXdv1*(1.-U-V) + dXdv2*U + dXdv3*V;
+ return Pair<SVector3, SVector3>(dXdu,dXdv);
+
}
+// Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
+// {
+// if(!oct) parametrize();
+
+// if(trivial())
+// return (*(_compound.begin()))->firstDer(param);
+
+// double U, V;
+// GFaceCompoundTriangle *lt;
+// getTriangle(param.x(), param.y(), <, U,V);
+// if(!lt && _mapping != RBF)
+// return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
+// else if (!lt && _mapping == RBF){
+// double x, y, z;
+// SVector3 dXdu, dXdv ;
+// _rbf->UVStoXYZ(param.x(), param.y(), x, y, z, dXdu, dXdv);
+// return Pair<SVector3, SVector3>(dXdu, dXdv);
+// }
+
+// double mat[2][2] = {{lt->p2.x() - lt->p1.x(), lt->p3.x() - lt->p1.x()},
+// {lt->p2.y() - lt->p1.y(), lt->p3.y() - lt->p1.y()}};
+// double inv[2][2];
+// double det = inv2x2(mat,inv);
+// if (!det && _mapping == RBF){
+// double x, y, z;
+// SVector3 dXdu, dXdv ;
+// _rbf->UVStoXYZ(param.x(), param.y(), x, y, z, dXdu, dXdv);
+// return Pair<SVector3, SVector3>(dXdu, dXdv);
+// }
+
+// SVector3 dXdxi(lt->v2 - lt->v1);
+// SVector3 dXdeta(lt->v3 - lt->v1);
+
+// SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]);
+// SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]);
+
+// return Pair<SVector3, SVector3>(dXdu, dXdv);
+// }
+
void GFaceCompound::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
@@ -2339,30 +2367,8 @@ void GFaceCompound::getTriangle(double u, double v,
double &_u, double &_v) const
{
double uv[3] = {u, v, 0};
- // if (_mapping == RBF){
- // std::list<void*> l;
- // Octree_SearchAll(uv, oct, &l);
- // if (l.size() > 1 || l.size() == 0){
- // GFaceCompoundTriangle *gfct = NULL;
- // *lt = gfct;
- // return;
- // }
- // else{
- // std::list<void*>::iterator it = l.begin();
- // *lt = (GFaceCompoundTriangle*)(*it);
- // }
- // }
-
*lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
- // if(!(*lt)) {
- // for(int i=0;i<nbT;i++){
- // if(GFaceCompoundInEle (&_gfct[i],uv)){
- // *lt = &_gfct[i];
- // break;
- // }
- // }
- // }
if(!(*lt)){
return;
}
@@ -2391,6 +2397,7 @@ void GFaceCompound::buildOct() const
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
+ //create bounding box
for(int j = 0; j < 3; j++){
std::map<MVertex*,SPoint3>::const_iterator itj = coordinates.find(t->getVertex(j));
_coordPoints.insert(std::make_pair(t->getVertex(j)->point(), itj->second));
@@ -2414,6 +2421,7 @@ void GFaceCompound::buildOct() const
const int maxElePerBucket = 15;
oct = Octree_Create(maxElePerBucket, origin, ssize, GFaceCompoundBB,
GFaceCompoundCentroid, GFaceCompoundInEle);
+ std::map<MElement*, Pair<SVector3,SVector3> > firstElemDerivatives;
it = _compound.begin();
count = 0;
@@ -2430,7 +2438,7 @@ void GFaceCompound::buildOct() const
_gfct[count].p1 = it0->second;
_gfct[count].p2 = it1->second;
_gfct[count].p3 = it2->second;
- if((*it)->geomType() != GEntity::DiscreteSurface){
+ if((*it)->geomType() != GEntity::DiscreteSurface){
// take care of the seam !!!!
if (t->getVertex(0)->onWhat()->dim() == 2){
reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(1), *it,
@@ -2464,6 +2472,20 @@ void GFaceCompound::buildOct() const
t->getVertex(2)->z());
_gfct[count].gf = *it;
_gfct[count].tri = t;
+
+ //compute first derivatives for every triangle
+ double mat[2][2] = {{_gfct[count].p2.x() - _gfct[count].p1.x(),
+ _gfct[count].p3.x() - _gfct[count].p1.x()},
+ {_gfct[count].p2.y() - _gfct[count].p1.y(),
+ _gfct[count].p3.y() - _gfct[count].p1.y()}};
+ double inv[2][2];
+ double det = inv2x2(mat,inv);
+ SVector3 dXdxi (_gfct[count].v2 - _gfct[count].v1);
+ SVector3 dXdeta(_gfct[count].v3 - _gfct[count].v1);
+ SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]);
+ SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]);
+ firstElemDerivatives[(MElement*)t] = Pair<SVector3,SVector3>(dXdu,dXdv);
+
Octree_Insert(&_gfct[count], oct);
count++;
}
@@ -2471,6 +2493,29 @@ void GFaceCompound::buildOct() const
nbT = count;
Octree_Arrange(oct);
+ //smooth first derivatives at vertices
+ if(adjv.size() == 0){
+ std::vector<MTriangle*> allTri;
+ std::list<GFace*>::const_iterator it = _compound.begin();
+ for( ; it != _compound.end(); ++it){
+ allTri.insert(allTri.end(), (*it)->triangles.begin(), (*it)->triangles.end() );
+ }
+ buildVertexToTriangle(allTri, adjv);
+ }
+ for(v2t_cont::iterator it = adjv.begin(); it!= adjv.end(); ++it){
+ MVertex *v = it->first;
+ std::vector<MElement*> vTri = it->second;
+ SVector3 dXdu(0.0), dXdv(0.0);
+ int nbTri = vTri.size();
+ for (unsigned int j = 0; j < nbTri; j++){
+ dXdu += firstElemDerivatives[vTri[j]].first();
+ dXdv += firstElemDerivatives[vTri[j]].second();
+ }
+ dXdu*= 1./nbTri;
+ dXdv*= 1./nbTri;
+ firstDerivatives[v] = Pair<SVector3, SVector3>(dXdu, dXdv);
+ }
+
printStuff();
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 77301d85fc..da913da8b4 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -86,6 +86,7 @@ class GFaceCompound : public GFace {
mutable v2t_cont adjv;
mutable bool mapv2Tri;
mutable std::map<MVertex*, SPoint3> coordinates;
+ mutable std::map<MVertex*, Pair<SVector3,SVector3> > firstDerivatives;
mutable std::map<MVertex*, SVector3> xuu;
mutable std::map<MVertex*, SVector3> xvv;
mutable std::map<MVertex*, SVector3> xuv;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 4bbd8a51f8..b1d8f0e097 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -545,7 +545,7 @@ backgroundMesh::backgroundMesh(GFace *_gf, bool cfd)
updateSizes(_gf);
// compute optimal mesh orientations
- propagatecrossField(_gf);
+ propagateCrossField(_gf);
_3Dto2D.clear();
_2Dto3D.clear();
@@ -770,7 +770,7 @@ inline double myAngle (const SVector3 &a, const SVector3 &b, const SVector3 &d){
}
-void backgroundMesh::propagatecrossField(GFace *_gf)
+void backgroundMesh::propagateCrossField(GFace *_gf)
{
std::map<MVertex*,double> _cosines4,_sines4;
@@ -790,24 +790,14 @@ void backgroundMesh::propagatecrossField(GFace *_gf)
reparamMeshEdgeOnFace(v[0],v[1],_gf,p1,p2);
Pair<SVector3, SVector3> der = _gf->firstDer((p1+p2)*.5);
SVector3 t1 = der.first();
- SVector3 s2 = der.second();
- SVector3 n = crossprod(t1,s2);
+ SVector3 t2 = der.second();
+ SVector3 n = crossprod(t1,t2);
n.normalize();
- SVector3 t2 (v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
+ SVector3 d1(v[1]->x()-v[0]->x(),v[1]->y()-v[0]->y(),v[1]->z()-v[0]->z());
t1.normalize();
- t2.normalize();
- double _angle = myAngle (t1,t2,n);
- // printf("GFACE %d %g %g %g %g\n",_gf->tag(),t1.x(),t1.y(),t1.z(),_angle*180/M_PI);
- // printf("angle = %12.5E\n",_angle);
- // angle_ = atan2 ( p1.y()-p2.y() , p1.x()-p2.x() );
- //double angle = atan2 ( v[0]->y()-v[1]->y() , v[0]->x()- v[1]->x() );
- //double angle = atan2 ( v0->y()-v1->y() , v0->x()- v1->x() );
+ d1.normalize();
+ double _angle = myAngle (t1,d1,n);
crossField2d::normalizeAngle (_angle);
- // SVector3 s2 = der.second();
- // s2.normalize();
- SVector3 x = t1 * cos (_angle) + crossprod(n,t1) * sin (_angle);
- // printf("angle = %g --> %g %g %g vs %g %g %g\n",_angle,x.x(),x.y(),x.z(),t2.x(),t2.y(),t2.z());
- // printf("GFACE %d GEDGE %d %g %g %g %g\n",_gf->tag(),(*it)->tag(),t1.x(),t1.y(),t1.z(),_angle*180/M_PI);
for (int i=0;i<2;i++){
std::map<MVertex*,double>::iterator itc = _cosines4.find(v[i]);
std::map<MVertex*,double>::iterator its = _sines4.find(v[i]);
@@ -1004,15 +994,6 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
v1->x(),v1->y(),v1->z(),
v2->x(),v2->y(),v2->z(),
v3->x(),v3->y(),v3->z(),itv1->second,itv2->second,itv3->second);
- /*
- fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- v1->x(),v1->y(),v1->z(),
- v2->x(),v2->y(),v2->z(),
- v3->x(),v3->y(),v3->z(),
- getAngle(v1->x(),v1->y(),v1->z()),
- getAngle(v2->x(),v2->y(),v2->z()),
- getAngle(v3->x(),v3->y(),v3->z()));
- */
}
else {
@@ -1023,13 +1004,6 @@ void backgroundMesh::print(const std::string &filename, GFace *gf,
p1.x(),p1.y(),p1.z(),
p2.x(),p2.y(),p2.z(),
p3.x(),p3.y(),p3.z(),itv1->second,itv2->second,itv3->second);
- /*
- fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
- v1->x(),v1->y(),v1->z(),
- v2->x(),v2->y(),v2->z(),
- v3->x(),v3->y(),v3->z(),
- itv1->second,itv2->second,itv3->second);
- */
}
}
fprintf(f,"};\n");
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index c8e01529a1..8ef51f6d8c 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -67,7 +67,7 @@ class backgroundMesh : public simpleFunction<double>
static void unset();
static backgroundMesh *current () { return _current; }
void propagate1dMesh(GFace *);
- void propagatecrossField(GFace *);
+ void propagateCrossField(GFace *);
void propagateCrossFieldByDistance(GFace *);
void updateSizes(GFace *);
double operator () (double u, double v, double w) const; // returns mesh size
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index aa8ec31257..445c732c30 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -294,7 +294,7 @@ static void cutTriangle(MTriangle *tri,
Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), nodesR(0)
{
- recombine = CTX::instance()->mesh.recombineAll;
+ recombine = (CTX::instance()->mesh.recombineAll) || (CTX::instance()->mesh.recombine3DAll);
index = new ANNidx[1];
dist = new ANNdist[1];
@@ -320,7 +320,7 @@ Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0)
index = new ANNidx[1];
dist = new ANNdist[1];
- recombine = CTX::instance()->mesh.recombineAll;
+ recombine = (CTX::instance()->mesh.recombineAll) || (CTX::instance()->mesh.recombine3DAll);
fileName = "centerlines.vtk";//default
nbPoints = 25;
hLayer = 0.3;
@@ -892,7 +892,7 @@ void Centerline::cutMesh()
double AR = L/D;
// printf("*** Centerline branch %d (AR=%.1f) \n", edges[i].tag, AR);
- int nbSplit = (int)floor(AR/2 + 0.9); //AR/2 + 0.9
+ int nbSplit = (int)floor(AR/2 + 1.9); //AR/2 + 0.9
if( nbSplit > 1 ){
//printf("->> cut branch in %d parts \n", nbSplit);
double li = L/nbSplit;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 2f6210170b..e2d9778a58 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1608,7 +1608,7 @@ static bool meshGeneratorElliptic(GFace *gf, bool debug = true)
center = dynamic_cast<Centerline*> (myField);
}
- bool recombine = (CTX::instance()->mesh.recombineAll);
+ bool recombine = (CTX::instance()->mesh.recombineAll) ;
int nbBoundaries = gf->edges().size();
if (center && recombine && nbBoundaries == 2) {
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index e825705e77..39d6fbbe44 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -349,7 +349,7 @@ void callback(const alglib::real_1d_array& x,double& func,alglib::real_1d_array&
}
if(start>0.0 && !error1 && !error2 && !error3){
- printf("%d %.3f\n",iteration,100.0*(start-energy)/start);
+ printf("Lloyd: %d %.3f\n",iteration,100.0*(start-energy)/start);
w->set_iteration(iteration+1);
}
else if(!error1 && !error2 && !error3){
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index d185e6ce9f..77944afc14 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -109,6 +109,7 @@ void meshMetric::updateMetrics()
_nodalMetrics[ver] = setOfMetrics[0][ver];
// _detMetric[ver] = setOfDetMetric[0][ver];
_nodalSizes[ver] = setOfSizes[0][ver];
+
if (setOfMetrics.size() > 1)
for (unsigned int i=1;i<setOfMetrics.size();i++){
_nodalMetrics[ver] = intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]);
@@ -824,6 +825,9 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
SMetric3 m4 = _nodalMetrics[e->getVertex(3)];
metr = interpolation(m1,m2,m3,m4,uvw[0],uvw[1],uvw[2]);
}
+ //if recomputeAnalyticalValues
+ //then compute exact metric
+
}
else{
Msg::Warning("point %g %g %g not found, looking for nearest node",x,y,z);
diff --git a/benchmarks/centerlines/aorta_centerlines.geo b/benchmarks/centerlines/aorta_centerlines.geo
index cd8a243f6d..5e93be28df 100644
--- a/benchmarks/centerlines/aorta_centerlines.geo
+++ b/benchmarks/centerlines/aorta_centerlines.geo
@@ -1,9 +1,12 @@
-Mesh.Algorithm = 5; //(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 = 1; //(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
+//Mesh.SmoothCrossField = 20;
+//Mesh.Smoothing=0;
+
Mesh.LcIntegrationPrecision = 1.e-3;
-//Mesh.RecombineAll = 1;
+Mesh.RecombineAll = 1;
//Mesh.Bunin = 60;
Merge "aorta.stl";
@@ -15,12 +18,12 @@ Field[1].nbPoints = 20;
Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2;//percent of vessel radius (lumen)
-Field[1].nbElemSecondLayer = 4;
-Field[1].hSecondLayer = 0.2;//percent of vessel (lumen+hlayer) radius
+//Field[1].nbElemSecondLayer = 4;
+//Field[1].hSecondLayer = 0.2;//percent of vessel (lumen+hlayer) radius
Field[1].closeVolume =1;
-Field[1].extrudeWall =1;
-//Field[1].reMesh =1;
+//Field[1].extrudeWall =1;
+Field[1].reMesh =1;
Field[1].run;
diff --git a/benchmarks/centerlines/carotid_centerlines.geo b/benchmarks/centerlines/carotid_centerlines.geo
index 51aca8bc7a..f3fa2d9754 100644
--- a/benchmarks/centerlines/carotid_centerlines.geo
+++ b/benchmarks/centerlines/carotid_centerlines.geo
@@ -1,9 +1,12 @@
-Mesh.Algorithm = 1; //(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 = 1;//(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
-Mesh.LcIntegrationPrecision = 1.e-2;
+//Mesh.SmoothCrossField = 20;
+Mesh.Smoothing=0;
-//Mesh.RecombineAll = 1;
+Mesh.LcIntegrationPrecision = 1.e-5;
+
+Mesh.RecombineAll = 1;
//Mesh.Bunin = 100;
Merge "carotid.stl";
diff --git a/benchmarks/centerlines/cyl_centerlines.geo b/benchmarks/centerlines/cyl_centerlines.geo
index b7f307b063..5bc1fc6cbf 100644
--- a/benchmarks/centerlines/cyl_centerlines.geo
+++ b/benchmarks/centerlines/cyl_centerlines.geo
@@ -1,16 +1,19 @@
-Mesh.Algorithm = 1; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
-Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 7=mmg3D
+Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad, 9=packing)
+Mesh.Algorithm3D = 9; //(1=tetgen, 4=netgen, 7=mmg3D, 9=Rtree)
+
+Mesh.RecombineAll=1;
+Mesh.Recombine3DAll=1;
+Mesh.Smoothing=0;
Mesh.LcIntegrationPrecision = 1.e-3;
Merge "cylemi.stl";
-//Mesh.RecombineAll = 1;
//Mesh.Bunin = 100;
Field[1] = Centerline;
Field[1].FileName = "centerlinesCYL.vtk";
-Field[1].nbPoints = 25;
+Field[1].nbPoints = 35;
Field[1].closeVolume =1;
Field[1].reMesh =1;
--
GitLab