diff --git a/Mesh/qualityMeasuresJacobian.cpp b/Mesh/qualityMeasuresJacobian.cpp
index 48ed8d53fbab83b89dfbbc76e8f452a7e045144c..d48d3206393c11771d9d9cb73de9eee98d131bab 100644
--- a/Mesh/qualityMeasuresJacobian.cpp
+++ b/Mesh/qualityMeasuresJacobian.cpp
@@ -37,7 +37,8 @@ double _CoeffDataScaledJac::cTri = 2/std::sqrt(3);
double _CoeffDataScaledJac::cTet = std::sqrt(2);
double _CoeffDataScaledJac::cPyr = std::sqrt(2)*4;
-void minMaxJacobianDeterminant(MElement *el, double &min, double &max)
+void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
+ const fullMatrix<double> *normals)
{
_CoeffDataAnisotropy::currentElement = el;
const JacobianBasis *jfs = el->getJacobianFuncSpace();
@@ -53,7 +54,7 @@ void minMaxJacobianDeterminant(MElement *el, double &min, double &max)
fullVector<double> coeffLag(jfs->getNumJacNodes());
fullVector<double> coeffBez(jfs->getNumJacNodes());
- jfs->getSignedJacobian(nodesXYZ, coeffLag);
+ jfs->getSignedJacobian(nodesXYZ, coeffLag, normals);
jfs->lag2Bez(coeffLag, coeffBez);
std::vector<_CoeffData*> domains;
diff --git a/Mesh/qualityMeasuresJacobian.h b/Mesh/qualityMeasuresJacobian.h
index dbf4eae0bfe72fd499cf9084373a972d22435367..fcc7d52a92e3274b8d7ec60246aaa9183e8736a3 100644
--- a/Mesh/qualityMeasuresJacobian.h
+++ b/Mesh/qualityMeasuresJacobian.h
@@ -15,7 +15,8 @@ class MElement;
namespace jacobianBasedQuality {
-void minMaxJacobianDeterminant(MElement *el, double &min, double &max);
+void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
+ const fullMatrix<double> *normals = NULL);
double minScaledJacobian(MElement *el,
bool knownValid = false,
bool reversedOk = false);
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index dc04a7bc72a44d4acce916bbcb06cba46d9f19a4..3bb6ac9d07badc028a600fdda16708a7d0a0c1be 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -446,7 +446,8 @@ void JacobianBasis::getJacobianGeneral(int nJacNodes,
const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ,
bool idealNorm, bool scaling,
- fullVector<double> &jacobian) const
+ fullVector<double> &jacobian,
+ const fullMatrix<double> *geomNormals) const
{
switch (_dim) {
@@ -458,11 +459,13 @@ void JacobianBasis::getJacobianGeneral(int nJacNodes,
case 1 : {
fullMatrix<double> normals(2,3);
const double invScale = getPrimNormals1D(nodesXYZ,normals);
- if (invScale == 0) {
- for (int i = 0; i < nJacNodes; i++) jacobian(i) = 0;
- return;
- }
+ if (geomNormals)
+ normals.setAll(*geomNormals);
if (scaling) {
+ if (invScale == 0) {
+ for (int i = 0; i < nJacNodes; i++) jacobian(i) = 0;
+ return;
+ }
const double scale = 1./invScale;
normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
}
@@ -482,6 +485,8 @@ void JacobianBasis::getJacobianGeneral(int nJacNodes,
case 2 : {
fullMatrix<double> normal(1,3);
const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
+ if (geomNormals)
+ normal.setAll(*geomNormals);
if (scaling) {
if (invScale == 0) {
for (int i = 0; i < nJacNodes; i++) jacobian(i) = 0;
@@ -541,7 +546,8 @@ void JacobianBasis::getJacobianGeneral(int nJacNodes,
const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ,
bool idealNorm, bool scaling,
- fullMatrix<double> &jacobian) const
+ fullMatrix<double> &jacobian,
+ const fullMatrix<double> *geomNormals) const
{
switch (_dim) {
@@ -565,11 +571,13 @@ void JacobianBasis::getJacobianGeneral(int nJacNodes,
}
fullMatrix<double> normals(2,3);
const double invScale = getPrimNormals1D(nodesXYZ,normals);
- if (invScale == 0) {
- for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 0;
- continue;
- }
+ if (geomNormals)
+ normals.setAll(*geomNormals);
if (scaling) {
+ if (invScale == 0) {
+ for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 0;
+ continue;
+ }
const double scale = 1./invScale;
normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
}
@@ -598,11 +606,13 @@ void JacobianBasis::getJacobianGeneral(int nJacNodes,
}
fullMatrix<double> normal(1,3);
const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
- if (invScale == 0) {
- for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 0;
- continue;
- }
+ if (geomNormals)
+ normal.setAll(*geomNormals);
if (scaling) {
+ if (invScale == 0) {
+ for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 0;
+ continue;
+ }
const double scale = 1./invScale;
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 07e797e04521bf1fa3053682a32637e3e4802224..25d14b0be30d357743384711f27602cef301a056 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -123,43 +123,47 @@ public:
fullVector<double> &lambdaJ,
fullMatrix<double> &gradLambdaJ) const;
inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ,
- fullVector<double> &jacobian) const {
+ fullVector<double> &jacobian,
+ const fullMatrix<double> *normals = NULL) const {
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
- nodesXYZ, false, false, jacobian);
+ nodesXYZ, false, false, jacobian, normals);
}
inline void getSignedJacobian(const fullMatrix<double> &nodesX,
const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ,
- fullMatrix<double> &jacobian) const {
+ fullMatrix<double> &jacobian,
+ const fullMatrix<double> *normals = NULL) const {
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
- nodesX, nodesY, nodesZ, false, false, jacobian);
+ nodesX, nodesY, nodesZ, false, false, jacobian, normals);
}
inline void getSignedIdealJacobian(const fullMatrix<double> &nodesXYZ,
- fullVector<double> &jacobian) const {
+ fullVector<double> &jacobian,
+ const fullMatrix<double> *normals = NULL) const {
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
_gradBasis->gradShapeIdealMatY,
_gradBasis->gradShapeIdealMatZ,
- nodesXYZ, true, false, jacobian);
+ nodesXYZ, true, false, jacobian, normals);
}
inline void getSignedIdealJacobian(const fullMatrix<double> &nodesX,
const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ,
- fullMatrix<double> &jacobian) const {
+ fullMatrix<double> &jacobian,
+ const fullMatrix<double> *normals = NULL) const {
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
_gradBasis->gradShapeIdealMatY,
_gradBasis->gradShapeIdealMatZ,
- nodesX, nodesY, nodesZ, true, false, jacobian);
+ nodesX, nodesY, nodesZ, true, false, jacobian, normals);
}
inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ,
fullVector<double> &jacobian) const {
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
- nodesXYZ, false, true, jacobian);
+ nodesXYZ, false, true, jacobian, NULL);
}
inline void getScaledJacobian(const fullMatrix<double> &nodesX,
const fullMatrix<double> &nodesY,
@@ -168,19 +172,21 @@ public:
getJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY,
_gradBasis->gradShapeMatZ,
- nodesX, nodesY, nodesZ, false, true, jacobian);
+ nodesX, nodesY, nodesZ, false, true, jacobian, NULL);
}
inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ,
- fullVector<double> &jacobian) const {
+ fullVector<double> &jacobian,
+ const fullMatrix<double> *normals = NULL) const {
getJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
gradShapeMatYFast, gradShapeMatZFast,
- nodesXYZ, false, false, jacobian);
+ nodesXYZ, false, false, jacobian, normals);
}
inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ,
- fullVector<double> &jacobian) const {
+ fullVector<double> &jacobian,
+ const fullMatrix<double> *normals = NULL) const {
getJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
gradShapeMatYFast, gradShapeMatZFast,
- nodesXYZ, false, true, jacobian);
+ nodesXYZ, false, true, jacobian, normals);
}
void lag2Bez(const fullVector<double> &lag, fullVector<double> &bez) const;
@@ -206,7 +212,8 @@ public:
const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ,
bool idealNorm, bool scaling,
- fullVector<double> &jacobian) const;
+ fullVector<double> &jacobian,
+ const fullMatrix<double> *normals) const;
void getJacobianGeneral(int nJacNodes,
const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY,
@@ -215,7 +222,8 @@ public:
const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ,
bool idealNorm, bool scaling,
- fullMatrix<double> &jacobian) const;
+ fullMatrix<double> &jacobian,
+ const fullMatrix<double> *normals) const;
void getSignedJacAndGradientsGeneral(int nJacNodes,
const fullMatrix<double> &gSMatX,
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 6e8eb1bfc71aa7353816f506f0bee9b4b6fff525..b2c14347c893a9f35031e1daba804a7e6b399ce9 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -248,6 +248,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
for (it = entities.begin(); it != entities.end(); ++it) {
GEntity *entity = *it;
unsigned num = entity->getNumMeshElements();
+ fullMatrix<double> *normals = NULL;
switch (dim) {
case 3:
Msg::StatusBar(true, "Volume %d: checking the Jacobian of %d elements",
@@ -256,10 +257,69 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
case 2:
Msg::StatusBar(true, "Surface %d: checking the Jacobian of %d elements",
entity->tag(), num);
+ if (entity->geomType() == GEntity::Plane &&
+ entity->haveParametrization()) {
+ double u = entity->parBounds(0).low();
+ double v = entity->parBounds(1).low();
+ SVector3 n = dynamic_cast<GFace*>(entity)->normal(SPoint2(u, v));
+ normals = new fullMatrix<double>(1, 3);
+ normals->set(0, 0, n[0]);
+ normals->set(0, 1, n[1]);
+ normals->set(0, 2, n[2]);
+ }
+ else if (entity->geomType() == GEntity::DiscreteSurface) {
+ SBoundingBox3d bb = entity->bounds();
+ // If we don't have the CAD, check if the mesh is 2D:
+ if (!bb.empty() && bb.max().z() - bb.min().z() == .0) {
+ normals = new fullMatrix<double>(1, 3);
+ normals->set(0, 0, 0);
+ normals->set(0, 1, 0);
+ normals->set(0, 2, 1);
+ Msg::Error("Discrete face is 2D!");
+ }
+ }
break;
case 1:
Msg::StatusBar(true, "Line %d: checking the Jacobian of %d elements",
entity->tag(), num);
+ if (entity->geomType() == GEntity::Line &&
+ entity->haveParametrization()) {
+ double u = entity->parBounds(0).low();
+ SVector3 t = dynamic_cast<GEdge*>(entity)->firstDer(u);
+ SVector3 dum = SVector3(0, 0, 0);
+ if(t[0] == 0.)
+ dum[0] = 1;
+ else if(t[1] == 0.)
+ dum[1] = 1;
+ else
+ dum[2] = 1;
+ SVector3 n1, n2;
+ n1 = crossprod(t, dum);
+ n1.normalize();
+ n2 = crossprod(t, n1);
+ n2.normalize();
+ normals = new fullMatrix<double>(2, 3);
+ normals->set(0, 0, n1[0]);
+ normals->set(0, 1, n1[1]);
+ normals->set(0, 2, n1[2]);
+ normals->set(1, 0, n2[0]);
+ normals->set(1, 1, n2[1]);
+ normals->set(1, 2, n2[2]);
+ }
+ else if (entity->geomType() == GEntity::DiscreteCurve) {
+ SBoundingBox3d bb = entity->bounds();
+ // If we don't have the CAD, check if the mesh is 1D:
+ if (!bb.empty() && bb.max().y() - bb.min().y() == .0 &&
+ bb.max().z() - bb.min().z() == .0) {
+ normals = new fullMatrix<double>(2, 3);
+ normals->set(0, 0, 0);
+ normals->set(0, 1, 1);
+ normals->set(0, 2, 0);
+ normals->set(1, 0, 0);
+ normals->set(1, 1, 0);
+ normals->set(1, 2, 1);
+ }
+ }
break;
default: break;
}
@@ -270,7 +330,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
for (unsigned i = 0; i < num; ++i) {
MElement *el = entity->getMeshElement(i);
double min, max;
- jacobianBasedQuality::minMaxJacobianDeterminant(el, min, max);
+ jacobianBasedQuality::minMaxJacobianDeterminant(el, min, max, normals);
_data.push_back(data_elementMinMax(el, min, max));
if (min < 0 && max < 0) {
Msg::Warning("Element %d is completely inverted", el->getNum());
@@ -297,6 +357,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(int dim)
}
}
}
+ delete normals;
}
_computedJ[dim-1] = true;
}