Migrate from dune-fufem to dune-functions bases

Migrate from dune-fufem to dune-functions bases.

dune/gfe/localintegralenergy.hh

+#ifndef DUNE_GFE_LOCALINTEGRALENERGY_HH
+#define DUNE_GFE_LOCALINTEGRALENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localenergy.hh>
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::GFE {
+
+  /**
+    \brief LocalIntegralEnergy: Assembles the elastic energy for a single element using the localdensity.
+           This class works similarly to the class Dune::Elasticity::LocalIntegralEnergy, where Dune::Elasticity::LocalIntegralEnergy works with
+           a power basis (for the displacement function) and Dune::GFE::LocalIntegralEnergy works with a CompositeBasis (for the displacement function
+           and the rotational field)
+  */
+template<class LocalView, class field_type=double>
+class LocalIntegralEnergy
+  : public Dune::Elasticity::LocalEnergy<LocalView,field_type>
+{
+  using GridView = typename LocalView::GridView;
+  using DT = typename GridView::Grid::ctype;
+
+  enum {gridDim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a local energy density
+    */
+  LocalIntegralEnergy(const std::shared_ptr<Dune::Elasticity::LocalDensity<gridDim,field_type,DT>>& ld)
+  : localDensity_(ld)
+  {}
+
+  /** \brief Virtual destructor */
+  virtual ~LocalIntegralEnergy()
+  {}
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy (const LocalView& localView,
+                     const std::vector<field_type>& localConfiguration) const;
+
+protected:
+  const std::shared_ptr<Dune::Elasticity::LocalDensity<gridDim,field_type,DT>> localDensity_ = nullptr;
+
+};
+
+template <class LocalView, class field_type>
+field_type
+LocalIntegralEnergy<LocalView, field_type>::
+energy(const LocalView& localView,
+       const std::vector<field_type>& localConfiguration) const
+{
+  using namespace Dune::Indices;
+  // powerBasis: grab the finite element of the first child
+  const auto& localFiniteElement = localView.tree().child(_0,0).finiteElement();
+  const auto& element = localView.element();
+
+  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());
+
+  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                               : localFiniteElement.localBasis().order() * gridDim;
+
+  const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+
+  for (const auto& qp : quad)
+  {
+    const DT integrationElement = element.geometry().integrationElement(qp.position());
+
+    const auto jacobianInverseTransposed = 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);
+
+    // compute gradients of base functions
+    for (size_t i=0; i<gradients.size(); ++i)
+      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+    // Deformation gradient
+    FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
+    for (size_t i=0; i<gradients.size(); i++)
+      for (int j=0; j<gridDim; j++)
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(_0,j).localIndex(i) ], gradients[i]);
+    // Integrate the energy density
+    energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
+  }
+
+  return energy;
+}
+
+}  // namespace Dune::GFE
+
+#endif   //#ifndef DUNE_GFE_LOCALINTEGRALENERGY_HH