diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 9f03c226d9626acb95aaa25f1e3399b46f900515..9a5f07bfe76cd206f99c747942a0ae075065c82a 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -112,19 +112,12 @@ void MElement::scaledJacRange(double &jmin, double &jmax)
{
jmin = jmax = 1.0;
#if defined(HAVE_MESH)
- if (getPolynomialOrder() == 1) return;
- const JacobianBasis *jac = getJacobianFuncSpace(),*jac1 = getJacobianFuncSpace(1); // Jac. and prim. Jac. basis
- const int numJacNodes = jac->getNumJacNodes(), numJac1Nodes = jac1->getNumJacNodes();
- const int numMapNodes = jac->getNumMapNodes(), numMap1Nodes = jac1->getNumMapNodes();
- fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZ1(numMap1Nodes,3);
+ const JacobianBasis *jac = getJacobianFuncSpace();
+ const int numJacNodes = jac->getNumJacNodes();
+ fullMatrix<double> nodesXYZ(jac->getNumMapNodes(),3);
getNodesCoord(nodesXYZ);
- nodesXYZ1.copy(nodesXYZ,0,numMap1Nodes,0,3,0,0);
- fullVector<double> Ji(numJacNodes), J1i(numJac1Nodes), J1Ji(numJacNodes);
- jac->getSignedJacobian(nodesXYZ,Ji);
- jac1->getSignedJacobian(nodesXYZ1,J1i);
- jac->primJac2Jac(J1i,J1Ji);
- fullVector<double> SJi(numJacNodes), Bi(numJacNodes); // Calc scaled Jacobian -> Bezier
- for (int i=0; i<numJacNodes; i++) SJi(i) = Ji(i)/J1Ji(i);
+ fullVector<double> SJi(numJacNodes), Bi(numJacNodes);
+ jac->getScaledJacobian(nodesXYZ,SJi);
jac->lag2Bez(SJi,Bi);
jmin = *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
jmax = *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 4c8f420a9ca1f9fa9b0a4d89bdaeda7905a8f6e5..29ea8fb90224a7b1d986974d87724eace6d35b6d 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -1179,8 +1179,8 @@ const JacobianBasis *JacobianBasis::find(int tag)
-// Computes (unit) normals to straight line element
-void JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
+// Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value)
+double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
{
fullVector<double> dxyzdXbar(3);
@@ -1206,11 +1206,13 @@ void JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMat
const double norm1 = sqrt(result(1,0)*result(1,0)+result(1,1)*result(1,1)+result(1,2)*result(1,2));
result(1,0) /= norm1; result(1,1) /= norm1; result(1,2) /= norm1;
+ return sqrt(dxyzdXbar(0)*dxyzdXbar(0)+dxyzdXbar(1)*dxyzdXbar(1)+dxyzdXbar(2)*dxyzdXbar(2));
+
}
-// Computes (unit) normal to straight surface element (with norm as return value)
+// Computes (unit) normal to straight surface element at barycenter (with norm as return value)
double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
{
@@ -1236,9 +1238,9 @@ double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMa
-inline double calcDet3D(double dxdX, double dydX, double dzdX,
- double dxdY, double dydY, double dzdY,
- double dxdZ, double dydZ, double dzdZ)
+static inline double calcDet3D(double dxdX, double dydX, double dzdX,
+ double dxdY, double dydY, double dzdY,
+ double dxdZ, double dydZ, double dzdZ)
{
return dxdX*dydY*dzdZ + dxdY*dydZ*dzdX + dydX*dzdY*dxdZ
- dxdZ*dydY*dzdX - dxdY*dydX*dzdZ - dydZ*dzdY*dxdX;
@@ -1246,6 +1248,29 @@ inline double calcDet3D(double dxdX, double dydX, double dzdX,
+// Returns absolute value of Jacobian of straight volume element at barycenter
+double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
+{
+
+ double dxdX = 0, dydX = 0, dzdX = 0, dxdY = 0, dydY = 0, dzdY = 0, dxdZ = 0, dydZ = 0, dzdZ = 0;
+ for (int j=0; j<numPrimMapNodes; j++) {
+ dxdX += primGradShapeBarX(j)*nodesXYZ(j,0);
+ dydX += primGradShapeBarX(j)*nodesXYZ(j,1);
+ dzdX += primGradShapeBarX(j)*nodesXYZ(j,2);
+ dxdY += primGradShapeBarY(j)*nodesXYZ(j,0);
+ dydY += primGradShapeBarY(j)*nodesXYZ(j,1);
+ dzdY += primGradShapeBarY(j)*nodesXYZ(j,2);
+ dxdZ += primGradShapeBarZ(j)*nodesXYZ(j,0);
+ dydZ += primGradShapeBarZ(j)*nodesXYZ(j,1);
+ dzdZ += primGradShapeBarZ(j)*nodesXYZ(j,2);
+ }
+
+ return fabs(calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ));
+
+}
+
+
+
// Calculate (signed) Jacobian at mapping's nodes for one element, with normal vectors to
// straight element for regularization
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
@@ -1280,6 +1305,44 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVe
+// Calculate scaled (signed) Jacobian at mapping's nodes for one element, with normal vectors to
+// straight element for regularization and scaling
+void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
+{
+
+ switch (bezier->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
+ getSignedJacobian(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
+ getSignedJacobian(nodesXYZ,normal,jacobian);
+ break;
+ }
+
+ case 0 :
+ case 3 : {
+ fullMatrix<double> dum;
+ const double scale = 1./getPrimJac3D(nodesXYZ);
+ getSignedJacobian(nodesXYZ,dum,jacobian);
+ jacobian.scale(scale);
+ break;
+ }
+
+ }
+
+}
+
+
+
// Calculate (signed) Jacobian at mapping's nodes for one element, given vectors for regularization
// of line and surface Jacobians in 3D
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index c8da26ec32b2b7e5e3ccc41eb13fda8ef237e415..8901c1d9a5a32f265916ec94abb4af1a70d68e43 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -46,8 +46,9 @@ class JacobianBasis {
inline int getNumDivisions() const { return bezier->numDivisions; }
inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
inline int getNumLagPts() const { return bezier->numLagPts; }
- void getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
+ double getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
double getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
+ double getPrimJac3D(const fullMatrix<double> &nodesXYZ) const;
void getSignedJacobian(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals, fullVector<double> &jacobian) const;
void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
@@ -58,6 +59,7 @@ class JacobianBasis {
void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
+ void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
bezier->matrixLag2Bez.mult(jac,bez);
}