Bugfix hessian

dune/functions/functionspacebases/test/test-bound-localfe.hh

@@ -77,6 +77,9 @@ class ShapeFunctionDerivativeAsCallable
 {
   typedef typename LFE::Traits::LocalBasisType::Traits::DomainType DomainType;
   typedef typename LFE::Traits::LocalBasisType::Traits::RangeType ValueType;
+  typedef typename LFE::Traits::LocalBasisType::Traits::RangeFieldType R;
+  static int constexpr dimRange = LFE::Traits::LocalBasisType::Traits::dimRange;
+
   typedef typename LFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
   static constexpr unsigned int dim = LFE::Traits::LocalBasisType::Traits::dimDomain;
 
@@ -100,7 +103,7 @@ public:
   {
     return [&t](DomainType const &x)
     {
-      Dune::FieldMatrix<ValueType, dim, dim> referenceHessian, hessian;
+      Dune::FieldVector<Dune::FieldMatrix<R, dim, dim>, dimRange> referenceHessian, hessian;
       std::vector<ValueType> partials;
       std::array<unsigned int, dim> dir;
       for (std::size_t i = 0; i < dim; ++i)
@@ -110,11 +113,13 @@ public:
           dir[i]++;
           dir[j]++;
           t.fe_.localBasis().partial(dir, x, partials);
-          referenceHessian[i][j] = partials[t.shapeFunction_];
+          for (std::size_t k = 0; k < dimRange; ++k)
+            referenceHessian[k][i][j] = partials[t.shapeFunction_][k];
         }
       hessian = 0;
       auto JIT = t.e_.geometry().jacobianInverseTransposed(x);
-      hessian = JIT * referenceHessian * transpose(JIT);
+      for (std::size_t k = 0; k < dimRange; ++k)
+        hessian[k] = JIT * referenceHessian[k] * transpose(JIT);
       return hessian;
     };
   }