diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 13320799a460b31a27b10ec0e42a206d684f0434..cba205218dcf19dc938960d3d9494cd8457279a0 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -277,27 +277,33 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
}
-namespace {
-
-// Calculate (signed) Jacobian for one element, given the gradients of the shape functions
-// and vectors for regularization of line and surface Jacobians in 3D.
-inline void computeSignedJacobian(int dim, int numJacNodes, const fullMatrix<double> &gradShapeMatX,
- const fullMatrix<double> &gradShapeMatY, const fullMatrix<double> &gradShapeMatZ,
- const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
- fullVector<double> &jacobian)
+// Calculate (signed, possibly scaled) Jacobian for one element, with normal vectors to straight element
+// for regularization. Evaluation points depend on the given matrices for shape function gradients.
+template<bool scaling>
+inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
+ const int dim = bezier->getDim();
+
switch (dim) {
case 0 : {
- for (int i = 0; i < numJacNodes; i++) jacobian(i) = 1.;
+ for (int i = 0; i < nJacNodes; i++) jacobian(i) = 1.;
break;
}
case 1 : {
- fullMatrix<double> dxyzdX(numJacNodes,3);
- gradShapeMatX.mult(nodesXYZ, dxyzdX);
- for (int i = 0; i < numJacNodes; i++) {
+ fullMatrix<double> normals(2,3);
+ const double invScale = getPrimNormals1D(nodesXYZ,normals);
+ if (scaling) {
+ const double scale = 1./invScale;
+ normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
+ }
+ fullMatrix<double> dxyzdX(nJacNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ for (int i = 0; i < nJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
@@ -307,29 +313,40 @@ inline void computeSignedJacobian(int dim, int numJacNodes, const fullMatrix<dou
}
case 2 : {
- fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3);
- gradShapeMatX.mult(nodesXYZ, dxyzdX);
- gradShapeMatY.mult(nodesXYZ, dxyzdY);
- for (int i = 0; i < numJacNodes; i++) {
+ fullMatrix<double> normal(1,3);
+ const double invScale = getPrimNormal2D(nodesXYZ,normal);
+ if (scaling) {
+ const double scale = 1./invScale;
+ normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
+ }
+ fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ for (int i = 0; i < nJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
- const double &dxdZ = normals(0,0), &dydZ = normals(0,1), &dzdZ = normals(0,2);
+ const double &dxdZ = normal(0,0), &dydZ = normal(0,1), &dzdZ = normal(0,2);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
}
break;
}
case 3 : {
- fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3);
- gradShapeMatX.mult(nodesXYZ, dxyzdX);
- gradShapeMatY.mult(nodesXYZ, dxyzdY);
- gradShapeMatZ.mult(nodesXYZ, dxyzdZ);
- for (int i = 0; i < numJacNodes; i++) {
+ fullMatrix<double> dum;
+ fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3), dxyzdZ(nJacNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ gSMatZ.mult(nodesXYZ, dxyzdZ);
+ for (int i = 0; i < nJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
}
+ if (scaling) {
+ const double scale = 1./getPrimJac3D(nodesXYZ);
+ jacobian.scale(scale);
+ }
break;
}
@@ -337,91 +354,32 @@ inline void computeSignedJacobian(int dim, int numJacNodes, const fullMatrix<dou
}
-} // namespace
-
-// Calculate (signed) Jacobian for one element, with normal vectors to straight element for regularization.
-// Evaluation points depend on the given matrices for shape function gradients.
+// Calculate signed Jacobian for one element, with normal vectors to straight element for
+// regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
-
- const int dim = bezier->getDim();
-
- switch (dim) {
-
- case 1 : {
- fullMatrix<double> normals(2,3);
- getPrimNormals1D(nodesXYZ,normals);
- computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
- break;
- }
-
- case 2 : {
- fullMatrix<double> normal(1,3);
- getPrimNormal2D(nodesXYZ,normal);
- computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
- break;
- }
-
- case 0 :
- case 3 : {
- fullMatrix<double> dum;
- computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
- break;
- }
-
- }
-
+ getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
}
-// Calculate scaled (signed) Jacobian for one element, with normal vectors to straight element for regularization
-// and scaling. Evaluation points depend on the given matrices for shape function gradients.
+// Calculate (signed) scaled Jacobian for one element, with normal vectors to straight element
+// for regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
-
- const int dim = bezier->getDim();
-
- switch (dim) {
-
- case 1 : {
- fullMatrix<double> normals(2,3);
- const double scale = 1./getPrimNormals1D(nodesXYZ,normals);
- normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
- computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
- break;
- }
-
- case 2 : {
- fullMatrix<double> normal(1,3);
- const double scale = 1./getPrimNormal2D(nodesXYZ,normal);
- normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
- computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
- break;
- }
-
- case 0 :
- case 3 : {
- fullMatrix<double> dum;
- const double scale = 1./getPrimJac3D(nodesXYZ);
- computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
- jacobian.scale(scale);
- break;
- }
-
- }
-
+ getJacobianGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
}
-// Calculate (signed) Jacobian for several elements.
-// Evaluation points depend on the given matrices for shape function gradients.
+// Calculate (signed, possibly scaled) Jacobian for several elements, with normal vectors to straight
+// elements for regularization. Evaluation points depend on the given matrices for shape function gradients.
// TODO: Optimize and test 1D & 2D
-void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
+template<bool scaling>
+inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
switch (bezier->getDim()) {
@@ -438,14 +396,18 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
for (int iEl = 0; iEl < numEl; iEl++) {
- fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
- for (int i = 0; i < numPrimJacNodes; i++) {
+ fullMatrix<double> nodesXYZ(numPrimMapNodes,3);
+ for (int i = 0; i < numPrimMapNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normals(2,3);
- getPrimNormals1D(nodesXYZ,normals);
+ const double invScale = getPrimNormals1D(nodesXYZ,normals);
+ if (scaling) {
+ const double scale = 1./invScale;
+ normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
+ }
const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
for (int i = 0; i < nJacNodes; i++)
@@ -463,14 +425,18 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
for (int iEl = 0; iEl < numEl; iEl++) {
- fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
- for (int i = 0; i < numPrimJacNodes; i++) {
+ fullMatrix<double> nodesXYZ(numPrimMapNodes,3);
+ for (int i = 0; i < numPrimMapNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normal(1,3);
- getPrimNormal2D(nodesXYZ,normal);
+ const double invScale = getPrimNormal2D(nodesXYZ,normal);
+ if (scaling) {
+ const double scale = 1./invScale;
+ normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
+ }
const double &dxdZ = normal(0,0), &dydZ = normal(0,1), &dzdZ = normal(0,2);
for (int i = 0; i < nJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
@@ -488,11 +454,22 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
gSMatZ.mult(nodesX, dxdZ); gSMatZ.mult(nodesY, dydZ); gSMatZ.mult(nodesZ, dzdZ);
- for (int iEl = 0; iEl < numEl; iEl++)
+ for (int iEl = 0; iEl < numEl; iEl++) {
for (int i = 0; i < nJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY(i,iEl),dydY(i,iEl),dzdY(i,iEl),
dxdZ(i,iEl),dydZ(i,iEl),dzdZ(i,iEl));
+ if (scaling) {
+ fullMatrix<double> nodesXYZ(numPrimMapNodes,3);
+ for (int i = 0; i < numPrimMapNodes; i++) {
+ nodesXYZ(i,0) = nodesX(i,iEl);
+ nodesXYZ(i,1) = nodesY(i,iEl);
+ nodesXYZ(i,2) = nodesZ(i,iEl);
+ }
+ const double scale = 1./getPrimJac3D(nodesXYZ);
+ for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) *= scale;
+ }
+ }
break;
}
@@ -500,6 +477,26 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
}
+// Calculate signed Jacobian for several elements, with normal vectors to straight elements for
+// regularization. Evaluation points depend on the given matrices for shape function gradients.
+void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
+{
+ getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
+}
+
+// Calculate (signed) scaled Jacobian for several elements, with normal vectors to straight elements
+// for regularization. Evaluation points depend on the given matrices for shape function gradients.
+void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
+{
+ getJacobianGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
+}
+
// Calculate (signed) Jacobian and its gradients for one element, with normal vectors to straight element
// for regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index db54b0c1021bd93718d9d29b32c60a0054db3d4c..cdcd1adab34bc8bed90e26ccce24618af7a0d65c 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -86,14 +86,14 @@ class JacobianBasis {
getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesXYZ,jacobian);
}
+ inline void getScaledJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
+ getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
+ _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
+ }
inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
}
- inline void getSignedJacobianFast(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
- getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
- gradShapeMatZFast,nodesX,nodesY,nodesZ,jacobian);
- }
inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
}
@@ -122,6 +122,15 @@ class JacobianBasis {
private :
+ template<bool scaling>
+ inline void getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
+ template<bool scaling>
+ inline void getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
@@ -132,6 +141,11 @@ class JacobianBasis {
void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
+ void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
+
void getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index f99be2624e06e9058e20412291bb58902fe50932..530d39f2193432abb36c33abf1291525b9a9c5c8 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -24,15 +24,36 @@
#endif
namespace {
- double sum(fullVector<double> &v)
- {
- double sum = .0;
- for (int i = 0; i < v.size(); ++i) {
- sum += v(i);
- }
- return sum;
- }
-}
+
+class BezierJacobian;
+struct lessMinB {
+ bool operator()(BezierJacobian*, BezierJacobian*) const;
+};
+struct lessMaxB {
+ bool operator()(BezierJacobian*, BezierJacobian*) const;
+};
+
+class BezierJacobian
+{
+private:
+ fullVector<double> _jacBez;
+ double _minJ, _maxJ, _minB, _maxB; //Extremum of Jac at corners and of bezier values
+ int _depthSub;
+ const JacobianBasis *_jfs;
+
+public:
+ BezierJacobian(fullVector<double> &, const JacobianBasis *, int depth);
+ void subDivisions(fullVector<double> &vect) const
+ {_jfs->subdivideBezierCoeff(_jacBez, vect);}
+
+ inline int depth() const {return _depthSub;}
+ inline double minJ() const {return _minJ;}
+ inline double maxJ() const {return _maxJ;}
+ inline double minB() const {return _minB;}
+ inline double maxB() const {return _maxB;}
+};
+
+} // namespace
//#define UNDEF_JAC_TAG -999
//#define _ANALYSECURVEDMESH_BLAS_
@@ -141,8 +162,6 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
break;
case 2 :
- Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
-
for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
GFace *f = *it;
@@ -159,7 +178,6 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
GFace *f = *it;
unsigned int numType[3] = {0, 0, 0};
f->getNumMeshElements(numType);
-
for (int type = 0; type < 3; type++) {
MElement *const *el = f->getStartElementType(type);
checkValidity(el, numType[type], invalids);
@@ -169,7 +187,6 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
break;
case 1 :
- Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
GEdge *e = *it;
unsigned int numElement = e->getNumMeshElements();
@@ -224,23 +241,20 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
if (numEl < 1)
return;
- const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1);
- const JacobianBasis *jfs1 = el[0]->getJacobianFuncSpace(1);
- if (!jfs || !jfs1) {
+ const JacobianBasis *jfs = el[0]->getJacobianFuncSpace();
+ if (!jfs) {
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
- const int numSamplingPt = jfs->getNumJacNodes(), numSamplingPt1 = jfs1->getNumJacNodes();
- const int numMapNodes = jfs->getNumMapNodes(), numMapNodes1 = jfs1->getNumMapNodes();
+ const int numSamplingPt = jfs->getNumJacNodes();
+ const int numMapNodes = jfs->getNumMapNodes();
#ifdef _ANALYSECURVEDMESH_BLAS_
fullMatrix<double> jacobianB(numSamplingPt, numEl);
fullMatrix<double> jacBezB(numSamplingPt, numEl);
- fullMatrix<double> jac1B(numSamplingPt1, numEl);
- fullVector<double> jacBez, jacobian, jac1;
+ fullVector<double> jacBez, jacobian;
fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
- fullMatrix<double> nodesX1(numMapNodes, numEl), nodesY1(numMapNodes, numEl), nodesZ1(numMapNodes, numEl);
for (int k = 0; k < numEl; ++k)
for (int i = 0; i < numMapNodes; ++i)
{
@@ -248,20 +262,13 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
nodesX(i,k) = v->x();
nodesY(i,k) = v->y();
nodesZ(i,k) = v->z();
- if (i < numMapNodes1) {
- nodesX1(i,k) = nodesX(i,k);
- nodesY1(i,k) = nodesY(i,k);
- nodesZ1(i,k) = nodesZ(i,k);
- }
}
- jfs->getSignedJacobian(nodesX, nodesY, nodesZ, jacobianB);
- jfs1->getSignedJacobian(nodesX1, nodesY1, nodesZ1, jac1B);
+ jfs->getScaledJacobian(nodesX, nodesY, nodesZ, jacobianB);
jfs->lag2Bez(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
- fullVector<double> jac1(numSamplingPt1);
#endif
for (int k = 0; k < numEl; ++k) {
@@ -269,25 +276,14 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
#ifdef _ANALYSECURVEDMESH_BLAS_
jacBez.setAsProxy(jacBezB, k);
jacobian.setAsProxy(jacobianB, k);
- jac1.setAsProxy(jac1B, k);
#else
- fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZ1(numMapNodes1,3);
+ fullMatrix<double> nodesXYZ(numMapNodes,3);
el[k]->getNodesCoord(nodesXYZ);
- nodesXYZ1.copy(nodesXYZ,0,numMapNodes1,0,3,0,0);
- jfs->getSignedJacobian(nodesXYZ,jacobian);
- jfs1->getSignedJacobian(nodesXYZ1,jac1);
+ jfs->getScaledJacobian(nodesXYZ,jacobian);
#endif
- // AmJ : avgJ is not the average Jac for quad, prism or hex
- double avgJ = sum(jac1) / jac1.size();
- if (avgJ < 0) {
- jacBez.scale(-1);
- jacobian.scale(-1);
- avgJ *= -1;
- }
-
int i;
- for (i = 0; i < numSamplingPt && jacobian(i) > _jacBreak * avgJ; ++i);
+ for (i = 0; i < numSamplingPt && jacobian(i) > _jacBreak; ++i);
if (i < numSamplingPt) {
invalids.push_back(el[k]);
++_numInvalid;
@@ -304,7 +300,7 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
jfs->lag2Bez(jacobian, jacBez);
#endif
- for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak * avgJ; ++i);
+ for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak; ++i);
if (i >= jacBez.size()) {
++_numValid;
continue;
@@ -386,7 +382,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
_numInvalid = 0;
_numValid = 0;
_numUncertain = 0;
- _min_Javg = .0, _max_Javg = .0, _avg_Javg = .0;
_min_pJmin = .0, _avg_pJmin = .0;
_min_ratioJ = .0, _avg_ratioJ = .0;
@@ -406,7 +401,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
break;
case 2 :
- Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
GFace *f = *it;
unsigned int numType[3] = {0, 0, 0};
@@ -421,7 +415,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
break;
case 1 :
- Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
GEdge *e = *it;
unsigned int numElement = e->getNumMeshElements();
@@ -435,7 +428,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
Msg::Error("I can't analyse any element.");
return;
}
- Msg::Info("Extrema of J_avg : %g, %g (avg: %g)", _min_Javg, _max_Javg, _avg_Javg/_numAnalysedEl);
Msg::Info("Minimum of min(~distortion) : %g", _min_pJmin);
Msg::Info("Average of min(~distortion) : %g", _avg_pJmin / _numAnalysedEl);
Msg::Info("Minimum of min(J) / max(J) : %g", _min_ratioJ);
@@ -447,24 +439,21 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
if (numEl < 1)
return;
- const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1);
- const JacobianBasis *jfs1 = el[0]->getJacobianFuncSpace(1);
- if (!jfs || !jfs1) {
+ const JacobianBasis *jfs = el[0]->getJacobianFuncSpace();
+ if (!jfs) {
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
- const int numSamplingPt = jfs->getNumJacNodes(), numSamplingPt1 = jfs1->getNumJacNodes();
- const int numMapNodes = jfs->getNumMapNodes(), numMapNodes1 = jfs1->getNumMapNodes();
+ const int numSamplingPt = jfs->getNumJacNodes();
+ const int numMapNodes = jfs->getNumMapNodes();
#ifdef _ANALYSECURVEDMESH_BLAS_
fullMatrix<double> jacobianB(numSamplingPt, numEl);
fullMatrix<double> jacBezB(numSamplingPt, numEl);
- fullMatrix<double> jac1B(numSamplingPt1, numEl);
fullVector<double> jacBez, jacobian, jac1;
fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
- fullMatrix<double> nodesX1(numMapNodes1, numEl), nodesY1(numMapNodes1, numEl), nodesZ1(numMapNodes1, numEl);
for (int k = 0; k < numEl; ++k)
for (int i = 0; i < numMapNodes; ++i)
{
@@ -472,25 +461,16 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
nodesX(i,k) = v->x();
nodesY(i,k) = v->y();
nodesZ(i,k) = v->z();
- if (i < numMapNodes1) {
- nodesX1(i,k) = nodesX(i,k);
- nodesY1(i,k) = nodesY(i,k);
- nodesZ1(i,k) = nodesZ(i,k);
- }
}
- jfs->getSignedJacobian(nodesX, nodesY, nodesZ, jacobianB);
- jfs1->getSignedJacobian(nodesX, nodesY, nodesZ, jac1B);
+ jfs->getScaledJacobian(nodesX, nodesY, nodesZ, jacobianB);
jfs->lag2Bez(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
- fullVector<double> jac1(numSamplingPt1);
#endif
fullVector<double> subJacBez(jfs->getNumSubNodes());
- _min_Javg = 1.7e308;
- _max_Javg = -1.7e308;
_min_pJmin = 1.7e308;
_min_ratioJ = 1.7e308;
@@ -502,24 +482,13 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
#ifdef _ANALYSECURVEDMESH_BLAS_
jacBez.setAsProxy(jacBezB, k);
jacobian.setAsProxy(jacobianB, k);
- jac1.setAsProxy(jac1B, k);
#else
- fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZ1(numMapNodes1,3);
+ fullMatrix<double> nodesXYZ(numMapNodes,3);
el[k]->getNodesCoord(nodesXYZ);
- nodesXYZ1.copy(nodesXYZ,0,numMapNodes1,0,3,0,0);
- jfs->getSignedJacobian(nodesXYZ,jacobian);
- jfs1->getSignedJacobian(nodesXYZ1,jac1);
+ jfs->getScaledJacobian(nodesXYZ,jacobian);
jfs->lag2Bez(jacobian, jacBez);
#endif
- // AmJ : avgJ is not the average Jac for quad, prism or hex
- double avgJ = sum(jac1) / jac1.size();
- if (avgJ < 0) {
- jacBez.scale(-1);
- jacobian.scale(-1);
- avgJ *= -1;
- }
-
double minJ, maxJ = minJ = jacobian(0);
for (int i = 1; i < numSamplingPt; ++i) {
@@ -527,18 +496,12 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
if (jacobian(i) > maxJ) maxJ = jacobian(i);
}
- double minB, maxB = minB = jacBez(0);//, avgJ = .0;
+ double minB, maxB = minB = jacBez(0);
for (int i = 1; i < numSamplingPt; ++i) {
if (jacBez(i) < minB) minB = jacBez(i);
if (jacBez(i) > maxB) maxB = jacBez(i);
- //avgJ += jacBez(i);
}
- //avgJ /= numSamplingPt;
-
- _avg_Javg += avgJ;
- _min_Javg = std::min(_min_Javg, avgJ);
- _max_Javg = std::max(_max_Javg, avgJ);
if (_maxDepth > 1 &&
(minJ - minB > _tol * (std::abs(minJ) + std::abs(minB)) / 2 ||
@@ -618,18 +581,18 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
delete bj;
}
}
- fwrite << minB/avgJ << " " << minB/maxB << "\r";
+ fwrite << minB << " " << minB/maxB << "\r";
if (data){
if (1-minB <= _tol * minJ && maxB-1 <= _tol * maxB)
(*data)[el[k]->getNum()].push_back(1.);
- else if (1-minB/avgJ <= 1e-8)
+ else if (1-minB <= 1e-8)
(*data)[el[k]->getNum()].push_back(1.);
else
- (*data)[el[k]->getNum()].push_back(minB/avgJ);
+ (*data)[el[k]->getNum()].push_back(minB);
}
- _min_pJmin = std::min(_min_pJmin, minB/avgJ);
- _avg_pJmin += minB/avgJ;
+ _min_pJmin = std::min(_min_pJmin, minB);
+ _avg_pJmin += minB;
_min_ratioJ = std::min(_min_ratioJ, minB/maxB);
_avg_ratioJ += minB/maxB;
}
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index a16c68c3c9687ad8e6e6938f77b8fca65009c9e5..b9153a68bde5d4e82cb3cfda94ebf3ce42da6e84 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -25,7 +25,6 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
int _numAnalysedEl;
int _numInvalid, _numValid, _numUncertain;
- double _min_Javg, _max_Javg, _avg_Javg;
double _min_pJmin, _avg_pJmin;
double _min_ratioJ, _avg_ratioJ;
@@ -53,32 +52,4 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
void hideValid_ShowInvalid(std::vector<MElement*> &invalids);
};
-class BezierJacobian;
-struct lessMinB {
- bool operator()(BezierJacobian*, BezierJacobian*) const;
-};
-struct lessMaxB {
- bool operator()(BezierJacobian*, BezierJacobian*) const;
-};
-
-class BezierJacobian
-{
-private:
- fullVector<double> _jacBez;
- double _minJ, _maxJ, _minB, _maxB; //Extremum of Jac at corners and of bezier values
- int _depthSub;
- const JacobianBasis *_jfs;
-
-public:
- BezierJacobian(fullVector<double> &, const JacobianBasis *, int depth);
- void subDivisions(fullVector<double> &vect) const
- {_jfs->subdivideBezierCoeff(_jacBez, vect);}
-
- inline int depth() const {return _depthSub;}
- inline double minJ() const {return _minJ;}
- inline double maxJ() const {return _maxJ;}
- inline double minB() const {return _minB;}
- inline double maxB() const {return _maxB;}
-};
-
#endif