diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 8248973bf92f86eeb01382f97670380f7a4e6a64..bd70e94b536f2041442a242e440fa9aed08f2148 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -271,7 +271,7 @@ double GEdgeCompound::curvatures(const double par, SVector3 *dirMax, SVector3 *d
double cMax[2];
SVector3 dMin[2];
SVector3 dMax[2];
- curvature.edgeNodalValuesAndDirections(mline, dMax, dMin, cMax, cMin, 0);
+ curvature.edgeNodalValuesAndDirections(mline, dMax, dMin, cMax, cMin, 1);
*dirMax = (1-tLocMLine)*dMax[0] + tLocMLine*dMax[1];
*dirMin = (1-tLocMLine)*dMin[0] + tLocMLine*dMin[1];
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index a7a0f011e8b3632c17d93a61342237823cd8ff02..d8a56c26025dd0c9da630ea07fb3305bf62e8cba 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -151,23 +151,23 @@ static double max_surf_curvature(const GEdge *ge, double u)
return val;
}
-/*
-static double max_surf_curvature(const GVertex *gv)
-{
- double val = 0;
- std::list<GEdge*> l_edges = gv->edges();
- for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
- ite != l_edges.end(); ++ite){
- GEdge *_myGEdge = *ite;
- Range<double> bounds = _myGEdge->parBounds(0);
- if (gv == _myGEdge->getBeginVertex())
- val = std::max(val, max_surf_curvature(_myGEdge, bounds.low()));
- else
- val = std::max(val, max_surf_curvature(_myGEdge, bounds.high()));
- }
- return val;
-}
-*/
+
+// static double max_surf_curvature_vertex(const GVertex *gv)
+// {
+// double val = 0;
+// std::list<GEdge*> l_edges = gv->edges();
+// for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
+// ite != l_edges.end(); ++ite){
+// GEdge *_myGEdge = *ite;
+// Range<double> bounds = _myGEdge->parBounds(0);
+// if (gv == _myGEdge->getBeginVertex())
+// val = std::max(val, max_surf_curvature(_myGEdge, bounds.low()));
+// else
+// val = std::max(val, max_surf_curvature(_myGEdge, bounds.high()));
+// }
+// return val;
+// }
+
SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
bool surface_isotropic,
@@ -180,10 +180,10 @@ SMetric3 metric_based_on_surface_curvature(const GFace *gf, double u, double v,
cmax = gf->curvatures(SPoint2(u, v),&dirMax, &dirMin, &cmax,&cmin);
if (cmin == 0)cmin =1.e-12;
if (cmax == 0)cmax =1.e-12;
- double lambda2 = ((2 * M_PI) /( fabs(cmax) * CTX::instance()->mesh.minCircPoints ) );
double lambda1 = ((2 * M_PI) /( fabs(cmin) * CTX::instance()->mesh.minCircPoints ) );
+ double lambda2 = ((2 * M_PI) /( fabs(cmax) * CTX::instance()->mesh.minCircPoints ) );
SVector3 Z = crossprod(dirMax,dirMin);
- if (surface_isotropic) lambda2 = lambda1 = std::min(lambda2,lambda1);
+ if (surface_isotropic) lambda2 = lambda1 = std::min(lambda2,lambda1);
dirMin.normalize();
dirMax.normalize();
Z.normalize();
@@ -209,8 +209,8 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u, boo
double cmax, cmin;
SVector3 dirMax,dirMin;
cmax = ptrCompoundEdge->curvatures(u,&dirMax, &dirMin, &cmax,&cmin);
- if (cmin == 0)cmin =1.e-5;
- if (cmax == 0)cmax =1.e-5;
+ if (cmin == 0)cmin =1.e-12;
+ if (cmax == 0)cmax =1.e-12;
double lambda2 = ((2 * M_PI) /( fabs(cmax) * CTX::instance()->mesh.minCircPoints ) );
double lambda1 = ((2 * M_PI) /( fabs(cmin) * CTX::instance()->mesh.minCircPoints ) );
SVector3 Z = crossprod(dirMax,dirMin);
@@ -221,7 +221,7 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u, boo
lambda2 = std::min(lambda2, CTX::instance()->mesh.lcMax);
SMetric3 curvMetric (1. / (lambda1 * lambda1), 1. / (lambda2 * lambda2),
- 1.e-5, dirMin, dirMax, Z);
+ 1.e-12, dirMin, dirMax, Z);
return curvMetric;
}
else{
@@ -240,6 +240,7 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u, boo
}
++it;
}
+
return curvMetric;
}
}
@@ -283,7 +284,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
switch(ge->dim()){
case 0:
Crv = max_edge_curvature((const GVertex *)ge);
- // Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
+ //Crv = std::max(max_surf_curvature_vertex((const GVertex *)ge), Crv);
// Crv = max_surf_curvature((const GVertex *)ge);
break;
case 1:
@@ -448,7 +449,8 @@ SMetric3 BGM_MeshMetric(GEntity *ge,
}
SMetric3 LC(1./(lc*lc));
- SMetric3 m = intersection_conserveM1(intersection_conserveM1 (l4, LC),l3);
+ //SMetric3 m = intersection_conserveM1(intersection_conserveM1 (l4, LC),l3);
+ SMetric3 m = intersection(intersection (l4, LC),l3);
//printf("%g %g %g %g %g %g\n",m(0,0),m(1,1),m(2,2),m(0,1),m(0,2),m(1,2));
{
fullMatrix<double> V(3,3);
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 88dbe9e8c45fab2176f7b751ae51727d17ba32c1..490e1dc09f7ec8b9ee3f628050c451d263e6145c 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -18,6 +18,7 @@
#include <fstream>
#include <sstream>
#include <iostream>
+#include <math.h>
#include "OS.h"
#include "GModel.h"
#include "MElement.h"
@@ -1078,7 +1079,7 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge)
void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge)
{
- if (update_needed){
+ if (update_needed){
std::ifstream input;
input.open(fileName.c_str());
if(StatFile(fileName))
@@ -1088,10 +1089,10 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
update_needed = false;
}
- double xyz[3] = {x,y,z};
- //take xyz = closest point on boundary in case we are on the planar IN/OUT
- //FACES or in VOLUME
+ //take xyz = closest point on boundary in case we are on
+ //the planar IN/OUT FACES or in VOLUME
+ double xyz[3] = {x,y,z};
bool onTubularSurface = true;
double ds = 0.0;
bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ?
@@ -1122,6 +1123,13 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
delete[]index2;
delete[]dist2;
+ //dir_a = direction along the centerline
+ //dir_n = normal direction of the disk
+ //dir_t = tangential direction of the disk
+ SVector3 dir_a = p1-p0; dir_a.normalize();
+ SVector3 dir_n(xyz[0]-p0.x(), xyz[1]-p0.y(), xyz[2]-p0.z()); dir_n.normalize();
+ SVector3 dir_cross = crossprod(dir_a,dir_n); dir_cross.normalize();
+
//find discrete curvature at closest vertex
Curvature& curvature = Curvature::getInstance();
if( !Curvature::valueAlreadyComputed() ) {
@@ -1131,60 +1139,75 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
}
double curv, cMin, cMax;
SVector3 dMin, dMax;
- curvature.vertexNodalValuesAndDirections(vertices[index[0]],&dMax, &dMin, &cMax, &cMin, 1);
- curvature.vertexNodalValues(vertices[index[0]], curv, 0);
-
- double sign = (curv > 0.0) ? -1.0: 1.0;
- double beta = CTX::instance()->mesh.smoothRatio; //beta = 1.25 better !
- double ratio = 1.1;
-
+ int isAbs = 1.0;
+ curvature.vertexNodalValuesAndDirections(vertices[index[0]],&dMax, &dMin, &cMax, &cMin, isAbs);
+ curvature.vertexNodalValues(vertices[index[0]], curv, 1);
+ if (cMin == 0) cMin =1.e-12;
+ if (cMax == 0) cMax =1.e-12;
+ double rhoMin = 1./cMin;
+ double rhoMax = 1./cMax;
+ double signMin = (rhoMin > 0.0) ? -1.0: 1.0;
+ double signMax = (rhoMax > 0.0) ? -1.0: 1.0;
+
+ //user-defined parameters
+ double beta = 1.2; //CTX::instance()->mesh.smoothRatio;
double thickness = radMax/3.;
- double rho = radMax;
- double hwall_n = thickness/20.;
- double hwall_t = 2*M_PI*rho/nbPoints;
- double hfar = radMax/4.;
- double lc_a = 3.*hwall_t;
-
- //dir_a = direction along the centerline
- //dir_n = normal direction of the disk
- //dir_t = tangential direction of the disk
- SVector3 dir_a = p1-p0; dir_a.normalize();
- SVector3 dir_n(xyz[0]-p0.x(), xyz[1]-p0.y(), xyz[2]-p0.z()); dir_n.normalize();
- SVector3 dir_cross = crossprod(dir_a,dir_n); dir_cross.normalize();
- SVector3 dir_t1, dir_t2, dir_a1, dir_a2;
- buildOrthoBasis2(dir_n, dir_t1, dir_t2);
- buildOrthoBasis2(dir_a, dir_a1, dir_a2);
-
- double ll1 = ds*(ratio-1) + hwall_n;
- double lc_n = std::min(ll1,hfar);
- double ll2 = hwall_t *(rho+sign*ds)/rho ; //sign * ds*(ratio-1) + hwall_t;
- if (ds > thickness) ll2 = hwall_t *(rho+sign*thickness)/rho ; //lc_a; //hfar
- double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, ll2); //std::min(ll2,hfar));
-
- lc_n = std::max(lc_n, CTX::instance()->mesh.lcMin);
- lc_n = std::min(lc_n, CTX::instance()->mesh.lcMax);
- lc_t = std::max(lc_t, CTX::instance()->mesh.lcMin);
- lc_t = std::min(lc_t, CTX::instance()->mesh.lcMax);
+ double hwall_t = 2*M_PI*radMax/nbPoints;
+ double h_far = radMax/2.;
+ double h_a = radMax; //3.*hwall_t;
+
+ //define h_n, h_t1, and h_t2
+ double h_n_0 = thickness/20.;
+ double h_n = std::min( (h_n_0+ds*log(beta)), h_far);
+ //double h_n = ds*(beta-1) + hwall_n;
+ //if (ds > thickness) h_n = hwall_t *(rho+sign*thickness)/rho ;
+
+ double oneOverD2_min = 1./(2.*rhoMin*rhoMin*(beta*beta-1)) *
+ (sqrt(1+ (4.*rhoMin*rhoMin*(beta*beta-1))/(h_n*h_n))-1.);
+ double oneOverD2_max = 1./(2.*rhoMax*rhoMax*(beta*beta-1)) *
+ (sqrt(1+ (4.*rhoMax*rhoMax*(beta*beta-1))/(h_n*h_n))-1.);
+ double h_t1_0 = sqrt(1./oneOverD2_min);
+ double h_t2_0 = sqrt(1./oneOverD2_max);
+ double h_t1 = (h_t1_0+signMin*ds*log(beta));
+ double h_t2 = (h_t2_0+signMax*ds*log(beta));
+ //double ll2 = hwall_t *(rho+sign*ds)/rho ; //sign * ds*(ratio-1) + hwall_t;
+ //if (ds > thickness) ll2 = hwall_t *(rho+sign*thickness)/rho ;
+
+ //printf("************* rhoMin =%g rhoMax=%g \n", rhoMin, rhoMax);
+ //printf("h_n_0 =%g h_n =%g\n", h_n_0 , h_n);
+ //printf("h_t1_0 =%g h_t2_0=%g ds=%g\n", h_t1_0, h_t2_0, ds);
+ //printf("h_t1 =%g h_t2=%g radMax=%g h_a=%g\n", h_t1, h_t2, radMax, h_a);
+
+ //length between min and max
+ double lcMin = ((2 * M_PI *radMax) /( 20*nbPoints )); //CTX::instance()->mesh.lcMin;
+ double lcMax = lcMin*2000.; //CTX::instance()->mesh.lcMax;
+ h_n = std::max(h_n, lcMin);
+ h_n = std::min(h_n, lcMax);
+ h_t1 = std::max(h_t1, lcMin);
+ h_t1 = std::min(h_t1, lcMax);
+ h_t2 = std::max(h_t2, lcMin);
+ h_t2 = std::min(h_t2, lcMax);
//curvature metric
SMetric3 curvMetric, curvMetric1, curvMetric2;
- if (ds <= thickness){
- metr = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, radMax,
- beta, lc_n, lc_t, hwall_t);
+ if (onInOutlets){
+ metr = buildMetricTangentToCurve(dir_n,h_n,h_t2);
}
- else if (ds > thickness && onInOutlets){
- metr = buildMetricTangentToCurve(dir_n,lc_n,lc_t);
- }
- else if (ds > thickness && !onInOutlets){
- //metr = buildMetricTangentToCurve(dir_n,lc_n,lc_t);
- //curvMetric = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, radMax, beta, lc_n, lc_t, hwall_t);
- curvMetric1 = buildMetricTangentToCurve(dir_n,lc_n,lc_a); //lc_t
- curvMetric2 = buildMetricTangentToCurve(dir_cross,lc_t,lc_a); //lc_t
- curvMetric = intersection(curvMetric1,curvMetric2);
- //metr = SMetric3(1./(lc_a*lc_a), 1./(hfar*hfar),1./(hfar*hfar), dir_a, dir_a1, dir_a2);
- metr = SMetric3(1./(lc_a*lc_a), 1./(lc_n*lc_n), 1./(lc_t*lc_t), dir_a, dir_n, dir_cross);
- metr = intersection_conserveM1(metr,curvMetric);
- //metr = intersection(metr,curvMetric);
+ else {
+ if (ds <= thickness ){
+ metr = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, h_n, h_t1, h_t2);
+ }
+ else {
+ //printf("in volume \n");
+ curvMetric = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, h_n, h_t1, h_t2);
+ //curvMetric1 = buildMetricTangentToCurve(dir_n,lc_n,lc_a); //lc_t
+ //curvMetric2 = buildMetricTangentToCurve(dir_cross,lc_t,lc_a); //lc_t
+ //curvMetric = intersection(curvMetric1,curvMetric2);
+ metr = SMetric3(1./(h_a*h_a), 1./(h_far*h_far), 1./(h_far*h_far), dir_a, dir_n, dir_cross);
+ //metr = intersection_conserveM1(metr,curvMetric);
+ metr = intersection_conserve_mostaniso(metr, curvMetric);
+ //metr = intersection(metr,curvMetric);
+ }
}
return;
@@ -1192,39 +1215,13 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
}
SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dirMax,
- double cmin, double cmax, double radMax,
- double beta, double lc_n, double lc_t,
- double hwall_t)
+ double cmin, double cmax,
+ double h_n, double h_t1, double h_t2)
{
- double lcMin = ((2 * M_PI *radMax) /( 20*nbPoints )); //CTX::instance()->mesh.lcMin;
- double lcMax = ((2 * M_PI *radMax) /( 0.05*nbPoints)); //CTX::instance()->mesh.lcMax;
- if (cmin == 0) cmin =1.e-12;
- if (cmax == 0) cmax =1.e-12;
- double lambda1 = ((2 * M_PI) /( fabs(cmin) * nbPoints ) );
- double lambda2 = ((2 * M_PI) /( fabs(cmax) * nbPoints ) );
-
- double betaM = 1*beta;
- double rhoMin = 1./cmin;
- double rhoMax = 1./cmax;
- double oneOverD2_min = 1./(2.*rhoMin*rhoMin*(betaM*betaM-1)) *
- (sqrt(1+ (4.*rhoMin*rhoMin*(betaM*betaM-1))/(lc_n*lc_n))-1.);
- double oneOverD2_max = 1./(2.*rhoMax*rhoMax*(beta*beta-1)) *
- (sqrt(1+ (4.*rhoMax*rhoMax*(beta*beta-1))/(lc_n*lc_n))-1.);
-
- lambda1 = sqrt(1./oneOverD2_min);
- lambda2 = sqrt(1./oneOverD2_max);
SVector3 dirNorm = crossprod(dirMax,dirMin);
- dirMin.normalize();
- dirMax.normalize();
- dirNorm.normalize();
- lambda1 = std::max(lambda1, lcMin);
- lambda2 = std::max(lambda2, lcMin);
- lambda1 = std::min(lambda1, lcMax);
- lambda2 = std::min(lambda2, lcMax);
-
- SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),
- 1./(lc_n*lc_n), dirMin, dirMax, dirNorm); // 1.e-6
+ SMetric3 curvMetric (1./(h_t1*h_t1),1./(h_t2*h_t2),
+ 1./(h_n*h_n), dirMin, dirMax, dirNorm);
return curvMetric;
}
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index ae0429b364f8db67ff0b727b8c10fd6570dbe2ec..cd1b6834c1cb52b77db996a0265011e86f851433 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -159,8 +159,8 @@ class Centerline : public Field{
//Print for debugging
void printSplit() const;
- SMetric3 metricBasedOnSurfaceCurvature(SVector3 dMin, SVector3 dMax, double cMin, double cMax, double radMax,
- double beta, double lc_n, double lc_t, double hwall_t);
+ SMetric3 metricBasedOnSurfaceCurvature(SVector3 dMin, SVector3 dMax, double cMin, double cMax,
+ double lc_n, double lc_t1, double lc_t2);
};
#else
diff --git a/benchmarks/centerlines/aorta_centerlines.geo b/benchmarks/centerlines/aorta_centerlines.geo
index 645bc2cf90969926931764f52ca9fcb76874a7de..99d8274a67ddfce7e6b466e9bc46ac8f67783559 100644
--- a/benchmarks/centerlines/aorta_centerlines.geo
+++ b/benchmarks/centerlines/aorta_centerlines.geo
@@ -1,16 +1,16 @@
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 7; //(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-3;
-Mesh.RecombineAll = 1;
-Mesh.Bunin = 60;
+//Mesh.RecombineAll = 1;
+//Mesh.Bunin = 60;
Merge "aorta2.stl";
Field[1] = Centerline;
Field[1].FileName = "centerlinesAORTA.vtk";
-Field[1].nbPoints = 25;
+Field[1].nbPoints = 20;
Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2;//percent of vessel radius
diff --git a/benchmarks/centerlines/bypass_centerlines.geo b/benchmarks/centerlines/bypass_centerlines.geo
index 916b8a5aaa5ecd0779f1e5e66309b0fa3ae287f8..c49694ad77b619c3888d67f3d4009f042edbd0d7 100644
--- a/benchmarks/centerlines/bypass_centerlines.geo
+++ b/benchmarks/centerlines/bypass_centerlines.geo
@@ -10,7 +10,7 @@ Merge "bypass.stl";
Field[1] = Centerline;
Field[1].FileName = "centerlinesBYPASS.msh";
-Field[1].nbPoints = 15;
+Field[1].nbPoints = 25;
Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2;//percent of vessel radius
diff --git a/benchmarks/centerlines/cyl_centerlines.geo b/benchmarks/centerlines/cyl_centerlines.geo
index 8e88efc01c467e184f870e54b97f2ecbd8dac5be..09309bff8c20d3dd7e6d8c70984327040b9aeda8 100644
--- a/benchmarks/centerlines/cyl_centerlines.geo
+++ b/benchmarks/centerlines/cyl_centerlines.geo
@@ -1,11 +1,11 @@
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 7; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 7=mmg3D
-Mesh.LcIntegrationPrecision = 1.e-2;
+Mesh.LcIntegrationPrecision = 1.e-3;
Merge "cylemi.stl";
-Mesh.RecombineAll = 1;
+//Mesh.RecombineAll = 1;
//Mesh.Bunin = 100;
Field[1] = Centerline;