From 89c065edc4b04d65f2c605a16ca010d8d2743e4d Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Thu, 31 Aug 2023 09:34:19 +0200
Subject: [PATCH] MxFiniteElement::Hessian

support for a separate geometry basis. checks length
of the passed vector, if greater than number of bases,
use the last entry for the geometry mapping
---
 src/ASM/FiniteElement.C | 28 ++++++++++++++--------------
 src/ASM/FiniteElement.h |  4 ++++
 2 files changed, 18 insertions(+), 14 deletions(-)

diff --git a/src/ASM/FiniteElement.C b/src/ASM/FiniteElement.C
index ebec2d2d..7f465431 100644
--- a/src/ASM/FiniteElement.C
+++ b/src/ASM/FiniteElement.C
@@ -110,22 +110,22 @@ bool MxFiniteElement::Hessian (Matrix3D& Hess, const Matrix& Jac,
                                const std::vector<const BasisFunctionVals*>* bf,
                                const std::vector<Matrix3D>* d2Nxdu2)
 {
-  bool ok = utl::Hessian(Hess,
-                         gBasis > 1 ? d2MdX2[gBasis-2] : d2NdX2, Jac, Xnod,
-                         bf ? (*bf)[gBasis-1]->d2Ndu2 : (*d2Nxdu2)[gBasis-1],
-                         gBasis > 1 ? dMdX[gBasis-2] : dNdX);
-
-  if (ok && gBasis > 1)
-    ok = utl::Hessian(Hess, d2NdX2, Jac, Xnod,
-                      bf ? bf->front()->d2Ndu2 : d2Nxdu2->front(),
-                      dNdX, false);
-
   size_t nBasis = bf ? bf->size() : d2Nxdu2->size();
-  for (size_t b = 2; b <= nBasis && ok; b++)
-    if (b != gBasis)
-      ok = utl::Hessian(Hess, d2MdX2[b-2], Jac, Xnod,
+  const bool separateGeometry = nBasis > this->getNoBasis();
+  bool ok;
+  if (separateGeometry)
+    ok = Hess.multiply(Xnod, (bf ? bf->back()->d2Ndu2 : d2Nxdu2->back()));
+  else {
+    ok = utl::Hessian(Hess, this->hess(gBasis), Jac, Xnod,
+                      bf ? (*bf)[gBasis-1]->d2Ndu2 : (*d2Nxdu2)[gBasis-1],
+                      this->grad(gBasis));
+  }
+
+  for (size_t b = 1; b <= this->getNoBasis() && ok; b++)
+    if (b != gBasis || separateGeometry)
+      ok = utl::Hessian(Hess, this->hess(b), Jac, Xnod,
                         bf ? (*bf)[b-1]->d2Ndu2 : (*d2Nxdu2)[b-1],
-                        dMdX[b-2], false);
+                        this->grad(b), false);
 
   return ok;
 }
diff --git a/src/ASM/FiniteElement.h b/src/ASM/FiniteElement.h
index 00ff0791..6c420ff0 100644
--- a/src/ASM/FiniteElement.h
+++ b/src/ASM/FiniteElement.h
@@ -51,6 +51,8 @@ public:
 
   //! \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; }
 
@@ -115,6 +117,8 @@ public:
 
   //! \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]; }