Some cleanup of mixed LR

src/ASM/ASMs2Dmx.C

@@ -215,8 +215,11 @@ bool ASMs2Dmx::generateFEMTopology ()
     else if (ASMmxBase::Type == ASMmxBase::SUBGRID) {
       projB = proj = m_basis.front()->clone();
       altProjBasis = ASMmxBase::raiseBasis(surf);
-    } else
-      projB = proj = m_basis[2-ASMmxBase::geoBasis]->clone();
+    }
+    else if (geoBasis < 3)
+      projB = proj = m_basis[2-geoBasis]->clone();
+    else
+      return false; // Logic error
   }
   delete surf;
   geomB = surf = m_basis[geoBasis-1]->clone();
@@ -657,25 +660,12 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
 
             // Compute Jacobian inverse of the coordinate mapping and
             // basis function derivatives w.r.t. Cartesian coordinates
-            fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,
-                                      dNxdu[geoBasis-1]);
-            if (fe.detJxW == 0.0) continue; // skip singular points
-
-            for (size_t b = 0; b < m_basis.size(); ++b)
-              if (b != (size_t)geoBasis-1)
-                fe.grad(b+1).multiply(dNxdu[b],Jac);
+            if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis))
+              continue; // skip singular points
 
             // Compute Hessian of coordinate mapping and 2nd order derivatives
-            if (use2ndDer) {
-              if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod,
-                                d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true))
-                ok = false;
-              for (size_t b = 0; b < m_basis.size() && ok; ++b)
-                if ((int)b != geoBasis)
-                  if (!utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod,
-                                    d2Nxdu2[b],fe.grad(b+1),false))
-                    ok = false;
-            }
+            if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis))
+              ok = false;
 
             // Compute G-matrix
             if (integrand.getIntegrandType() & Integrand::G_MATRIX)
@@ -1197,12 +1187,8 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
 
     // Compute Jacobian inverse of the coordinate mapping and
     // basis function derivatives w.r.t. Cartesian coordinates
-    fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xtmp,dNxdu[geoBasis-1]);
-    if (fe.detJxW == 0.0) continue; // skip singular points
-
-    for (size_t b = 1; b <= m_basis.size(); b++)
-      if (b != (size_t)geoBasis)
-        fe.grad(b).multiply(dNxdu[b-1],Jac);
+    if (!fe.Jacobian(Jac,Xtmp,dNxdu,geoBasis))
+      continue; // skip singular points
 
     // Cartesian coordinates of current integration point
     utl::Point X4(Xtmp*fe.basis(geoBasis),{fe.u,fe.v});

src/ASM/ASMs2DmxLag.C

@@ -311,12 +311,8 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
                 ok = false;
 
             // Compute Jacobian inverse of coordinate mapping and derivatives
-            fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,
-                                      dNxdu[geoBasis-1]);
-            if (fe.detJxW == 0.0) continue; // skip singular points
-            for (size_t b = 0; b < nxx.size(); ++b)
-              if (b != (size_t)geoBasis-1)
-                fe.grad(b+1).multiply(dNxdu[b],Jac);
+            if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis))
+              continue; // skip singular points
 
             // Cartesian coordinates of current integration point
             X.assign(Xnod * fe.basis(geoBasis));
@@ -423,10 +419,11 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
                                       elem_sizes[b][1],xi[1]))
             ok = false;
 
-	// Compute basis function derivatives and the edge normal
-	fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
+        // Compute basis function derivatives and the edge normal
+        fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
                                   dNxdu[geoBasis-1],t1,t2);
-	if (fe.detJxW == 0.0) continue; // skip singular points
+        if (fe.detJxW == 0.0) continue; // skip singular points
+
         for (size_t b = 0; b < nxx.size(); ++b)
           if (b != (size_t)geoBasis-1)
             fe.grad(b+1).multiply(dNxdu[b],Jac);
@@ -526,13 +523,10 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
                                       elem_sizes[b][1],eta))
 	  return false;
 
-	// Compute the Jacobian inverse
-        fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]);
-        if (fe.detJxW == 0.0) continue; // skip singular points
-
-        for (size_t b = 1; b <= nxx.size(); b++)
-          if (b != (size_t)geoBasis)
-            fe.grad(b).multiply(dNxdu[b-1],Jac);
+        // Compute Jacobian inverse of the coordinate mapping and
+        // basis function derivatives w.r.t. Cartesian coordinates
+        if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis))
+          continue; // skip singular points
 
 	// Now evaluate the solution field
 	if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),

src/ASM/ASMs3Dmx.C

@@ -216,8 +216,11 @@ bool ASMs3Dmx::generateFEMTopology ()
     else if (ASMmxBase::Type == ASMmxBase::SUBGRID) {
       projB = proj = m_basis.front()->clone();
       altProjBasis = ASMmxBase::raiseBasis(svol);
-    } else
-      projB = proj = m_basis[2-ASMmxBase::geoBasis]->clone();
+    }
+    else if (geoBasis < 3)
+      projB = proj = m_basis[2-geoBasis]->clone();
+    else
+      return false; // Logic error
   }
   delete svol;
   geomB = svol = m_basis[geoBasis-1]->clone();
@@ -620,25 +623,12 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
 
               // Compute Jacobian inverse of the coordinate mapping and
               // basis function derivatives w.r.t. Cartesian coordinates
-              fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,
-                                        dNxdu[geoBasis-1]);
-              if (fe.detJxW == 0.0) continue; // skip singular points
-
-              for (size_t b = 0; b < m_basis.size(); ++b)
-                if (b != (size_t)geoBasis-1)
-                  fe.grad(b+1).multiply(dNxdu[b],Jac);
+              if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis))
+                continue; // skip singular points
 
               // Compute Hessian of coordinate mapping and 2nd order derivatives
-              if (use2ndDer) {
-                if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod,
-                                  d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true))
-                  ok = false;
-                for (size_t b = 0; b < m_basis.size() && ok; ++b)
-                  if ((int)b != geoBasis)
-                    if (!utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod,
-                                      d2Nxdu2[b],fe.grad(b+1),false))
-                      ok = false;
-              }
+              if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis))
+                ok = false;
 
               // Compute G-matrix
               if (integrand.getIntegrandType() & Integrand::G_MATRIX)
@@ -1246,13 +1236,9 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
       SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),dNxdu[b]);
 
     // Compute Jacobian inverse of the coordinate mapping and
-    // basis function derivatives w.r.t. Cartesian coordinates
-    fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xtmp,dNxdu[geoBasis-1]);
-    if (fe.detJxW == 0.0) continue; // skip singular points
-
-    for (size_t b = 1; b <= m_basis.size(); b++)
-      if (b != (size_t)geoBasis)
-        fe.grad(b).multiply(dNxdu[b-1],Jac);
+    // basis function derivatives w.r.t. Cartesian coordinate
+    if (!fe.Jacobian(Jac,Xtmp,dNxdu,geoBasis))
+      continue; // skip singular points
 
     // Cartesian coordinates of current integration point
     utl::Point X4(Xtmp * fe.basis(geoBasis),{fe.u,fe.v,fe.w});

src/ASM/ASMs3DmxLag.C

@@ -348,12 +348,8 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
                   ok = false;
 
               // Compute Jacobian inverse of coordinate mapping and derivatives
-              fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,
-                                        dNxdu[geoBasis-1]);
-              if (fe.detJxW == 0.0) continue; // skip singular points
-              for (size_t b = 0; b < nxx.size(); ++b)
-                if (b != (size_t)geoBasis-1)
-                  fe.grad(b+1).multiply(dNxdu[b],Jac);
+              if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis))
+                continue; // skip singular points
 
               // Cartesian coordinates of current integration point
               X.assign(Xnod * fe.basis(geoBasis));
@@ -494,6 +490,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
 	    fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod,
                                       dNxdu[geoBasis-1],t1,t2);
 	    if (fe.detJxW == 0.0) continue; // skip singular points
+
             for (size_t b = 0; b < nxx.size(); ++b)
               if (b != (size_t)geoBasis-1)
                 fe.grad(b+1).multiply(dNxdu[b],Jac);
@@ -597,14 +594,10 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
                                         elem_sizes[b][2],fe.zeta))
               return false;
 
-	  // Compute the Jacobian inverse
-          fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,
-                                    dNxdu[geoBasis-1]);
-          if (fe.detJxW == 0.0) continue; // skip singular points
-
-          for (size_t b = 1; b <= nxx.size(); b++)
-            if (b != (size_t)geoBasis)
-              fe.grad(b).multiply(dNxdu[b-1],Jac);
+          // Compute Jacobian inverse of the coordinate mapping and
+          // basis function derivatives w.r.t. Cartesian coordinates
+          if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis))
+            continue; // skip singular points
 
 	  // Now evaluate the solution field
 	  if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis),

src/ASM/FiniteElement.C

@@ -12,6 +12,7 @@
 //==============================================================================
 
 #include "FiniteElement.h"
+#include "CoordinateMapping.h"
 #include "Vec3Oper.h"
 
 
@@ -65,3 +66,58 @@ std::ostream& MxFiniteElement::write (std::ostream& os) const
   }
   return os;
 }
+
+
+/*!
+  This method also calculates the first-derivatives of the basis functions
+  with respect to the Cartesian coordinates, using the same geometry mapping
+  for all bases.
+*/
+
+bool MxFiniteElement::Jacobian (Matrix& Jac, const Matrix& Xnod,
+                                const std::vector<Matrix>& dNxdu,
+                                unsigned short int gBasis)
+{
+  detJxW = utl::Jacobian(Jac,
+                         gBasis > 1 ? dMdX[gBasis-2] : dNdX,
+                         Xnod, dNxdu[gBasis-1]);
+  if (detJxW == 0.0) return false; // singular point
+
+  if (gBasis > 1)
+    dNdX.multiply(dNxdu.front(),Jac);
+
+  for (size_t basis = 2; basis <= dNxdu.size(); basis++)
+    if (basis != gBasis)
+      dMdX[basis-2].multiply(dNxdu[basis-1],Jac);
+
+  return true;
+}
+
+
+/*!
+  This method also calculates the second-derivatives of the basis functions
+  with respect to the Cartesian coordinates, using the same geometry mapping
+  for all bases.
+*/
+
+bool MxFiniteElement::Hessian (Matrix3D& Hess, const Matrix& Jac,
+                               const Matrix& Xnod,
+                               const std::vector<Matrix3D>& d2Nxdu2,
+                               unsigned short int gBasis)
+{
+  bool ok = utl::Hessian(Hess,
+                         gBasis > 1 ? d2MdX2[gBasis-2] : d2NdX2,
+                         Jac, Xnod, d2Nxdu2[gBasis-1],
+                         gBasis > 1 ? dMdX[gBasis-2] : dNdX);
+
+  if (ok && gBasis > 1)
+    ok = utl::Hessian(Hess, d2NdX2, Jac, Xnod,
+                      d2Nxdu2.front(), dNdX, false);
+
+  for (size_t basis = 2; basis <= d2Nxdu2.size() && ok; basis++)
+    if (basis != gBasis)
+      ok = utl::Hessian(Hess, d2MdX2[basis-2], Jac, Xnod,
+                        d2Nxdu2[basis-1], dMdX[basis-2], false);
+
+  return ok;
+}

src/ASM/FiniteElement.h

@@ -49,13 +49,8 @@ class FiniteElement : public ItgPoint
 
   //! \brief Returns a const reference to the basis function 2nd-derivatives.
   virtual const Matrix3D& hess(char) const { return d2NdX2; }
-  //! \brief Returns a reference to the basis function 2nd-derivatives.
-  virtual Matrix3D& hess(char) { return d2NdX2; }
-
   //! \brief Returns a const reference to the basis function 3nd-derivatives.
   virtual const Matrix4D& hess2(char) const { return d3NdX3; }
-  //! \brief Returns a reference to the basis function 3nd-derivatives.
-  virtual Matrix4D& hess2(char) { return d3NdX3; }
 
 protected:
   //! \brief Writes the finite element object to the given output stream.
@@ -118,13 +113,25 @@ class MxFiniteElement : public FiniteElement
 
   //! \brief Returns a const reference to the basis function 2nd-derivatives.
   virtual const Matrix3D& hess(char b) const { return b < 2 ? d2NdX2 : d2MdX2[b-2]; }
-  //! \brief Returns a reference to the basis function 2nd-derivatives.
-  virtual Matrix3D& hess(char b) { return b < 2 ? d2NdX2 : d2MdX2[b-2]; }
-
   //! \brief Returns a const reference to the basis function 3rd-derivatives.
   virtual const Matrix4D& hess2(char b) const { return b < 2 ? d3NdX3 : d3MdX3[b-2]; }
-  //! \brief Returns a reference to the basis function 3rd-derivatives.
-  virtual Matrix4D& hess2(char b) { return b < 2 ? d3NdX3 : d3MdX3[b-2]; }
+
+  //! \brief Sets up the Jacobian matrix of the coordinate mapping.
+  //! \param[out] Jac The inverse of the Jacobian matrix
+  //! \param[in] Xnod Matrix of element nodal coordinates
+  //! \param[in] dNxdu First order derivatives of basis functions
+  //! \param[in] gBasis 1-based index of basis representing the geometry
+  bool Jacobian(Matrix& Jac, const Matrix& Xnod,
+                const std::vector<Matrix>& dNxdu, unsigned short int gBasis);
+
+  //! \brief Sets up the Hessian matrix of the coordinate mapping.
+  //! \param[out] Hess The Hessian matrix
+  //! \param[in] Jac The inverse of the Jacobian matrix
+  //! \param[in] Xnod Matrix of element nodal coordinates
+  //! \param[in] d2Nxdu2 Second order derivatives of basis functions
+  //! \param[in] gBasis 1-based index of basis representing the geometry
+  bool Hessian(Matrix3D& Hess, const Matrix& Jac, const Matrix& Xnod,
+               const std::vector<Matrix3D>& d2Nxdu2, unsigned short int gBasis);
 
 protected:
   //! \brief Writes the finite element object to the given output stream.