Changed definition of scaled Jacobian (affects disto measure)

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());

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,7 +1238,7 @@ double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMa
 
 
 
-inline double calcDet3D(double dxdX, double dydX, double dzdX,
+static inline double calcDet3D(double dxdX, double dydX, double dzdX,
                                double dxdY, double dydY, double dzdY,
                                double dxdZ, double dydZ, double dzdZ)
 {
@@ -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,

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);
   }