diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp index 4e7204aaba7ad4fea953b9e089ff9eea22927e48..ecfdb186914400ff970a82d48bb0d904c80a1ebc 100644 --- a/Geo/STensor3.cpp +++ b/Geo/STensor3.cpp @@ -36,6 +36,7 @@ SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2) double l1 = std::max(dot(v1,m1,v1),dot(v1,m2,v1)); double l2 = std::max(dot(v2,m1,v2),dot(v2,m2,v2)); SMetric3 iv(l0,l1,l2,v0,v1,v2); + return iv; } @@ -62,13 +63,11 @@ SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2 fullVector<double> S1(3); m1.eig(V1,S1,true); double lambda1_min = std::min(std::min(fabs(S1(0)),fabs(S1(1))),fabs(S1(2))); - //double lambda1_max = std::max(std::max(fabs(S1(0)),fabs(S1(1))),fabs(S1(2))); fullMatrix<double> V2(3,3); fullVector<double> S2(3); m2.eig(V2,S2,true); double lambda2_min = std::min(std::min(fabs(S2(0)),fabs(S2(1))),fabs(S2(2))); - //double lambda2_max = std::max(std::max(fabs(S2(0)),fabs(S2(1))),fabs(S2(2))); - + if (lambda1_min < lambda2_min) return intersection_conserveM1(m1, m2); else @@ -77,8 +76,8 @@ SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2 // (1-t) * m1 + t * m2 SMetric3 interpolation (const SMetric3 &m1, - const SMetric3 &m2, - const double t) + const SMetric3 &m2, + const double t) { SMetric3 im1 = m1.invert(); SMetric3 im2 = m2.invert(); diff --git a/Geo/STensor3.h b/Geo/STensor3.h index 92243848ed60894e46bf6f625b817b03e05403e6..be344b5b1b695fabe01e9c632554182c2fcdb25e 100644 --- a/Geo/STensor3.h +++ b/Geo/STensor3.h @@ -62,7 +62,6 @@ class SMetric3 { e(1,0) = t2(0); e(1,1) = t2(1); e(1,2) = t2(2); e(2,0) = t3(0); e(2,1) = t3(1); e(2,2) = t3(2); e.transposeInPlace(); - // e.invertInPlace(); fullMatrix<double> tmp(3,3); tmp(0,0) = l1 * e(0,0); @@ -83,6 +82,7 @@ class SMetric3 { _val[3] = tmp(2,0) * e(0,0) + tmp(2,1) * e(1,0) + tmp(2,2) * e(2,0); _val[4] = tmp(2,0) * e(0,1) + tmp(2,1) * e(1,1) + tmp(2,2) * e(2,1); _val[5] = tmp(2,0) * e(0,2) + tmp(2,1) * e(1,2) + tmp(2,2) * e(2,2); + } inline double &operator()(int i, int j) { @@ -166,6 +166,7 @@ inline double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b) // preserve orientation of m1 SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2); + // preserve orientation of the most anisotropic metric SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2); // compute the largest inscribed ellipsoid... diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp index 384d9b33039d124bac615a77cde1dff718adfe41..d11e585726fa150b67b1cb42a25b4b044d5b9782 100644 --- a/Mesh/CenterlineField.cpp +++ b/Mesh/CenterlineField.cpp @@ -1150,16 +1150,14 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt double signMax = (rhoMax > 0.0) ? -1.0: 1.0; //user-defined parameters - double thickness = radMax/3.; - double hwall_t = 2*M_PI*radMax/nbPoints; - double h_far = radMax/4.; - //define h_n, h_t1, and h_t2 + double thickness = radMax/2.; + double h_far = radMax/5.; double beta = (ds <= thickness) ? 1.2 : 2.8; //CTX::instance()->mesh.smoothRatio; double dist = (ds <= thickness) ? ds: thickness; + double h_n_0 = thickness/20.; double h_n = std::min( (h_n_0+ds*log(beta)), h_far); - double h_a_0 = radMax; double betaMin = 10.; double betaMax = 3.2; @@ -1167,23 +1165,21 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt (sqrt(1+ (4.*rhoMin*rhoMin*(betaMin*betaMin-1))/(h_n*h_n))-1.); double oneOverD2_max = 1./(2.*rhoMax*rhoMax*(betaMax*betaMax-1)) * (sqrt(1+ (4.*rhoMax*rhoMax*(betaMax*betaMax-1))/(h_n*h_n))-1.); + double l_t1 = ((2 * M_PI) /(cMin*nbPoints)); + double l_t2 = ((2 * M_PI) /(cMax*nbPoints)); double h_t1_0 = sqrt(1./oneOverD2_min); double h_t2_0 = sqrt(1./oneOverD2_max); - double h_t1 = h_t1_0*(rhoMin+signMin*dist)/rhoMin ; - double h_t2 = h_t2_0*(rhoMax+signMax*dist)/rhoMax ; - //double h_t1 = std::min((h_t1_0+dist*log(beta)), radMax); - //double h_t2 = std::min((h_t2_0+dist*log(beta)), radMax); - - double dCenter = radMax -ds; - double h_a = h_a_0 - (h_a_0-h_t1)/(radMax-thickness)*dCenter; + //double h_t1 = h_t1_0*(rhoMin+signMin*dist)/rhoMin ; + //double h_t2 = h_t2_0*(rhoMax+signMax*dist)/rhoMax ; + double h_t1 = std::min( (h_t1_0+(dist*log(beta))), radMax); + double h_t2 = std::min( (h_t2_0+(dist*log(beta))), h_far); + + double dCenter = radMax-ds; + double h_a_0 = 0.5*radMax; + double h_a = h_a_0 - (h_a_0-h_t1_0)/(radMax)*dCenter; - //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 lcMin = ((2 * M_PI *radMax) /( 50*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); @@ -1191,19 +1187,24 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt //curvature metric SMetric3 curvMetric, curvMetric1, curvMetric2; + SMetric3 centMetric1, centMetric2, centMetric; if (onInOutlets){ metr = buildMetricTangentToCurve(dir_n,h_n,h_t2); } else { - if (ds < thickness ){ + //on surface and in volume boundary layer + if ( ds < thickness ){ metr = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, h_n, h_t1, h_t2); } - else { - curvMetric = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, h_n, h_t1, h_t2); - metr = SMetric3( 1./(h_a_0*h_a_0), 1./(h_n*h_n), 1./(h_n*h_n), dir_a, dir_n, dir_cross); + //in volume + else { + //curvMetric = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, h_n, h_t1, h_t2); + metr = SMetric3( 1./(h_a*h_a), 1./(h_n*h_n), 1./(h_n*h_n), dir_a, dir_n, dir_cross); + + //metr = intersection_conserveM1_bis(metr, curvMetric); //metr = intersection_conserveM1(metr,curvMetric); //metr = intersection_conserve_mostaniso(metr, curvMetric); - metr = intersection(metr,curvMetric); + //metr = intersection(metr,curvMetric); } } diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp index 643ee81d298c5ae8672c4d9ef7979daf16ebbfd5..4ae56efedea1c5a3fd769f51a9dcfe1f099f07c7 100644 --- a/Mesh/meshGRegionMMG3D.cpp +++ b/Mesh/meshGRegionMMG3D.cpp @@ -226,7 +226,7 @@ static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, if (itv != LCS.end()){ double LL = itv->second.first/itv->second.second; SMetric3 l4(1./(LL*LL)); - // printf("adding a size %g\n",LL); + //printf("adding a size %g\n",LL); SMetric3 MM = intersection_conserve_mostaniso (l4, m); m = MM; }