Merge branch 'hesseboundary' into 'master'

Hesseboundary

See merge request jschoeberl/ngsolve!279

comp/h1hofespace.cpp

@@ -187,7 +187,8 @@ namespace ngcomp
       case 2:
         additional_evaluators.Set ("hesse", make_shared<T_DifferentialOperator<DiffOpHesse<2>>> ()); break;
       case 3:
-        additional_evaluators.Set ("hesse", make_shared<T_DifferentialOperator<DiffOpHesse<3>>> ()); break;
+        additional_evaluators.Set ("hesse", make_shared<T_DifferentialOperator<DiffOpHesse<3>>> ());
+	additional_evaluators.Set ("hesseboundary", make_shared<T_DifferentialOperator<DiffOpHesseBoundary<3>>> ());break;
       default:
         ;
       }

fem/bdbequations.hpp

@@ -168,6 +168,8 @@ namespace ngfem
 
   };
 
+
+
 
 
 
@@ -535,15 +537,92 @@ namespace ngfem
   };
 
 
+  template <int D, typename FEL = ScalarFiniteElement<D-1> >
+  class DiffOpHesseBoundary : public DiffOp<DiffOpHesseBoundary<D, FEL> >
+  {
+  public:
+    enum { DIM = 1 };
+    enum { DIM_SPACE = D };
+    enum { DIM_ELEMENT = D-1 };
+    enum { DIM_DMAT = D*D };
+    enum { DIFFORDER = 2 };
 
+    static string Name() { return "hesseboundary"; }
+
+    static auto & Cast (const FiniteElement & fel) 
+    { return static_cast<const FEL&> (fel); }
+
+    ///
+    template <typename AFEL, typename MIP, typename MAT>
+    static void GenerateMatrix (const AFEL & fel, const MIP & mip,
+				MAT & mat, LocalHeap & lh)
+    {
+      HeapReset hr(lh);
+      int nd = Cast(fel).GetNDof();
+
+      auto ddshape = Cast(fel).GetDDShape(mip.IP(), lh);
+      auto jinv = mip.GetJacobianInverse();
+      auto dshape = Cast(fel).GetDShape (mip.IP(), lh);
+
 
+      for( int n = 0; n < nd; n++)
+	{
+	  for( int i = 0; i < D; i++)
+	    {
+	      for( int j = 0; j < D; j++)
+		{
+		  mat(i*D+j,n) = 0.0;
+		  for( int k = 0; k < D-1; k++)
+		    {
+		      for( int l = 0; l < D-1; l++)
+			{
+			  mat(i*D+j,n) += jinv(k,i)*ddshape(n,k*(D-1)+l)*jinv(l,j);
+			}
+		    }
+		}		  
+	    }
+	}
+
+      //for non-curved elements, we are finished, otherwise derivatives of Jacobian have to be computed...
+      if (!mip.GetTransformation().IsCurvedElement()) return;
 
+
+      Mat<2,2> hesse[3];
+      mip.CalcHesse (hesse[0], hesse[1], hesse[2]);
 
+      FlatMatrix<> tmp(D, (D-1)*(D-1), lh);
 
+      for (int i = 0; i < D; i++)
+        {
+          for (int j = 0; j < D-1; j++)
+            {
+              for (int k = 0; k < D-1; k++)
+                tmp(i,j*(D-1)+k) = hesse[i](j,k);
+            }
+        }
 
+      FlatMatrix<> tmpmat(nd, (D-1)*(D-1), lh);
+      tmpmat = dshape*jinv*tmp;
 
+      for( int n = 0; n < nd; n++)
+	{
+	  for( int i = 0; i < D; i++)
+	    {
+	      for( int j = 0; j < D; j++)
+		{
+		  for( int k = 0; k < D-1; k++)
+		    {
+		      for( int l = 0; l < D-1; l++)
+			{
+			  mat(i*D+j,n) -= jinv(k,i)*tmpmat(n,k*(D-1)+l)*jinv(l,j);
+			}
+		    }
+		}		  
+	    }
+	}
+    }
 
-
+  };