Modernize LocalIntegralEnergy

dune/elasticity/materials/localintegralenergy.hh

@@ -2,6 +2,7 @@
 #define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
 
 #include <dune/common/fmatrix.hh>
+#include <dune/common/transpose.hh>
 
 #include <dune/geometry/quadraturerules.hh>
 
@@ -53,8 +54,9 @@ energy(const LocalView& localView,
   field_type energy = 0;
 
   // store gradients of shape functions and base functions
-  std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-  std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+  using FiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
+  using Jacobian = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
+  std::vector<Jacobian> jacobians(localFiniteElement.size());
 
   int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
                                                : localFiniteElement.localBasis().order() * gridDim;
@@ -65,23 +67,23 @@ energy(const LocalView& localView,
   {
     const DT integrationElement = element.geometry().integrationElement(qp.position());
 
-    const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
+    const auto geometryJacobianIT = element.geometry().jacobianInverseTransposed(qp.position());
 
     // Global position
     auto x = element.geometry().global(qp.position());
 
     // Get gradients of shape functions
-    localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
+    localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
 
     // compute gradients of base functions
-    for (size_t i=0; i<gradients.size(); ++i)
-      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+    for (size_t i=0; i<jacobians.size(); ++i)
+      jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
 
     // Deformation gradient
     FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
-    for (size_t i=0; i<gradients.size(); i++)
+    for (size_t i=0; i<jacobians.size(); i++)
       for (int j=0; j<gridDim; j++)
-        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], gradients[i]);
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], jacobians[i][0]);
     // Integrate the energy density
     energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
   }