Added: Hessian calculation for shells, stored in FiniteElement::H.

Added: Store polynomial degree in each direction in FiniteElement::p,q,r.

src/ASM/ASMs2D.C

@@ -1544,6 +1544,8 @@ bool ASMs2D::integrate (Integrand& integrand,
     for (size_t t = 0; t < threadGroups[g].size(); t++)
     {
       FiniteElement fe(p1*p2);
+      fe.p = p1 - 1;
+      fe.q = p2 - 1;
       Matrix   dNdu, Xnod, Jac;
       Matrix3D d2Ndu2, Hess;
       double   dXidu[2];
@@ -1707,8 +1709,12 @@ bool ASMs2D::integrate (Integrand& integrand,
 
             // Compute Hessian of coordinate mapping and 2nd order derivatives
             if (use2ndDer)
+            {
               if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,fe.dNdX))
                 ok = false;
+              else if (nsd > 2)
+                utl::Hessian(Hess,fe.H);
+            }
 
             // Compute G-matrix
             if (integrand.getIntegrandType() & Integrand::G_MATRIX)
@@ -1725,6 +1731,8 @@ bool ASMs2D::integrate (Integrand& integrand,
                 std::cout <<"d2NdX2 ="<< fe.d2NdX2;
               if (!fe.G.empty())
                 std::cout <<"G ="<< fe.G;
+              if (!fe.H.empty())
+                std::cout <<"H ="<< fe.H;
             }
 #endif
 
@@ -1816,6 +1824,8 @@ bool ASMs2D::integrate (Integrand& integrand,
     for (size_t t = 0; t < threadGroups[g].size(); t++)
     {
       FiniteElement fe(p1*p2);
+      fe.p = p1 - 1;
+      fe.q = p2 - 1;
       Matrix   dNdu, Xnod, Jac;
       Matrix3D d2Ndu2, Hess;
       double   dXidu[2];
@@ -1901,8 +1911,12 @@ bool ASMs2D::integrate (Integrand& integrand,
 
           // Compute Hessian of coordinate mapping and 2nd order derivatives
           if (use2ndDer)
+          {
             if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,fe.dNdX))
               ok = false;
+            else if (nsd > 2)
+              utl::Hessian(Hess,fe.H);
+          }
 
           // Compute G-matrix
           if (integrand.getIntegrandType() & Integrand::G_MATRIX)
@@ -2175,6 +2189,8 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
   size_t firstp = iit == firstBp.end() ? 0 : iit->second;
 
   FiniteElement fe(p1*p2);
+  fe.p = p1 - 1;
+  fe.q = p2 - 1;
   fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
   fe.u = gpar[0](1,1);
   fe.v = gpar[1](1,1);
@@ -2617,9 +2633,11 @@ bool ASMs2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
   this->getNodalCoordinates(Xnod);
 
   FiniteElement fe(p1*p2,firstIp);
-  Vector        solPt;
+  fe.p = p1 - 1;
-  Matrix        dNdu, Jac;
+  fe.q = p2 - 1;
-  Matrix3D      d2Ndu2, Hess;
+  Vector   solPt;
+  Matrix   dNdu, Jac;
+  Matrix3D d2Ndu2, Hess;
 
   // Evaluate the secondary solution field at each point
   for (size_t i = 0; i < nPoints; i++, fe.iGP++)
@@ -2654,8 +2672,12 @@ bool ASMs2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
 
     // Compute Hessian of coordinate mapping and 2nd order derivatives
     if (use2ndDer)
+    {
       if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xtmp,d2Ndu2,fe.dNdX))
         continue;
+      else if (nsd > 2)
+        utl::Hessian(Hess,fe.H);
+    }
 
     // Store tangent vectors in fe.G for shells
     if (nsd > 2) fe.G = Jac;

src/ASM/FiniteElement.C

@@ -16,7 +16,8 @@
 
 std::ostream& FiniteElement::write (std::ostream& os) const
 {
-  os <<"FiniteElement: iel="<< iel <<" iGP="<< iGP <<" p="<< p
+  os <<"FiniteElement: iel="<< iel <<" iGP="<< iGP
+     <<"\n               p, q, r: "<< p <<" "<< q <<" "<< r
      <<"\n               u, v, w: "<< u <<" "<< v <<" "<< w
      <<"\n               xi, eta, zeta: "<< xi <<" "<< eta <<" "<< zeta
      <<"\n               h, detJxW: "<< h <<" "<< detJxW << std::endl;
@@ -24,6 +25,7 @@ std::ostream& FiniteElement::write (std::ostream& os) const
   if (!dNdX.empty())   os <<"dNdX:"<< dNdX;
   if (!d2NdX2.empty()) os <<"d2NdX2:"<< d2NdX2;
   if (!G.empty())      os <<"G:"<< G;
+  if (!H.empty())      os <<"H:"<< H;
   if (!Navg.empty())   os <<"Navg:"<< Navg;
   if (!Xn.empty())     os <<"Xn:"<< Xn;
   if (!Te.isZero(0.0)) os <<"Te:\n"<< Te;

src/ASM/FiniteElement.h

@@ -28,19 +28,19 @@ class FiniteElement
 public:
   //! \brief Default constructor.
   explicit FiniteElement(size_t n = 0, size_t i = 0) : iGP(i), N(n), Te(3)
-  { iel = p = 0; u = v = w = xi = eta = zeta = h = 0.0; detJxW = 1.0; }
+  { iel = p = q = r = 0; u = v = w = xi = eta = zeta = h = 0.0; detJxW = 1.0; }
 
   //! \brief Empty destructor.
   virtual ~FiniteElement() {}
 
+  //! \brief Returns the number of bases.
+  virtual size_t getNoBasis() const { return 1; }
+
   //! \brief Returns a const reference to the basis function values.
   virtual const Vector& basis(char) const { return N; }
   //! \brief Returns a reference to the basis function values.
   virtual Vector& basis(char) { return N; }
 
-  //! \brief Returns number of bases
-  virtual size_t getNoBasis() const { return 1; }
-
   //! \brief Returns a const reference to the basis function derivatives.
   virtual const Matrix& grad(char) const { return dNdX; }
   //! \brief Returns a reference to the basis function derivatives.
@@ -74,11 +74,14 @@ public:
   Vector     N;    //!< Basis function values
   Matrix    dNdX;  //!< First derivatives (gradient) of the basis functions
   Matrix3D d2NdX2; //!< Second derivatives of the basis functions
-  Matrix     G;    //!< Matrix used for stabilized methods
+  Matrix     G;    //!< Covariant basis / Matrix used for stabilized methods
+  Matrix     H;    //!< Hessian
 
   // Element quantities
   int                 iel;  //!< Element identifier
-  short int           p;    //!< Polynomial order of the basis functions
+  short int           p;    //!< Polynomial order of the basis in u-direction
+  short int           q;    //!< Polynomial order of the basis in v-direction
+  short int           r;    //!< Polynomial order of the basis in r-direction
   double              h;    //!< Characteristic element size/diameter
   Vec3Vec             XC;   //!< Array with element corner coordinate vectors
   Vector              Navg; //!< Volume-averaged basis function values
@@ -101,14 +104,14 @@ public:
   //! \brief Empty destructor.
   virtual ~MxFiniteElement() {}
 
+  //! \brief Returns the number of bases.
+  virtual size_t getNoBasis() const { return 1+Nx.size(); }
+
   //! \brief Returns a const reference to the basis function values.
   virtual const Vector& basis(char b) const { return b == 1 ? N : Nx[b-2]; }
   //! \brief Returns a reference to the basis function values.
   virtual Vector& basis(char b) { return b == 1 ? N : Nx[b-2]; }
 
-  //! \brief Returns number of bases
-  virtual size_t getNoBasis() const { return 1+Nx.size(); }
-
   //! \brief Returns a const reference to the basis function derivatives.
   virtual const Matrix& grad(char b) const { return b == 1 ? dNdX : dNxdX[b-2]; }
   //! \brief Returns a reference to the basis function derivatives.

src/ASM/LR/ASMu2D.C

@@ -1043,6 +1043,8 @@ bool ASMu2D::integrate (Integrand& integrand,
 
       FiniteElement fe;
       fe.iel = MLGE[iel-1];
+      fe.p   = lrspline->order(0) - 1;
+      fe.q   = lrspline->order(1) - 1;
       Matrix   dNdu, Xnod, Jac;
       Matrix3D d2Ndu2, Hess;
       double   dXidu[2];
@@ -1190,11 +1192,15 @@ bool ASMu2D::integrate (Integrand& integrand,
 
           // Compute Hessian of coordinate mapping and 2nd order derivatives
           if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
+          {
             if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
             {
               ok = false;
               continue;
             }
+            else if (nsd > 2)
+              utl::Hessian(Hess,fe.H);
+          }
 
           // Compute G-matrix
           if (integrand.getIntegrandType() & Integrand::G_MATRIX)
@@ -1281,6 +1287,8 @@ bool ASMu2D::integrate (Integrand& integrand,
 
       FiniteElement fe;
       fe.iel = MLGE[iel-1];
+      fe.p   = lrspline->order(0) - 1;
+      fe.q   = lrspline->order(1) - 1;
       Matrix   dNdu, Xnod, Jac;
       Matrix3D d2Ndu2, Hess;
       Vec4     X;
@@ -1348,11 +1356,15 @@ bool ASMu2D::integrate (Integrand& integrand,
 
         // Compute Hessian of coordinate mapping and 2nd order derivatives
         if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
+        {
           if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
           {
             ok = false;
             continue;
           }
+          else if (nsd > 2)
+            utl::Hessian(Hess,fe.H);
+        }
 
         // Store tangent vectors in fe.G for shells
         if (nsd > 2) fe.G = Jac;
@@ -1473,6 +1485,8 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
 
     // Initialize element quantities
     fe.iel = MLGE[iel-1];
+    fe.p   = lrspline->order(0) - 1;
+    fe.q   = lrspline->order(1) - 1;
     fe.xi  = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
     LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel-1].size(),
                                                   fe.iel,true);
@@ -1742,6 +1756,8 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
     return false;
 
   FiniteElement fe;
+  fe.p = lrspline->order(0) - 1;
+  fe.q = lrspline->order(1) - 1;
   Vector   ptSol;
   Matrix   dNdu, dNdX, Jac, Xnod, eSol, ptDer;
   Matrix3D d2Ndu2, d2NdX2, Hess, ptDer2;
@@ -1882,6 +1898,8 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
     return false;
 
   FiniteElement fe(0,firstIp);
+  fe.p = lrspline->order(0) - 1;
+  fe.q = lrspline->order(1) - 1;
   Vector   solPt;
   Matrix   dNdu, Jac, Xnod;
   Matrix3D d2Ndu2, Hess;
@@ -1915,8 +1933,12 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
 
     // Compute Hessian of coordinate mapping and 2nd order derivatives
     if (use2ndDer)
+    {
       if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
         continue;
+      else if (nsd > 2)
+        utl::Hessian(Hess,fe.H);
+    }
 
     // Store tangent vectors in fe.G for shells
     if (nsd > 2) fe.G = Jac;

src/Utility/CoordinateMapping.C

@@ -41,14 +41,7 @@ Real utl::Jacobian (matrix<Real>& J, matrix<Real>& dNdX,
     // Special for two-parametric elements in 3-dimensional space (shell)
     dNdX = dNdu;
     // Compute the length of the normal vector {X},u x {X},v
-    Vec3 a1 = J.getColumn(1);
+    detJ = Vec3(J.getColumn(1),J.getColumn(2)).length();
-    Vec3 a2 = J.getColumn(2);
-    detJ = Vec3(a1,a2).length();
-    // Store the two tangent vectors in J
-    a1.normalize();
-    a2.normalize();
-    J.fillColumn(1,a1.ptr());
-    J.fillColumn(2,a2.ptr());
   }
   else
   {
@@ -115,9 +108,6 @@ Real utl::Jacobian (matrix<Real>& J, Vec3& n, matrix<Real>& dNdX,
     {
       // Special treatment for two-parametric elements in 3D space (shells)
       dNdX = dNdu;
-      v1.normalize();
-      J.fillColumn(t1,v1.ptr());
-      J.fillColumn(t2,v2.ptr());
       return dS;
     }
   }
@@ -161,6 +151,7 @@ bool utl::Hessian (matrix3d<Real>& H, matrix3d<Real>& d2NdX2,
   else if (Ji.cols() <= 2 && Ji.rows() > Ji.cols())
   {
     // Special treatment for one-parametric elements in multi-dimension space
+    // as well as two-parametric elements in 3D space (shells)
     d2NdX2 = d2Ndu2;
     return true;
   }
@@ -199,6 +190,24 @@ bool utl::Hessian (matrix3d<Real>& H, matrix3d<Real>& d2NdX2,
 }
 
 
+void utl::Hessian (const matrix3d<Real>& Hess, matrix<Real>& H)
+{
+  size_t i, n = Hess.dim(2);
+  if (n == Hess.dim(3))
+  {
+    H.resize(Hess.dim(1),n*(n+1)/2);
+    for (i = 1; i <= n; i++)
+      H.fillColumn(i,Hess.getColumn(i,i));
+    for (i = 2; i <= n; i++)
+      H.fillColumn(n+i-1,Hess.getColumn(1,i));
+    if (n > 2)
+      H.fillColumn(n+n,Hess.getColumn(2,3));
+  }
+  else
+    H.clear();
+}
+
+
 void utl::getGmat (const matrix<Real>& Ji, const Real* du, matrix<Real>& G)
 {
   size_t nsd = Ji.cols();

src/Utility/CoordinateMapping.h

@@ -64,13 +64,17 @@ namespace utl
   //! \param[in] X Matrix of element nodal coordinates
   //! \param[in] d2Ndu2 Second order derivatives of basis functions
   //! \param[in] dNdX First order derivatives of basis functions
-  //! \param[in] geoMapping If true, calculate geometry mapping
+  //! \param[in] geoMapping If \e true, calculate geometry mapping
   //! \return \e false if matrix dimensions are incompatible, otherwise \e true
-  //! \details If geoMapping is true, H is output, else H is input and assumed to be already calculated.
+  //!
+  //! \details If geoMapping is \e true, H is output,
+  //! else H is input and assumed to be already calculated in a previous call.
   bool Hessian(matrix3d<Real>& H, matrix3d<Real>& d2NdX2,
                const matrix<Real>& Ji, const matrix<Real>& X,
                const matrix3d<Real>& d2Ndu2, const matrix<Real>& dNdX,
-               bool geoMapping=true);
+               bool geoMapping = true);
+  //! \brief Convert a Hessian from a matrix3d to a matrix assuming symmetry.
+  void Hessian(const matrix3d<Real>& Hess, matrix<Real>& H);
 
   //! \brief Compute the stabilization matrix \b G from the Jacobian inverse.
   //! \param[in] Ji The inverse of the Jacobian matrix

src/Utility/Test/TestCoordinateMapping.C

@@ -221,35 +221,58 @@ TEST(TestCoordinateMapping, JacobianShell)
   SIMShell sim;
   ASSERT_TRUE(sim.createDefaultModel());
   ASMs2D& p = static_cast<ASMs2D&>(*sim.getPatch(1));
+  ASSERT_TRUE(p.raiseOrder(1,1));
   ASSERT_TRUE(p.uniformRefine(0,3));
   ASSERT_TRUE(p.uniformRefine(1,3));
   ASSERT_TRUE(sim.createFEMmodel());
 
   Vector N;
   Matrix X, dNdU;
-  p.extractBasis(0.25, 0.25, N, dNdU);
+  Matrix3D d2Ndu2;
+  p.extractBasis(0.25, 0.25, N, dNdU, d2Ndu2);
   p.getElementCoordinates(X, 2);
 
-  Matrix J, dNdX;
+  Matrix J, dNdX, Hs;
+  Matrix3D H, d2NdX2;
   Real det = utl::Jacobian(J, dNdX, X, dNdU, true);
+  utl::Hessian(H, d2NdX2, J, X, d2Ndu2, dNdX);
+  utl::Hessian(H, Hs);
 
-  const double dndx[] = {
+  const double j[] = {
-    -4.0, 4.0, 0.0, 0.0,
+    1.0, 0.0, 0.0,
-    -4.0, 0.0, 4.0, 0.0
+    0.0, 0.5, 0.0
   };
 
-  EXPECT_FLOAT_EQ(det, 1.0);
+  const double h[] = {
-  EXPECT_EQ(J.rows(), 3U);
+    0.0, 0.0, 0.0,
-  EXPECT_EQ(J.cols(), 2U);
+    0.0, 2.0, 0.0,
-  for (size_t i = 1; i <= J.rows(); i++)
+    0.0, 0.0, 0.0
-    for (size_t j = 1; j <= J.cols(); j++)
+  };
-      EXPECT_FLOAT_EQ(J(i,j), i == j ? 1.0 : 0.0);
+
+  const double dndx[] = {
+    -2.0, 2.0, 0.0,-2.0, 2.0, 0.0, 0.0, 0.0, 0.0,
+    -2.0,-2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0
+  };
+
+  EXPECT_FLOAT_EQ(det, 0.5);
 
   size_t i = 0;
-  EXPECT_EQ(dNdX.rows(), 4U);
+  EXPECT_EQ(J.rows(), 3U);
+  EXPECT_EQ(J.cols(), 2U);
+  for (double v : J)
+    EXPECT_FLOAT_EQ(j[i++], v);
+
+  i = 0;
+  EXPECT_EQ(dNdX.rows(), 9U);
   EXPECT_EQ(dNdX.cols(), 2U);
   for (double v : dNdX)
     EXPECT_FLOAT_EQ(dndx[i++], v);
+
+  i = 0;
+  EXPECT_EQ(Hs.rows(), 3U);
+  EXPECT_EQ(Hs.cols(), 3U);
+  for (double v : Hs)
+    EXPECT_FLOAT_EQ(h[i++], v);
 }