From 4e09e144575cd542c4584612848fe96ed5984406 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Wed, 28 Nov 2012 13:16:17 +0000
Subject: [PATCH] Fixed generic 3D distortion computation
---
Geo/MElement.cpp | 18 ++++++++----
Mesh/HighOrder.cpp | 4 +--
Mesh/qualityMeasures.cpp | 60 ++++++++++++++++++++--------------------
3 files changed, 44 insertions(+), 38 deletions(-)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index d5a221c8bd..5f0d00f67f 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -112,19 +112,25 @@ void MElement::scaledJacRange(double &jmin, double &jmax)
{
jmin = jmax = 1.0;
#if defined(HAVE_MESH)
-// if (getType() == TYPE_PYR) return;
- extern double mesh_functional_distorsion(MElement*,double,double);
+ extern double mesh_functional_distorsion_2D(MElement*,double,double);
+ extern double mesh_functional_distorsion_3D(MElement*,double,double,double);
if (getPolynomialOrder() == 1) return;
const bezierBasis *jac = getJacobianFuncSpace()->bezier;
fullVector<double>Ji(jac->points.size1());
for (int i=0;i<jac->points.size1();i++){
double u = jac->points(i,0);
double v = jac->points(i,1);
- if (getType() == TYPE_QUA){
- u = -1 + 2*u;
- v = -1 + 2*v;
+ if (getDim() == 2) {
+ if (getType() == TYPE_QUA) {
+ u = -1 + 2*u;
+ v = -1 + 2*v;
+ }
+ Ji(i) = mesh_functional_distorsion_2D(this,u,v);
+ }
+ else {
+ double w = jac->points(i,2);
+ Ji(i) = mesh_functional_distorsion_3D(this,u,v,w);
}
- Ji(i) = mesh_functional_distorsion(this,u,v);
}
fullVector<double> Bi( jac->matrixLag2Bez.size1() );
jac->matrixLag2Bez.mult(Ji,Bi);
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 8439e342d9..07e7077804 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1325,7 +1325,7 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m,
avg / (count ? count : 1));
}
-extern double mesh_functional_distorsion(MElement *t, double u, double v);
+extern double mesh_functional_distorsion_2D(MElement *t, double u, double v);
void printJacobians(GModel *m, const char *nm)
{
@@ -1343,7 +1343,7 @@ void printJacobians(GModel *m, const char *nm)
double u = (double)i / (n - 1);
double v = (double)k / (n - 1);
t->pnt(u, v, 0, pt);
- D[i][k] = mesh_functional_distorsion(t, u, v);
+ D[i][k] = mesh_functional_distorsion_2D(t, u, v);
//X[i][k] = u;
//Y[i][k] = v;
//Z[i][k] = 0.0;
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 4b1a646757..37b5f13bf4 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -208,7 +208,7 @@ double conditionNumberAndDerivativeOfTet(const double &x1, const double &y1, con
}
*/
-double mesh_functional_distorsion(MElement *t, double u, double v)
+double mesh_functional_distorsion_2D(MElement *t, double u, double v)
{
// compute uncurved element jacobian d_u x and d_v x
double mat[3][3];
@@ -253,25 +253,25 @@ static double MINQ (double a, double b, double c){
double mesh_functional_distorsion_p2_bezier_refined(MTriangle *t)
{
- double J1 =mesh_functional_distorsion(t,0.0,0.0);
- double J2 =mesh_functional_distorsion(t,1.0,0.0);
- double J3 =mesh_functional_distorsion(t,0.0,1.0);
- double J4 =mesh_functional_distorsion(t,0.5,0.0);
- double J5 =mesh_functional_distorsion(t,0.5,0.5);
- double J6 =mesh_functional_distorsion(t,0.0,0.5);
+ double J1 =mesh_functional_distorsion_2D(t,0.0,0.0);
+ double J2 =mesh_functional_distorsion_2D(t,1.0,0.0);
+ double J3 =mesh_functional_distorsion_2D(t,0.0,1.0);
+ double J4 =mesh_functional_distorsion_2D(t,0.5,0.0);
+ double J5 =mesh_functional_distorsion_2D(t,0.5,0.5);
+ double J6 =mesh_functional_distorsion_2D(t,0.0,0.5);
- double J36 =mesh_functional_distorsion(t,0.0,.75);
- double J35 =mesh_functional_distorsion(t,0.25,.75);
- double J56 =mesh_functional_distorsion(t,0.25,.5);
+ double J36 =mesh_functional_distorsion_2D(t,0.0,.75);
+ double J35 =mesh_functional_distorsion_2D(t,0.25,.75);
+ double J56 =mesh_functional_distorsion_2D(t,0.25,.5);
- double J16 =mesh_functional_distorsion(t,0.0,.25);
- double J14 =mesh_functional_distorsion(t,0.25,.0);
- double J46 =mesh_functional_distorsion(t,0.25,.25);
+ double J16 =mesh_functional_distorsion_2D(t,0.0,.25);
+ double J14 =mesh_functional_distorsion_2D(t,0.25,.0);
+ double J46 =mesh_functional_distorsion_2D(t,0.25,.25);
- double J45 =mesh_functional_distorsion(t,0.5,.25);
- double J52 =mesh_functional_distorsion(t,0.75,.25);
- double J24 =mesh_functional_distorsion(t,0.75,.0);
+ double J45 =mesh_functional_distorsion_2D(t,0.5,.25);
+ double J52 =mesh_functional_distorsion_2D(t,0.75,.25);
+ double J24 =mesh_functional_distorsion_2D(t,0.75,.0);
double d[15] = {
J1,J6,J4,2*J16-0.5*(J1+J6),2*J14-0.5*(J1+J4),2*J46-0.5*(J6+J4),
@@ -282,12 +282,12 @@ double mesh_functional_distorsion_p2_bezier_refined(MTriangle *t)
double mesh_functional_distorsion_p2_exact(MTriangle *t)
{
- double J1 =mesh_functional_distorsion(t,0.0,0.0);
- double J2 =mesh_functional_distorsion(t,1.0,0.0);
- double J3 =mesh_functional_distorsion(t,0.0,1.0);
- double J4 =mesh_functional_distorsion(t,0.5,0.0);
- double J5 =mesh_functional_distorsion(t,0.5,0.5);
- double J6 =mesh_functional_distorsion(t,0.0,0.5);
+ double J1 =mesh_functional_distorsion_2D(t,0.0,0.0);
+ double J2 =mesh_functional_distorsion_2D(t,1.0,0.0);
+ double J3 =mesh_functional_distorsion_2D(t,0.0,1.0);
+ double J4 =mesh_functional_distorsion_2D(t,0.5,0.0);
+ double J5 =mesh_functional_distorsion_2D(t,0.5,0.5);
+ double J6 =mesh_functional_distorsion_2D(t,0.0,0.5);
const double a = J1;
const double b = -3*J1-J2+4*J4;
@@ -347,7 +347,7 @@ double mesh_functional_distorsion_pN(MElement *t)
v = -1 + 2*v;
}
- Ji(i) = mesh_functional_distorsion(t,u,v);
+ Ji(i) = mesh_functional_distorsion_2D(t,u,v);
}
fullVector<double> Bi( jac->matrixLag2Bez.size1() );
@@ -406,16 +406,16 @@ double qmDistorsionOfMapping (MQuadrangle *e)
}
-static double mesh_functional_distorsion(MTetrahedron *t, double u, double v, double w)
+double mesh_functional_distorsion_3D(MElement *t, double u, double v, double w)
{
// compute uncurved element jacobian d_u x and d_v x
double mat[3][3];
- t->getPrimaryJacobian(u, v, w, mat);
- const double det1 = det3x3(mat);
- t->getJacobian(u, v, w, mat);
- const double detN = det3x3(mat);
+ const double det1 = t->getPrimaryJacobian(u, v, w, mat);
+// const double det1 = det3x3(mat);
+ const double detN = t->getJacobian(u, v, w, mat);
+// const double detN = det3x3(mat);
if (det1 == 0 || detN == 0) return 0;
- double dist = det1 / detN;
+ double dist = detN / det1;
return dist;
}
@@ -427,7 +427,7 @@ double qmDistorsionOfMapping(MTetrahedron *t)
const double u = jac->points(i,0);
const double v = jac->points(i,1);
const double w = jac->points(i,2);
- Ji(i) = mesh_functional_distorsion(t,u,v,w);
+ Ji(i) = mesh_functional_distorsion_3D(t,u,v,w);
}
fullVector<double> Bi( jac->matrixLag2Bez.size1() );
--
GitLab