From 77b6b3cffbdc70a4baef6526aa50d14be7b38415 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Wed, 22 Aug 2012 11:46:02 +0000
Subject: [PATCH] Fixed 2D curvature calculation for "Eigendirections" mesh
metric (3D still to be implemented)
---
Mesh/meshMetric.cpp | 81 ++++++++++++++++++---------------------------
1 file changed, 32 insertions(+), 49 deletions(-)
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 12e9f7d006..e52fd40b5c 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -482,62 +482,45 @@ void meshMetric::computeMetricEigenDir()
hessian(1,0) = gradudy(0); hessian(1,1) = gradudy(1); hessian(1,2) = gradudy(2);
hessian(2,0) = gradudz(0); hessian(2,1) = gradudz(1); hessian(2,2) = gradudz(2);
- double metric_value_hmax = 1./hmax/hmax;
- SVector3 gr = grads[ver];
- double norm = gr.normalize();
+ const double metric_value_hmax = 1./(hmax*hmax);
+ const SVector3 gVec = grads[ver]; // Gradient vector
+ const double norm = gVec.norm();
SMetric3 metric;
if (signed_dist < _E && signed_dist > _E_moins && norm != 0.0){
- // first, compute the eigenvectors of the hessian, for a curvature-based metric
- fullMatrix<double> eigvec_hessian(3,3);
- fullVector<double> lambda_hessian(3);
- hessian.eig(eigvec_hessian,lambda_hessian);
- std::vector<SVector3> ti;
- ti.push_back(SVector3(eigvec_hessian(0,0),eigvec_hessian(1,0),eigvec_hessian(2,0)));
- ti.push_back(SVector3(eigvec_hessian(0,1),eigvec_hessian(1,1),eigvec_hessian(2,1)));
- ti.push_back(SVector3(eigvec_hessian(0,2),eigvec_hessian(1,2),eigvec_hessian(2,2)));
-
- // compute the required characteristic length relative to the curvature and the distance to iso
- double maxkappa = std::max(std::max(fabs(lambda_hessian(0)),fabs(lambda_hessian(1))),fabs(lambda_hessian(2)));
- // epsilon is set such that the characteristic element size in the tangent direction is h = 2pi R_min/_Np = sqrt(1/metric_maxmax) = sqrt(epsilon/kappa_max)
- double un_sur_epsilon = maxkappa/((2.*3.14/_Np)*(2.*3.14/_Np));
- // metric curvature-based eigenvalues
- std::vector<double> eigenvals_curvature;
- eigenvals_curvature.push_back(fabs(lambda_hessian(0))*un_sur_epsilon);
- eigenvals_curvature.push_back(fabs(lambda_hessian(1))*un_sur_epsilon);
- eigenvals_curvature.push_back(fabs(lambda_hessian(2))*un_sur_epsilon);
- // metric distance-based eigenvalue, in gradient direction
- double metric_value_hmin = 1./hmin/hmin;
- double eigenval_direction;
+ const double metric_value_hmin = 1./(hmin*hmin);
+ const SVector3 nVec = (1./norm)*gVec; // Unit normal vector
+ double lambda_n; // Eigenvalues of metric for normal & tangential directions
if (_technique==meshMetric::EIGENDIRECTIONS_LINEARINTERP_H){
- double h_dist = hmin + ((hmax-hmin)/_E)*dist;// the charcteristic element size in the normal direction - linear interp between hmin and hmax
- eigenval_direction = 1./h_dist/h_dist;
+ const double h_dist = hmin + ((hmax-hmin)/_E)*dist; // Characteristic element size in the normal direction - linear interp between hmin and hmax
+ lambda_n = 1./(h_dist*h_dist);
}
else if(_technique==meshMetric::EIGENDIRECTIONS){
- // ... or linear interpolation between 1/h_min^2 and 1/h_max^2
- double maximum_distance = (signed_dist>0.) ? _E : fabs(_E_moins);
- eigenval_direction = metric_value_hmin + ((metric_value_hmax-metric_value_hmin)/maximum_distance)*dist;
+ const double maximum_distance = (signed_dist>0.) ? _E : fabs(_E_moins); // ... or linear interpolation between 1/h_min^2 and 1/h_max^2
+ lambda_n = metric_value_hmin + ((metric_value_hmax-metric_value_hmin)/maximum_distance)*dist;
}
-
- // now, consider only three eigenvectors (i.e. (grad(LS), u1, u2) with u1 and u2 orthogonal vectors in the tangent plane to the LS)
- // and imposing directly eigenvalues and directions, instead of intersecting the two metrics based on direction and curvature
- // need to find out which one of ti's is the closest to gr, to known which eigenvalue to discard...
- std::vector<double> grad_dot_ti;
- for (int i=0;i<3;i++)
- grad_dot_ti.push_back(fabs(dot(gr,ti[i])));
- std::vector<double>::iterator itbegin = grad_dot_ti.begin();
- std::vector<double>::iterator itmax = std::max_element(itbegin,grad_dot_ti.end());
- int grad_index = std::distance(itbegin,itmax);
- std::vector<int> ti_index;
- for (int i=0;i<3;i++)
- if (i!=grad_index) ti_index.push_back(i);
- // finally, creating the metric
- std::vector<double> eigenvals;
- eigenvals.push_back(std::min(std::max(eigenval_direction,metric_value_hmax),metric_value_hmin));// in gradient direction
- eigenvals.push_back(std::min(std::max(eigenvals_curvature[ti_index[0]],metric_value_hmax),metric_value_hmin));
- eigenvals.push_back(std::min(std::max(eigenvals_curvature[ti_index[1]],metric_value_hmax),metric_value_hmin));
- metric = SMetric3(eigenvals[0],eigenvals[1],eigenvals[2],gr,ti[ti_index[0]],ti[ti_index[1]]);
- setOfSizes[metricNumber].insert(std::make_pair(ver, std::min(std::min(1/sqrt(eigenvals[0]),1/sqrt(eigenvals[1])),1/sqrt(eigenvals[2]))));
+ SVector3 tVec0, tVec1 = SPoint3(0.,0.,1.); // Unit tangential vectors
+ double kappa0, kappa1 = 0.; // Curvatures
+ if (_dim == 2) { // Curvature formulas: cf. R. Goldman, "Curvature formulas for implicit curves and surfaces", Computer Aided Geometric Design 22 (2005), pp. 632–658
+ tVec0 = SPoint3(-nVec(1),nVec(0),0.);
+ kappa0 = fabs(-gVec(1)*(-gVec(1)*hessian(0,0)+gVec(0)*hessian(0,1))+
+ gVec(0)*(-gVec(1)*hessian(1,0)+gVec(0)*hessian(1,1)))/pow(norm,3);
+ }
+ else {
+ // TODO: Implement 3D curvatures
+ tVec0 = SPoint3(0.,0.,0.);
+ tVec1 = SPoint3(0.,0.,0.);
+ kappa0 = 0.;
+ kappa1 = 0.;
+ }
+ const double invh_t0 = (_Np*kappa0)/6.283185307, invh_t1 = (_Np*kappa1)/6.283185307; // Inverse of tangential size 0
+ const double lambda_t0 = invh_t0*invh_t0, lambda_t1 = invh_t1*invh_t1;
+ const double lambdaP_n = std::min(std::max(lambda_n,metric_value_hmax),metric_value_hmin); // Truncate eigenvalues
+ const double lambdaP_t0 = std::min(std::max(lambda_t0,metric_value_hmax),metric_value_hmin);
+ const double lambdaP_t1 = (_dim == 2) ? 1. : std::min(std::max(lambda_t1,metric_value_hmax),metric_value_hmin);
+ metric = SMetric3(lambdaP_n,lambdaP_t0,lambdaP_t1,nVec,tVec0,tVec1);
+ const double h_n = 1./sqrt(lambdaP_n), h_t0 = 1./sqrt(lambdaP_t0), h_t1 = 1./sqrt(lambdaP_t1);
+ setOfSizes[metricNumber].insert(std::make_pair(ver,std::min(std::min(h_n,h_t0),h_t1)));
}
else{// isotropic metric !
SMetric3 mymetric(metric_value_hmax);
--
GitLab