From 52ecae224ce07c05015238dbf14c5f482861526b Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Sat, 4 May 2013 19:19:35 +0000
Subject: [PATCH] Changed "most aniso" mesh metric intersection
---
Geo/STensor3.cpp | 62 ++++++++++++++++++++++++++++++++++++++++++---
Geo/STensor3.h | 1 +
Mesh/meshMetric.cpp | 43 +++++++++++++++++++++----------
3 files changed, 89 insertions(+), 17 deletions(-)
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 4a61ce2f6b..a282ad4f04 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -64,22 +64,76 @@ SMetric3 intersection_conserveM1 (const SMetric3 &m1, const SMetric3 &m2)
return iv;
}
+namespace {
+
+double anisoRatio2D (double v0x, double v0y, double v0z,
+ double v1x, double v1y, double v1z,
+ double v2x, double v2y, double v2z,
+ double s0, double s1, double s2) {
+
+ static const double eps = 1.e-8;
+
+ double ratio;
+ if (v0x*v0x+v0y*v0y+(v0z-1.)*(v0z-1.) < eps) // If 1st eigenvect corresponds to z dir. ...
+ ratio = std::min(fabs(s1/s2),fabs(s2/s1)); // ... consider ratio of 2nd and 3rd eigenval.
+ else if (v1x*v1x+v1y*v1y+(v1z-1.)*(v1z-1.) < eps) // If 2nd eigenvect corresponds to z dir. ...
+ ratio = std::min(fabs(s0/s2),fabs(s2/s0)); // ... consider ratio of 1st and 3rd eigenval.
+ else if (v2x*v2x+v2y*v2y+(v2z-1.)*(v2z-1.) < eps) // If 3rd eigenvect corresponds to z dir. ...
+ ratio = std::min(fabs(s0/s1),fabs(s1/s0)); // ... consider ratio of 1st and 2nd eigenval.
+ else {
+ printf("Error in anisoRatio2D: z direction not found!\n");
+ ratio = 0.;
+ }
+
+ return ratio;
+
+}
+
+}
+
// preserve orientation of the most anisotropic metric !!!
SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2)
{
fullMatrix<double> V1(3,3);
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 ratio1 = fabs(S1(0)/S1(2)); // Minimum ratio because we take sorted eigenvalues
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)));
-
- if (lambda1_min < lambda2_min)
+ double ratio2 = fabs(S2(0)/S2(2)); // Minimum ratio because we take sorted eigenvalues
+
+ if (ratio1 < ratio2)
+ return intersection_conserveM1(m1, m2);
+ else
+ return intersection_conserveM1(m2, m1);
+}
+
+// preserve orientation of the most anisotropic metric in 2D!!!
+SMetric3 intersection_conserve_mostaniso_2d (const SMetric3 &m1, const SMetric3 &m2)
+{
+
+ fullMatrix<double> V1(3,3);
+ fullVector<double> S1(3);
+ m1.eig(V1,S1,false);
+ double ratio1 = anisoRatio2D(V1(0,0),V1(1,0),V1(2,0),
+ V1(0,1),V1(1,1),V1(2,1),
+ V1(0,2),V1(1,2),V1(2,2),
+ S1(0),S1(1),S1(2));
+
+ fullMatrix<double> V2(3,3);
+ fullVector<double> S2(3);
+ m2.eig(V2,S2,false);
+ double ratio2 = anisoRatio2D(V2(0,0),V2(1,0),V2(2,0),
+ V2(0,1),V2(1,1),V2(2,1),
+ V2(0,2),V2(1,2),V2(2,2),
+ S2(0),S2(1),S2(2));
+
+ if (ratio1 < ratio2)
return intersection_conserveM1(m1, m2);
else
return intersection_conserveM1(m2, m1);
+
}
// (1-t) * m1 + t * m2
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index d54a2b8c95..cb427ab811 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -169,6 +169,7 @@ SMetric3 intersection_conserveM1 (const SMetric3 &m1,
// preserve orientation of the most anisotropic metric
SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2);
+SMetric3 intersection_conserve_mostaniso_2d (const SMetric3 &m1, const SMetric3 &m2);
// compute the largest inscribed ellipsoid...
SMetric3 intersection (const SMetric3 &m1,
const SMetric3 &m2);
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index f07ac65360..4a41e97ec2 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -95,7 +95,9 @@ void meshMetric::updateMetrics()
if (setOfMetrics.size() > 1)
for (unsigned int i=1;i<setOfMetrics.size();i++){
- _nodalMetrics[ver] = intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]);
+ _nodalMetrics[ver] = (_dim == 3) ?
+ intersection_conserve_mostaniso(_nodalMetrics[ver],setOfMetrics[i][ver]) :
+ intersection_conserve_mostaniso_2d(_nodalMetrics[ver],setOfMetrics[i][ver]);
_nodalSizes[ver] = std::min(_nodalSizes[ver],setOfSizes[i][ver]);
}
}
@@ -105,19 +107,22 @@ void meshMetric::updateMetrics()
void meshMetric::exportInfo(const char * fileendname)
{
if (needMetricUpdate) updateMetrics();
- std::stringstream sg,sm,sl,sh;
+ std::stringstream sg,sm,sl,sh,shm;
sg << "meshmetric_gradients_" << fileendname;
sm << "meshmetric_metric_" << fileendname;
sl << "meshmetric_levelset_" << fileendname;
sh << "meshmetric_hessian_" << fileendname;
+ shm << "meshmetric_hessianmat_" << fileendname;
std::ofstream out_grad(sg.str().c_str());
std::ofstream out_metric(sm.str().c_str());
std::ofstream out_ls(sl.str().c_str());
std::ofstream out_hess(sh.str().c_str());
+ std::ofstream out_hessmat(shm.str().c_str());
out_grad << "View \"ls_gradient\"{" << std::endl;
out_metric << "View \"metric\"{" << std::endl;
out_ls << "View \"ls\"{" << std::endl;
out_hess << "View \"hessian\"{" << std::endl;
+ out_hessmat << "View \"hessian_mat\"{" << std::endl;
std::vector<MElement*>::iterator itelem = _elements.begin();
std::vector<MElement*>::iterator itelemen = _elements.end();
for (;itelem!=itelemen;itelem++){
@@ -127,12 +132,14 @@ void meshMetric::exportInfo(const char * fileendname)
out_grad << "VT(";
out_ls << "ST(";
out_hess << "ST(";
+ out_hessmat << "TT(";
}
else {
out_metric << "TS(";
out_grad << "VS(";
out_ls << "SS(";
out_hess << "SS(";
+ out_hessmat << "TS(";
}
for ( int i = 0; i < e->getNumVertices(); i++) {
MVertex *ver = e->getVertex(i);
@@ -140,17 +147,20 @@ void meshMetric::exportInfo(const char * fileendname)
out_grad << ver->x() << "," << ver->y() << "," << ver->z();
out_ls << ver->x() << "," << ver->y() << "," << ver->z();
out_hess << ver->x() << "," << ver->y() << "," << ver->z();
+ out_hessmat << ver->x() << "," << ver->y() << "," << ver->z();
if (i!=e->getNumVertices()-1){
out_metric << ",";
out_grad << ",";
out_ls << ",";
out_hess << ",";
+ out_hessmat << ",";
}
else{
out_metric << "){";
out_grad << "){";
out_ls << "){";
- out_hess << "){";
+ out_hess << "){";
+ out_hessmat << "){";
}
}
for ( int i = 0; i < e->getNumVertices(); i++) {
@@ -161,12 +171,12 @@ void meshMetric::exportInfo(const char * fileendname)
SVector3 gradudz = dgrads[2][ver];
out_hess << (gradudx(0) + gradudy(1) + gradudz(2));
if (i == (e->getNumVertices() - 1)){
- out_ls << "};" << std::endl;
- out_hess << "};" << std::endl;
+ out_ls << "};" << std::endl;
+ out_hess << "};" << std::endl;
}
else {
- out_ls << ",";
- out_hess << ",";
+ out_ls << ",";
+ out_hess << ",";
}
for (int k=0;k<3;k++){
out_grad << grads[ver](k);
@@ -174,8 +184,15 @@ void meshMetric::exportInfo(const char * fileendname)
else out_grad << ",";
for (int l=0;l<3;l++){
out_metric << _nodalMetrics[ver](k,l);
- if ((k == 2) && (l == 2) && (i == (e->getNumVertices() - 1))) out_metric << "};" << std::endl;
- else out_metric << ",";
+ out_hessmat << dgrads[k][ver](l);
+ if ((k == 2) && (l == 2) && (i == (e->getNumVertices() - 1))) {
+ out_metric << "};" << std::endl;
+ out_hessmat << "};" << std::endl;
+ }
+ else {
+ out_metric << ",";
+ out_hessmat << ",";
+ }
}
}
}
@@ -184,10 +201,12 @@ void meshMetric::exportInfo(const char * fileendname)
out_metric << "};" << std::endl;
out_ls << "};" << std::endl;
out_hess << "};" << std::endl;
+ out_hessmat << "};" << std::endl;
out_grad.close();
out_metric.close();
out_ls.close();
out_hess.close();
+ out_hessmat.close();
}
meshMetric::~meshMetric()
@@ -360,10 +379,9 @@ void meshMetric::computeMetricLevelSet(MVertex *ver, SMetric3 &hessian, SMetric
void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian, SMetric3 &metric, double &size, double x, double y, double z)
{
- double signed_dist;
+
SVector3 gradudx, gradudy,gradudz,gr;
if(ver != NULL){
- signed_dist = vals[ver];
gr = grads[ver];
gradudx = dgrads[0][ver];
gradudy = dgrads[1][ver];
@@ -373,7 +391,6 @@ void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian, SMetric3
hessian(2,0) = gradudz(0); hessian(2,1) = gradudz(1); hessian(2,2) = gradudz(2);
}
else if (ver == NULL){
- signed_dist = (*_fct)(x,y,z);
_fct->gradient(x,y,z,gr(0),gr(1),gr(2));
_fct->hessian(x,y,z, hessian(0,0),hessian(0,1),hessian(0,2),
hessian(1,0),hessian(1,1),hessian(1,2),
@@ -394,7 +411,7 @@ void meshMetric::computeMetricHessian(MVertex *ver, SMetric3 &hessian, SMetric3
SVector3 t1 (V(0,0),V(1,0),V(2,0));
SVector3 t2 (V(0,1),V(1,1),V(2,1));
- SVector3 t3 (V(0,2),V(1,2),V(2,2));
+ SVector3 t3 = (_dim == 3) ? SVector3(V(0,2),V(1,2),V(2,2)) : SVector3(0.,0.,1.);
size = std::min(std::min(1/sqrt(lambda1),1/sqrt(lambda2)),1/sqrt(lambda3));
metric = SMetric3(lambda1,lambda2,lambda3,t1,t2,t3);
--
GitLab