diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000000000000000000000000000000000000..93f9ddf1880a4ba9a5dd845abcabd872ce3b29f6
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,8 @@
+# Master (will become release 2.8)
+
+- Introduce class `LocalDensity` for material-specific implementations
+- Introduce class `LocalIntegralEnergy` to work with the densities
+
+## Deprecations and removals
+
+- Redundant implementations of `LocalEnergy` classes are now replaced by combining `LocalDensity` and `LocalIntegralEnergy`
diff --git a/dune/elasticity/materials/CMakeLists.txt b/dune/elasticity/materials/CMakeLists.txt
index 2255928d0944a70b17845df20c5385193dc2dec9..fced7f1f79078f0ed1812d72112c3e7794d7587e 100644
--- a/dune/elasticity/materials/CMakeLists.txt
+++ b/dune/elasticity/materials/CMakeLists.txt
@@ -1,13 +1,21 @@
 install(FILES
     adolcmaterial.hh
     adolcneohookeanmaterial.hh
+    exphenckydensity.hh
     exphenckyenergy.hh
     geomexactstvenantkirchhoffmaterial.hh
+    henckydensity.hh
     henckyenergy.hh
+    localdensity.hh
+    localintegralenergy.hh
+    mooneyrivlindensity.hh
+    mooneyrivlinenergy.hh
     neohookeanmaterial.hh
+    neohookedensity.hh
     neohookeenergy.hh
     neumannenergy.hh
     ogdenmaterial.hh
+    stvenantkirchhoffdensity.hh
     stvenantkirchhoffenergy.hh
     sumenergy.hh
     DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/elasticity/materials)
diff --git a/dune/elasticity/materials/exphenckydensity.hh b/dune/elasticity/materials/exphenckydensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b20a440fb061768f3c2cc6216dddeb23b0eedb6a
--- /dev/null
+++ b/dune/elasticity/materials/exphenckydensity.hh
@@ -0,0 +1,95 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYDENSITY_HH
+#define DUNE_ELASTICITY_MATERIALS_EXPHENCKYDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::Elasticity {
+
+template<int dim, class field_type = double>
+class ExpHenckyDensity final
+  : public LocalDensity<dim,field_type>
+{
+public:
+
+  /** \brief Constructor with a set of material parameters
+  * \param parameters The material parameters
+  */
+  ExpHenckyDensity(const Dune::ParameterTree& parameters)
+  {
+    // Lame constants
+    mu_     = parameters.get<double>("mu");
+    lambda_ = parameters.get<double>("lambda");
+    kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
+
+    // Exponential Hencky constants with naming convention from Neff, Ghiba and Lankeit,
+    // "The exponentiated Hencky-logarithmic strain energy. Part I: Constitutive issues
+    //  and rank one convexity"
+    k_      = parameters.get<double>("k");
+    khat_   = parameters.get<double>("khat");
+
+    // weights for distortion and dilatation parts of the energy
+    alpha_  = mu_ / k_;
+    beta_   = kappa_ / (2.0 * khat_);
+  }
+
+  /** \brief Evaluation with the deformation gradient
+  *
+  * \param x Position of the gradient
+  * \param gradient Deformation gradient
+  */
+  field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
+  {
+    /////////////////////////////////////////////////////////
+    // compute C = F^TF
+    /////////////////////////////////////////////////////////
+
+    Dune::FieldMatrix<field_type,dim,dim> C(0);
+    for (int i=0; i<dim; i++)
+      for (int j=0; j<dim; j++)
+        for (int k=0; k<dim; k++)
+          C[i][j] += gradient[k][i] * gradient[k][j];
+
+    //////////////////////////////////////////////////////////
+    //  Eigenvalues of FTF
+    //////////////////////////////////////////////////////////
+
+    // compute eigenvalues of C, i.e., squared singular values \sigma_i of F
+    Dune::FieldVector<field_type, dim> sigmaSquared;
+    FMatrixHelp::eigenValues(C, sigmaSquared);
+
+    // logarithms of the singular values of F, i.e., eigenvalues of U
+    std::array<field_type, dim> lnSigma;
+    for (int i = 0; i < dim; i++)
+      lnSigma[i] = 0.5 * log(sigmaSquared[i]);
+
+    field_type trLogU = 0;
+    for (int i = 0; i < dim; i++)
+      trLogU += lnSigma[i];
+
+    field_type normDevLogUSquared = 0;
+    for (int i = 0; i < dim; i++)
+    {
+      // the deviator shifts the
+      field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
+      normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
+    }
+
+    // Add the local energy density
+    auto distortionEnergy = alpha_ * exp(khat_ * normDevLogUSquared);
+    auto dilatationEnergy = beta_  * exp(k_ * (trLogU * trLogU));
+
+    return distortionEnergy + dilatationEnergy;
+  }
+
+  /** \brief Lame constants */
+  field_type mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
+};
+
+} // namespace Dune::Elasticity
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGYDENSITY_HH
+
+
diff --git a/dune/elasticity/materials/exphenckyenergy.hh b/dune/elasticity/materials/exphenckyenergy.hh
index 3fd3ce61655667f05e6ac3ce0807cb3acd69c2af..e3a9f5be2c1b8321303c43f70a319365814346e6 100644
--- a/dune/elasticity/materials/exphenckyenergy.hh
+++ b/dune/elasticity/materials/exphenckyenergy.hh
@@ -1,148 +1,26 @@
 #ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
 #define DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_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/exphenckydensity.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
 
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class ExpHenckyEnergy
-  : public Elasticity::LocalEnergy<GridView,
-                            LocalFiniteElement,
-                            std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
+class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::ExpHenckyDensity")
+ExpHenckyEnergy
+  : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {
-  // grid types
-  typedef typename GridView::Grid::ctype DT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
-
-  // some other sizes
-  enum {gridDim=GridView::dimension};
-  enum {dim=GridView::dimension};
 
+  using Base = Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>;
 public:
 
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
+  // backwards compatibility: forward directly to integral energy
   ExpHenckyEnergy(const Dune::ParameterTree& parameters)
-  {
-    // Lame constants
-    mu_     = parameters.get<double>("mu");
-    lambda_ = parameters.get<double>("lambda");
-    kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
-
-    // Exponential Hencky constants with naming convention from Neff, Ghiba and Lankeit,
-    // "The exponentiated Hencky-logarithmic strain energy. Part I: Constitutive issues
-    //  and rank one convexity"
-    k_      = parameters.get<double>("k");
-    khat_   = parameters.get<double>("khat");
-
-    // weights for distortion and dilatation parts of the energy
-    alpha_  = mu_ / k_;
-    beta_   = kappa_ / (2.0 * khat_);
-  }
-
-  /** \brief Assemble the energy for a single element */
-  field_type energy (const Entity& e,
-                     const LocalFiniteElement& localFiniteElement,
-                     const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
-
-  /** \brief Lame constants */
-  field_type mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
+  : Base(std::make_shared<Elasticity::ExpHenckyDensity<GridView::dimension,field_type>>(parameters))
+  {}
 };
 
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-ExpHenckyEnergy<GridView, LocalFiniteElement, field_type>::
-energy(const Entity& element,
-       const LocalFiniteElement& localFiniteElement,
-       const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
-{
-  assert(element.type() == localFiniteElement.type());
-  typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
-  // store gradients of shape functions and base functions
-  std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-  std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
-
-  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                               : localFiniteElement.localBasis().order() * gridDim;
-
-  const Dune::QuadratureRule<DT, gridDim>& quad
-      = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
-  field_type distortionEnergy = 0, dilatationEnergy = 0;
-  for (size_t pt=0; pt<quad.size(); pt++)
-  {
-    // Local position of the quadrature point
-    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
-    const DT integrationElement = element.geometry().integrationElement(quadPos);
-
-    const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
-    DT weight = quad[pt].weight() * integrationElement;
-
-    // get gradients of shape functions
-    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-
-    // compute gradients of base functions
-    for (size_t i=0; i<gradients.size(); ++i)
-      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-    Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
-    for (size_t i=0; i<gradients.size(); i++)
-      for (int j=0; j<gridDim; j++)
-        derivative[j].axpy(localConfiguration[i][j], gradients[i]);
-
-    /////////////////////////////////////////////////////////
-    // compute C = F^T F
-    /////////////////////////////////////////////////////////
-
-    Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
-    for (int i=0; i<gridDim; i++)
-      for (int j=0; j<gridDim; j++)
-        for (int k=0; k<gridDim; k++)
-          C[i][j] += derivative[k][i] * derivative[k][j];
-
-    //////////////////////////////////////////////////////////
-    //  Eigenvalues of FTF
-    //////////////////////////////////////////////////////////
-
-    // compute eigenvalues of C, i.e., squared singular values \sigma_i of F
-    Dune::FieldVector<field_type, dim> sigmaSquared;
-    FMatrixHelp::eigenValues(C, sigmaSquared);
-
-    // logarithms of the singular values of F, i.e., eigenvalues of U
-    std::array<field_type, dim> lnSigma;
-    for (int i = 0; i < dim; i++)
-      lnSigma[i] = 0.5 * log(sigmaSquared[i]);
-
-    field_type trLogU = 0;
-    for (int i = 0; i < dim; i++)
-      trLogU += lnSigma[i];
-
-    field_type normDevLogUSquared = 0;
-    for (int i = 0; i < dim; i++)
-    {
-      // the deviator shifts the
-      field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
-      normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
-    }
-
-    // Add the local energy density
-    distortionEnergy += weight * alpha_ * exp(khat_ * normDevLogUSquared);
-    dilatationEnergy += weight * beta_  * exp(k_ * (trLogU * trLogU));
-  }
-
-  return distortionEnergy + dilatationEnergy;
-}
-
 }  // namespace Dune
 
 #endif   //#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
diff --git a/dune/elasticity/materials/henckydensity.hh b/dune/elasticity/materials/henckydensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..81a69562b711dd280c96a5f6e88ee0f46d5d1acd
--- /dev/null
+++ b/dune/elasticity/materials/henckydensity.hh
@@ -0,0 +1,84 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYDENSITY_HH
+#define DUNE_ELASTICITY_MATERIALS_HENCKYDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::Elasticity {
+
+template<int dim, class field_type = double>
+class HenckyDensity final
+  : public LocalDensity<dim,field_type>
+{
+public:
+  /** \brief Constructor with a set of material parameters
+    *
+    * \param parameters The material parameters
+    */
+  HenckyDensity(const Dune::ParameterTree& parameters)
+  {
+    // Lame constants
+    mu_ = parameters.template get<double>("mu");
+    lambda_ = parameters.template get<double>("lambda");
+    kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
+  }
+
+  /** \brief Evaluation with the deformation gradient
+  *
+  * \param x Position of the gradient
+  * \param gradient Deformation gradient
+  */
+  field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
+  {
+    /////////////////////////////////////////////////////////
+    // compute C = F^TF
+    /////////////////////////////////////////////////////////
+
+    Dune::FieldMatrix<field_type,dim,dim> C(0);
+    for (int i=0; i<dim; i++)
+      for (int j=0; j<dim; j++)
+        for (int k=0; k<dim; k++)
+          C[i][j] += gradient[k][i] * gradient[k][j];
+
+    //////////////////////////////////////////////////////////
+    //  Eigenvalues of C = F^TF
+    //////////////////////////////////////////////////////////
+
+    Dune::FieldVector<field_type,dim> sigmaSquared;
+    FMatrixHelp::eigenValues(C, sigmaSquared);
+
+    // logarithms of the singular values of F, i.e., eigenvalues of U
+    std::array<field_type, dim> lnSigma;
+    for (int i = 0; i < dim; i++)
+      lnSigma[i] = 0.5 * log(sigmaSquared[i]);
+
+    field_type trLogU = 0;
+    for (int i = 0; i < dim; i++)
+      trLogU += lnSigma[i];
+
+    field_type normDevLogUSquared = 0;
+    for (int i = 0; i < dim; i++)
+    {
+      // the deviator shifts the spectrum by the trace
+      field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
+      normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
+    }
+
+    // Add the local energy density
+    auto distortionEnergy = mu_ * normDevLogUSquared;
+    auto dilatationEnergy = 0.5 * kappa_ * (trLogU * trLogU);
+
+    return dilatationEnergy + distortionEnergy;
+  }
+
+  /** \brief Lame constants */
+  double mu_, lambda_, kappa_;
+};
+
+} // namespace Dune::Elasticity
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYDENSITY_HH
+
+
diff --git a/dune/elasticity/materials/henckyenergy.hh b/dune/elasticity/materials/henckyenergy.hh
index 830bb909107ede74523b6863f603b0895a15d0c0..e85f6dfeb483ea8f3a98638d8c10886aee274806 100644
--- a/dune/elasticity/materials/henckyenergy.hh
+++ b/dune/elasticity/materials/henckyenergy.hh
@@ -1,137 +1,26 @@
 #ifndef DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
 #define DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_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/henckydensity.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
 
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class HenckyEnergy
-    : public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension> > >
+class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::HenckyDensity")
+HenckyEnergy
+  : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {
-    // grid types
-    typedef typename GridView::Grid::ctype DT;
-    typedef typename GridView::template Codim<0>::Entity Entity;
-
-    // some other sizes
-    enum {gridDim=GridView::dimension};
-    enum {dim=GridView::dimension};
-
-public:  // for testing
 
-    /** \brief Constructor with a set of material parameters
-     * \param parameters The material parameters
-     */
-    HenckyEnergy(const Dune::ParameterTree& parameters)
-    {
-      // Lame constants
-      mu_ = parameters.template get<double>("mu");
-      lambda_ = parameters.template get<double>("lambda");
-      kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
-    }
+  using Base = Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>;
+public:
 
-    /** \brief Assemble the energy for a single element */
-    field_type energy (const Entity& e,
-               const LocalFiniteElement& localFiniteElement,
-               const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
-
-    /** \brief Lame constants */
-    double mu_, lambda_, kappa_;
+  // backwards compatibility: forward directly to integral energy
+  HenckyEnergy(const Dune::ParameterTree& parameters)
+  : Base(std::make_shared<Elasticity::HenckyDensity<GridView::dimension,field_type>>(parameters))
+  {}
 };
 
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-HenckyEnergy<GridView,LocalFiniteElement,field_type>::
-energy(const Entity& element,
-       const LocalFiniteElement& localFiniteElement,
-       const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
-{
-    assert(element.type() == localFiniteElement.type());
-    typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
-    field_type energy = 0;
-
-    // store gradients of shape functions and base functions
-    std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-    std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
-
-    int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                                 : localFiniteElement.localBasis().order() * gridDim;
-
-    const Dune::QuadratureRule<DT, gridDim>& quad
-        = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
-    field_type distortionEnergy = 0, dilatationEnergy = 0;
-    for (size_t pt=0; pt<quad.size(); pt++) {
-
-        // Local position of the quadrature point
-        const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
-        const DT integrationElement = element.geometry().integrationElement(quadPos);
-
-        const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
-        DT weight = quad[pt].weight() * integrationElement;
-
-        // get gradients of shape functions
-        localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-
-        // compute gradients of base functions
-        for (size_t i=0; i<gradients.size(); ++i)
-          jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-        Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
-        for (size_t i=0; i<gradients.size(); i++)
-          for (int j=0; j<gridDim; j++)
-            derivative[j].axpy(localConfiguration[i][j], gradients[i]);
-
-        /////////////////////////////////////////////////////////
-        // compute C = F^TF
-        /////////////////////////////////////////////////////////
-
-        Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
-        for (int i=0; i<gridDim; i++)
-          for (int j=0; j<gridDim; j++)
-            for (int k=0; k<gridDim; k++)
-              C[i][j] += derivative[k][i] * derivative[k][j];
-
-        //////////////////////////////////////////////////////////
-        //  Eigenvalues of C = F^TF
-        //////////////////////////////////////////////////////////
-
-        Dune::FieldVector<field_type,dim> sigmaSquared;
-        FMatrixHelp::eigenValues(C, sigmaSquared);
-
-        // logarithms of the singular values of F, i.e., eigenvalues of U
-        std::array<field_type, dim> lnSigma;
-        for (int i = 0; i < dim; i++)
-          lnSigma[i] = 0.5 * log(sigmaSquared[i]);
-
-        field_type trLogU = 0;
-        for (int i = 0; i < dim; i++)
-          trLogU += lnSigma[i];
-
-        field_type normDevLogUSquared = 0;
-        for (int i = 0; i < dim; i++)
-        {
-          // the deviator shifts the spectrum by the trace
-          field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
-          normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
-        }
-
-        // Add the local energy density
-        distortionEnergy += weight * mu_ * normDevLogUSquared;
-        dilatationEnergy += weight * 0.5 * kappa_ * (trLogU * trLogU);
-    }
-
-    return distortionEnergy + dilatationEnergy;
-}
-
 }  // namespace Dune
 
 #endif   //#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
diff --git a/dune/elasticity/materials/localdensity.hh b/dune/elasticity/materials/localdensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c091b4d7f84961c6e93c4c6870ae7c8d1a328f7e
--- /dev/null
+++ b/dune/elasticity/materials/localdensity.hh
@@ -0,0 +1,26 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_LOCALDENSITY_HH
+#define DUNE_ELASTICITY_MATERIALS_LOCALDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+
+namespace Dune::Elasticity {
+
+/** \brief A base class for energy densities to be evaluated in an integral energy */
+template<int dim, class field_type = double>
+class LocalDensity
+{
+
+public:
+
+  /** \brief Evaluation with the deformation gradient
+   *
+   * \param x Position of the gradient
+   * \param gradient Deformation gradient
+   */
+  virtual field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const = 0;
+
+};
+
+}  // namespace Dune::Elasticity
+
+#endif  // DUNE_ELASTICITY_MATERIALS_LOCALDENSITY_HH
diff --git a/dune/elasticity/materials/localintegralenergy.hh b/dune/elasticity/materials/localintegralenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1307c2503ed2e882538987db600e9b8df4e6307f
--- /dev/null
+++ b/dune/elasticity/materials/localintegralenergy.hh
@@ -0,0 +1,99 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_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::Elasticity {
+
+template<class GridView, class LocalFiniteElement, class field_type=double>
+class LocalIntegralEnergy
+  : public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension>>>
+{
+  // grid types
+  using DT = typename GridView::Grid::ctype;
+  using Entity = typename GridView::template Codim<0>::Entity;
+
+  // some other sizes
+  enum {gridDim=GridView::dimension};
+  enum {dim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a local energy density
+    */
+  LocalIntegralEnergy(const std::shared_ptr<LocalDensity<dim,field_type>>& ld)
+  : localDensity_(ld)
+  {}
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy (const Entity& e,
+                     const LocalFiniteElement& localFiniteElement,
+                     const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
+
+protected:
+  const std::shared_ptr<LocalDensity<dim,field_type>> localDensity_ = nullptr;
+
+};
+
+
+
+template <class GridView, class LocalFiniteElement, class field_type>
+field_type
+LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>::
+energy(const Entity& element,
+       const LocalFiniteElement& localFiniteElement,
+       const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
+{
+  assert(element.type() == localFiniteElement.type());
+  using Geometry = typename GridView::template Codim<0>::Entity::Geometry;
+
+  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[i][j], gradients[i]);
+
+    // Integrate the energy density
+    energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
+  }
+
+  return energy;
+}
+
+}  // namespace Dune::Elasticity
+
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
diff --git a/dune/elasticity/materials/mooneyrivlindensity.hh b/dune/elasticity/materials/mooneyrivlindensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4dcacdcde00188d0af641449690e09371149f99d
--- /dev/null
+++ b/dune/elasticity/materials/mooneyrivlindensity.hh
@@ -0,0 +1,132 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINDENSITY_HH
+#define DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::Elasticity {
+
+template<int dim, class field_type = double>
+class MooneyRivlinDensity final
+    : public LocalDensity<dim,field_type>
+{
+public:
+
+  /** \brief Constructor with a set of material parameters
+  * \param parameters The material parameters
+  */
+  MooneyRivlinDensity(const Dune::ParameterTree& parameters)
+  {
+    // Mooneyrivlin constants
+    mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
+    mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
+    mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
+
+    mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10");
+    mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01");
+    mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20");
+    mooneyrivlin_02 = parameters.get<double>("mooneyrivlin_02");
+    mooneyrivlin_11 = parameters.get<double>("mooneyrivlin_11");
+    mooneyrivlin_30 = parameters.get<double>("mooneyrivlin_30");
+    mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03");
+    mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21");
+    mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12");
+    mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k");
+
+    mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
+  }
+
+  /** \brief Evaluation with the deformation gradient
+  *
+  * \param x Position of the gradient
+  * \param gradient Deformation gradient
+  */
+  field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
+  {
+    /////////////////////////////////////////////////////////
+    // compute C = F^TF
+    /////////////////////////////////////////////////////////
+
+    Dune::FieldMatrix<field_type,dim,dim> C(0);
+    for (int i=0; i<dim; i++)
+      for (int j=0; j<dim; j++)
+        for (int k=0; k<dim; k++)
+          C[i][j] += gradient[k][i] * gradient[k][j];
+
+    //////////////////////////////////////////////////////////
+    //  Eigenvalues of C = F^TF
+    //////////////////////////////////////////////////////////
+
+
+    // eigenvalues of C, i.e., squared singular values \sigma_i of F
+    Dune::FieldVector<field_type, dim> sigmaSquared;
+    FMatrixHelp::eigenValues(C, sigmaSquared);
+
+    field_type normFSquared = gradient.frobenius_norm2();
+    field_type detF = gradient.determinant();
+
+    field_type normFinvSquared = 0;
+
+    field_type c2Tilde = 0;
+    for (int i = 0; i < dim; i++) {
+      normFinvSquared += 1/sigmaSquared[i];
+      // compute D, which is the sum of the squared eigenvalues
+      for (int j = i+1; j < dim; j++)
+        c2Tilde += sigmaSquared[j]*sigmaSquared[i];
+      }
+
+    field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
+    // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
+    c2Tilde /= pow(detF, 4.0/dim);
+    field_type c2TildeMinus3 = c2Tilde - 3;
+    field_type detFMinus1 = detF - 1;
+
+    // Add the local energy density
+    auto strainEnergyCiarlet = (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF));
+    auto strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
+                        mooneyrivlin_01 * c2TildeMinus3  +
+                        mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 +
+                        mooneyrivlin_02 * c2TildeMinus3  * c2TildeMinus3  +
+                        mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3  +
+                        mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 +
+                        mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3  +
+                        mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3  * c2TildeMinus3  +
+                        mooneyrivlin_03 * c2TildeMinus3  * c2TildeMinus3  * c2TildeMinus3;
+    field_type logDetF = std::log(detF);
+    auto strainEnergyWithLog =  ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF );
+    auto strainEnergyWithSquare =  ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1);
+
+
+    std::cout << std::scientific;
+    if (mooneyrivlin_energy == "log") {
+      return strainEnergyWithLog;
+    } else if (mooneyrivlin_energy == "square") {
+      return strainEnergyWithSquare;
+    }
+    std::cout << std::fixed;
+
+    return strainEnergyCiarlet;
+  }
+
+  /** \brief Lame constants */
+  field_type mooneyrivlin_a,
+             mooneyrivlin_b,
+             mooneyrivlin_c,
+             mooneyrivlin_10,
+             mooneyrivlin_01,
+             mooneyrivlin_20,
+             mooneyrivlin_02,
+             mooneyrivlin_11,
+             mooneyrivlin_30,
+             mooneyrivlin_21,
+             mooneyrivlin_12,
+             mooneyrivlin_03,
+             mooneyrivlin_k;
+  std::string mooneyrivlin_energy;
+};
+
+} // namespace Dune::Elasticity
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINDENSITY_HH
diff --git a/dune/elasticity/materials/mooneyrivlinenergy.hh b/dune/elasticity/materials/mooneyrivlinenergy.hh
index 508112e0b8f59644d05fa60f8eb1fd61dd95f690..d24d8e24c1340468a54031fea104b0c2cdccb247 100644
--- a/dune/elasticity/materials/mooneyrivlinenergy.hh
+++ b/dune/elasticity/materials/mooneyrivlinenergy.hh
@@ -1,192 +1,26 @@
 #ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH
 #define DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH
 
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
+#include <dune/elasticity/materials/mooneyrivlindensity.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
 
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class MooneyRivlinEnergy
-  : public Elasticity::LocalEnergy<GridView,
-                            LocalFiniteElement,
-                            std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
+class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::MooneyRivlinDensity")
+MooneyRivlinEnergy
+  : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {
-  // grid types
-  typedef typename GridView::Grid::ctype DT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
-
-  // some other sizes
-  enum {gridDim=GridView::dimension};
-  enum {dim=GridView::dimension};
 
+  using Base = Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>;
 public:
 
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
+  // backwards compatibility: forward directly to integral energy
   MooneyRivlinEnergy(const Dune::ParameterTree& parameters)
-  {
-    // Mooneyrivlin constants
-    mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
-    mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
-    mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
-    
-    mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10");
-    mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01");
-    mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20");
-    mooneyrivlin_02 = parameters.get<double>("mooneyrivlin_02");
-    mooneyrivlin_11 = parameters.get<double>("mooneyrivlin_11");
-    mooneyrivlin_30 = parameters.get<double>("mooneyrivlin_30");
-    mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03");
-    mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21");
-    mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12");
-    mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k");
-    
-    mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
-
-  }
-
-  /** \brief Assemble the energy for a single element */
-  field_type energy (const Entity& e,
-                     const LocalFiniteElement& localFiniteElement,
-                     const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
-
-  /** \brief Lame constants */
-  field_type mooneyrivlin_a,
-              mooneyrivlin_b,
-              mooneyrivlin_c,
-              mooneyrivlin_10,
-              mooneyrivlin_01,
-              mooneyrivlin_20,
-              mooneyrivlin_02,
-              mooneyrivlin_11, 
-              mooneyrivlin_30, 
-              mooneyrivlin_21, 
-              mooneyrivlin_12, 
-              mooneyrivlin_03,
-              mooneyrivlin_k;
-  std::string mooneyrivlin_energy;
+  : Base(std::make_shared<Elasticity::MooneyRivlinDensity<GridView::dimension,field_type>>(parameters))
+  {}
 };
 
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-MooneyRivlinEnergy<GridView, LocalFiniteElement, field_type>::
-energy(const Entity& element,
-       const LocalFiniteElement& localFiniteElement,
-       const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
-{
-
-  assert(element.type() == localFiniteElement.type());
-  typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
-  field_type strainEnergyCiarlet = 0;
-  field_type strainEnergy = 0;
-  field_type strainEnergyWithLog = 0;
-  field_type strainEnergyWithSquare = 0;
-
-  // store gradients of shape functions and base functions
-  std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-  std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
-
-  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                               : localFiniteElement.localBasis().order() * gridDim;
-
-  const Dune::QuadratureRule<DT, gridDim>& quad
-      = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
-  for (size_t pt=0; pt<quad.size(); pt++)
-  {
-    // Local position of the quadrature point
-    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
-    const DT integrationElement = element.geometry().integrationElement(quadPos);
-
-    const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
-    DT weight = quad[pt].weight() * integrationElement;
-
-    // get gradients of shape functions
-    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-
-    // compute gradients of base functions
-    for (size_t i=0; i<gradients.size(); ++i)
-      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-    Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
-    for (size_t i=0; i<gradients.size(); i++)
-      for (int j=0; j<gridDim; j++)
-        derivative[j].axpy(localConfiguration[i][j], gradients[i]);
-
-
-    Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
-
-    for (int i=0; i<gridDim; i++) {
-      for (int j=0; j<gridDim; j++) {
-        for (int k=0; k<gridDim; k++)
-          C[i][j] += derivative[k][i] * derivative[k][j];
-      }
-    }
-
-    //////////////////////////////////////////////////////////
-    //  Eigenvalues of FTF
-    //////////////////////////////////////////////////////////
-
-    // eigenvalues of C, i.e., squared singular values \sigma_i of F
-    Dune::FieldVector<field_type, dim> sigmaSquared;
-    FMatrixHelp::eigenValues(C, sigmaSquared);
-
-    field_type normFSquared = derivative.frobenius_norm2();
-    field_type detF = derivative.determinant();
-
-    field_type normFinvSquared = 0;
-
-    field_type c2Tilde = 0;
-    for (int i = 0; i < dim; i++) {
-      normFinvSquared += 1/sigmaSquared[i];
-      // compute D, which is the sum of the squared eigenvalues  
-      for (int j = i+1; j < dim; j++)
-        c2Tilde += sigmaSquared[j]*sigmaSquared[i];
-    }
-
-    field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
-    // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
-    c2Tilde /= pow(detF, 4.0/dim);
-    field_type c2TildeMinus3 = c2Tilde - 3;
-    field_type detFMinus1 = detF - 1; 
-
-    // Add the local energy density
-    strainEnergyCiarlet += weight * (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF));
-    strainEnergy = 0;
-    strainEnergy = mooneyrivlin_10 * trCTildeMinus3 + 
-                  mooneyrivlin_01 * c2TildeMinus3 + 
-                  mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 + 
-                  mooneyrivlin_02 * c2TildeMinus3 * c2TildeMinus3 + 
-                  mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3 +
-                  mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 + 
-                  mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3 + 
-                  mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3 * c2TildeMinus3 + 
-                  mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3;
-    field_type logDetF = std::log(detF);
-    strainEnergyWithLog += weight * ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF );
-    strainEnergyWithSquare += weight * ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1);
-  }
-
-  std::cout << std::scientific;
-  if (mooneyrivlin_energy == "log") {
-    return strainEnergyWithLog;
-  } else if (mooneyrivlin_energy == "square") {
-    return strainEnergyWithSquare;
-  }
-  std::cout << std::fixed;
-
-  return strainEnergyCiarlet;
-}
-
 }  // namespace Dune
 
 #endif   //#ifndef DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINENERGY_HH
diff --git a/dune/elasticity/materials/neohookedensity.hh b/dune/elasticity/materials/neohookedensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..47aa6d61ab2f71f6ba7df83e4ec9f64264522baa
--- /dev/null
+++ b/dune/elasticity/materials/neohookedensity.hh
@@ -0,0 +1,82 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_NEOHOOKEDENSITY_HH
+#define DUNE_ELASTICITY_MATERIALS_NEOHOOKEDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::Elasticity {
+
+template<int dim, class field_type = double>
+class NeoHookeDensity final
+  : public LocalDensity<dim,field_type>
+{
+public:
+
+  /** \brief Constructor with a set of material parameters
+  * \param parameters The material parameters
+  */
+  NeoHookeDensity(const Dune::ParameterTree& parameters)
+  {
+      // Lame constants
+      mu_     = parameters.get<double>("mu");
+      lambda_ = parameters.get<double>("lambda");
+      kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
+  }
+
+  /** \brief Evaluation with the deformation gradient
+  *
+  * \param x Position of the gradient
+  * \param gradient Deformation gradient
+  */
+  field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
+  {
+    /////////////////////////////////////////////////////////
+    // compute C = F^TF
+    /////////////////////////////////////////////////////////
+
+    Dune::FieldMatrix<field_type,dim,dim> C(0);
+    for (int i=0; i<dim; i++)
+      for (int j=0; j<dim; j++)
+        for (int k=0; k<dim; k++)
+          C[i][j] += gradient[k][i] * gradient[k][j];
+
+    //////////////////////////////////////////////////////////
+    //  Eigenvalues of FTF
+    //////////////////////////////////////////////////////////
+
+    // eigenvalues of C, i.e., squared singular values \sigma_i of F
+    Dune::FieldVector<field_type, dim> sigmaSquared;
+    FMatrixHelp::eigenValues(C, sigmaSquared);
+
+    // singular values of F, i.e., eigenvalues of U
+    std::array<field_type, dim> sigma;
+    for (int i = 0; i < dim; i++)
+      sigma[i] = std::sqrt(sigmaSquared[i]);
+
+    field_type detC = 1.0;
+    for (int i = 0; i < dim; i++)
+      detC *= sigmaSquared[i];
+    field_type detF = std::sqrt(detC);
+
+    // \tilde{C} = \tilde{F}^T\tilde{F} = \frac{1}{\det{F}^{2/3}}C
+    field_type trCTilde = 0;
+    for (int i = 0; i < dim; i++)
+      trCTilde += sigmaSquared[i];
+    trCTilde /= pow(detC, 1.0/3.0);
+
+    // Add the local energy density
+    auto distortionEnergy = (0.5 * mu_)    * (trCTilde - 3);
+    auto dilatationEnergy = (0.5 * kappa_) * ((detF - 1) * (detF - 1));
+
+    return distortionEnergy + dilatationEnergy;
+  }
+
+  /** \brief Lame constants */
+  field_type mu_, lambda_, kappa_;
+};
+
+} // namespace Dune::Elasticity
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEOHOOKEDENSITY_HH
diff --git a/dune/elasticity/materials/neohookeenergy.hh b/dune/elasticity/materials/neohookeenergy.hh
index a7a43325a9f05c853818d3733991764074b671b0..f45fa059d088ff42ac89675b5cea4f7ff5c4bd54 100644
--- a/dune/elasticity/materials/neohookeenergy.hh
+++ b/dune/elasticity/materials/neohookeenergy.hh
@@ -1,138 +1,26 @@
 #ifndef DUNE_ELASTICITY_MATERIALS_NEOHOOKEENERGY_HH
 #define DUNE_ELASTICITY_MATERIALS_NEOHOOKEENERGY_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/neohookedensity.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
 
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeoHookeEnergy
-  : public Elasticity::LocalEnergy<GridView,
-                            LocalFiniteElement,
-                            std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
+class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::NeoHookeDensity")
+NeoHookeEnergy
+  : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {
-  // grid types
-  typedef typename GridView::Grid::ctype DT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
-
-  // some other sizes
-  enum {gridDim=GridView::dimension};
-  enum {dim=GridView::dimension};
 
+  using Base = Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>;
 public:
 
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
+  // backwards compatibility: forward directly to integral energy
   NeoHookeEnergy(const Dune::ParameterTree& parameters)
-  {
-    // Lame constants
-    mu_     = parameters.get<double>("mu");
-    lambda_ = parameters.get<double>("lambda");
-    kappa_  = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
-  }
-
-  /** \brief Assemble the energy for a single element */
-  field_type energy (const Entity& e,
-                     const LocalFiniteElement& localFiniteElement,
-                     const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
-
-  /** \brief Lame constants */
-  field_type mu_, lambda_, kappa_;
+  : Base(std::make_shared<Elasticity::NeoHookeDensity<GridView::dimension,field_type>>(parameters))
+  {}
 };
 
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeoHookeEnergy<GridView, LocalFiniteElement, field_type>::
-energy(const Entity& element,
-       const LocalFiniteElement& localFiniteElement,
-       const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
-{
-  assert(element.type() == localFiniteElement.type());
-  typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
-  field_type distortionEnergy = 0, dilatationEnergy = 0;
-
-  // store gradients of shape functions and base functions
-  std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-  std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
-
-  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                               : localFiniteElement.localBasis().order() * gridDim;
-
-  const Dune::QuadratureRule<DT, gridDim>& quad
-      = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
-  for (size_t pt=0; pt<quad.size(); pt++)
-  {
-    // Local position of the quadrature point
-    const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
-    const DT integrationElement = element.geometry().integrationElement(quadPos);
-
-    const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
-    DT weight = quad[pt].weight() * integrationElement;
-
-    // get gradients of shape functions
-    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-
-    // compute gradients of base functions
-    for (size_t i=0; i<gradients.size(); ++i)
-      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-    Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
-    for (size_t i=0; i<gradients.size(); i++)
-      for (int j=0; j<gridDim; j++)
-        derivative[j].axpy(localConfiguration[i][j], gradients[i]);
-
-    /////////////////////////////////////////////////////////
-    // compute C = F^T F
-    /////////////////////////////////////////////////////////
-
-    Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
-    for (int i=0; i<gridDim; i++)
-      for (int j=0; j<gridDim; j++)
-        for (int k=0; k<gridDim; k++)
-          C[i][j] += derivative[k][i] * derivative[k][j];
-
-    //////////////////////////////////////////////////////////
-    //  Eigenvalues of FTF
-    //////////////////////////////////////////////////////////
-
-    // eigenvalues of C, i.e., squared singular values \sigma_i of F
-    Dune::FieldVector<field_type, dim> sigmaSquared;
-    FMatrixHelp::eigenValues(C, sigmaSquared);
-
-    // singular values of F, i.e., eigenvalues of U
-    std::array<field_type, dim> sigma;
-    for (int i = 0; i < dim; i++)
-      sigma[i] = std::sqrt(sigmaSquared[i]);
-
-    field_type detC = 1.0;
-    for (int i = 0; i < dim; i++)
-      detC *= sigmaSquared[i];
-    field_type detF = std::sqrt(detC);
-
-    // \tilde{C} = \tilde{F}^T\tilde{F} = \frac{1}{\det{F}^{2/3}}C
-    field_type trCTilde = 0;
-    for (int i = 0; i < dim; i++)
-      trCTilde += sigmaSquared[i];
-    trCTilde /= pow(detC, 1.0/3.0);
-
-    // Add the local energy density
-    distortionEnergy += weight * (0.5 * mu_)    * (trCTilde - 3);
-    dilatationEnergy += weight * (0.5 * kappa_) * ((detF - 1) * (detF - 1));
-  }
-
-  return distortionEnergy + dilatationEnergy;
-}
-
 }  // namespace Dune
 
 #endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEOHOOKEENERGY_HH
diff --git a/dune/elasticity/materials/stvenantkirchhoffdensity.hh b/dune/elasticity/materials/stvenantkirchhoffdensity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7d7bfbea9d43af7d3791154ad3dc9750f1fedd12
--- /dev/null
+++ b/dune/elasticity/materials/stvenantkirchhoffdensity.hh
@@ -0,0 +1,68 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFDENSITY_HH
+#define DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::Elasticity {
+
+template<int dim, class field_type = double>
+class StVenantKirchhoffDensity final
+  : public LocalDensity<dim,field_type>
+{
+public:
+
+  /** \brief Constructor with a set of material parameters
+    * \param parameters The material parameters
+    */
+  StVenantKirchhoffDensity(const Dune::ParameterTree& parameters)
+  {
+    // Lame constants
+    mu_ = parameters.template get<double>("mu");
+    lambda_ = parameters.template get<double>("lambda");
+  }
+
+  /** \brief Evaluation with the deformation gradient
+  *
+  * \param x Position of the gradient
+  * \param gradient Deformation gradient
+  */
+  field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
+  {
+    //////////////////////////////////////////////////////////
+    // compute strain E = 1/2 *( F^T F - I)
+    /////////////////////////////////////////////////////////
+
+    Dune::FieldMatrix<field_type,dim,dim> FTF(0);
+    for (int i=0; i<dim; i++)
+      for (int j=0; j<dim; j++)
+        for (int k=0; k<dim; k++)
+          FTF[i][j] += gradient[k][i] * gradient[k][j];
+
+    Dune::FieldMatrix<field_type,dim,dim> E = FTF;
+    for (int i=0; i<dim; i++)
+      E[i][i] -= 1.0;
+    E *= 0.5;
+
+    /////////////////////////////////////////////////////////
+    //  Compute energy
+    /////////////////////////////////////////////////////////
+
+    field_type trE = field_type(0);
+    for (int i=0; i<dim; i++)
+      trE += E[i][i];
+
+    // The trace of E^2, happens to be its squared Frobenius norm
+    field_type trESquared = E.frobenius_norm2();
+
+    return mu_ * trESquared + 0.5 * lambda_ * trE * trE;
+  }
+
+  /** \brief Lame constants */
+  double mu_, lambda_;
+};
+
+} // namespace Dune::Elasticity
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFDENSITY_HH
diff --git a/dune/elasticity/materials/stvenantkirchhoffenergy.hh b/dune/elasticity/materials/stvenantkirchhoffenergy.hh
index 8bb9c67b3347bb8cd1c0447ecf62223ab2827c7b..8b9229c7beb5f3cd7e80fd95768af668bdaddf8b 100644
--- a/dune/elasticity/materials/stvenantkirchhoffenergy.hh
+++ b/dune/elasticity/materials/stvenantkirchhoffenergy.hh
@@ -1,128 +1,26 @@
 #ifndef DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFENERGY_HH
 #define DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFENERGY_HH
 
-#include <dune/common/fmatrix.hh>
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localenergy.hh>
+#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
 
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class StVenantKirchhoffEnergy
-  : public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension> > >
+class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::StVenantKirchhoffDensity")
+StVenantKirchhoffEnergy
+  : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {
-    // grid types
-    typedef typename GridView::Grid::ctype DT;
-    typedef typename GridView::template Codim<0>::Entity Entity;
-
-    // some other sizes
-    enum {gridDim=GridView::dimension};
-    enum {dim=GridView::dimension};
 
+  using Base = Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>;
 public:
 
-    /** \brief Constructor with a set of material parameters
-     * \param parameters The material parameters
-     */
-    StVenantKirchhoffEnergy(const Dune::ParameterTree& parameters)
-    {
-      // Lame constants
-      mu_ = parameters.template get<double>("mu");
-      lambda_ = parameters.template get<double>("lambda");
-    }
-
-    /** \brief Assemble the energy for a single element */
-    field_type energy (const Entity& e,
-               const LocalFiniteElement& localFiniteElement,
-               const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
-
-    /** \brief Lame constants */
-    double mu_, lambda_;
+  // backwards compatibility: forward directly to integral energy
+  StVenantKirchhoffEnergy(const Dune::ParameterTree& parameters)
+  : Base(std::make_shared<Elasticity::StVenantKirchhoffDensity<GridView::dimension,field_type>>(parameters))
+  {}
 };
 
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-StVenantKirchhoffEnergy<GridView,LocalFiniteElement,field_type>::
-energy(const Entity& element,
-       const LocalFiniteElement& localFiniteElement,
-       const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
-{
-    assert(element.type() == localFiniteElement.type());
-    typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
-
-    field_type energy = 0;
-
-    // store gradients of shape functions and base functions
-    std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.localBasis().size());
-    std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.localBasis().size());
-
-    int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                                 : localFiniteElement.localBasis().order() * gridDim;
-
-    const Dune::QuadratureRule<DT, gridDim>& quad
-        = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
-    for (size_t pt=0; pt<quad.size(); pt++) {
-
-        // Local position of the quadrature point
-        const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
-
-        const DT integrationElement = element.geometry().integrationElement(quadPos);
-
-        const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
-
-        DT weight = quad[pt].weight() * integrationElement;
-
-        // get gradients of shape functions
-        localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-
-        // compute gradients of base functions
-        for (size_t i=0; i<gradients.size(); ++i) {
-
-          // transform gradients
-          jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-        }
-
-        Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
-        for (size_t i=0; i<gradients.size(); i++)
-          for (int j=0; j<gridDim; j++)
-            derivative[j].axpy(localConfiguration[i][j], gradients[i]);
-
-        /////////////////////////////////////////////////////////
-        // compute strain E = 1/2 *( F^T F - I)
-        /////////////////////////////////////////////////////////
-
-        Dune::FieldMatrix<field_type,gridDim,gridDim> FTF(0);
-        for (int i=0; i<gridDim; i++)
-          for (int j=0; j<gridDim; j++)
-            for (int k=0; k<gridDim; k++)
-              FTF[i][j] += derivative[k][i] * derivative[k][j];
-
-        Dune::FieldMatrix<field_type,dim,dim> E = FTF;
-        for (int i=0; i<dim; i++)
-          E[i][i] -= 1.0;
-        E *= 0.5;
-
-        /////////////////////////////////////////////////////////
-        //  Compute energy
-        /////////////////////////////////////////////////////////
-
-        field_type trE = field_type(0);
-        for (int i=0; i<dim; i++)
-          trE += E[i][i];
-
-        // The trace of E^2, happens to be its squared Frobenius norm
-        field_type trESquared = E.frobenius_norm2();
-
-        energy += weight * mu_ * trESquared + weight * 0.5 * lambda_ * trE * trE;
-
-    }
-
-    return energy;
-}
-
 }  // namespace Dune
 
 #endif   //#ifndef DUNE_ELASTICITY_MATERIALS_STVENANTKIRCHHOFFENERGY_HH
diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index 28e49701200c4b641e56be92074f5de5b12f3f7b..848c1e5624fdbc5f5a0fa6a4e820bb92efb7541a 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -34,11 +34,12 @@
 #include <dune/elasticity/common/trustregionsolver.hh>
 #include <dune/elasticity/assemblers/localadolcstiffness.hh>
 #include <dune/elasticity/assemblers/feassembler.hh>
-#include <dune/elasticity/materials/stvenantkirchhoffenergy.hh>
-#include <dune/elasticity/materials/henckyenergy.hh>
-#include <dune/elasticity/materials/exphenckyenergy.hh>
-#include <dune/elasticity/materials/mooneyrivlinenergy.hh>
-#include <dune/elasticity/materials/neohookeenergy.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
+#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
+#include <dune/elasticity/materials/henckydensity.hh>
+#include <dune/elasticity/materials/exphenckydensity.hh>
+#include <dune/elasticity/materials/mooneyrivlindensity.hh>
+#include <dune/elasticity/materials/neohookedensity.hh>
 #include <dune/elasticity/materials/neumannenergy.hh>
 #include <dune/elasticity/materials/sumenergy.hh>
 
@@ -270,36 +271,27 @@ int main (int argc, char *argv[]) try
     // Assembler using ADOL-C
     if (mpiHelper.rank()==0)
     std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-    std::shared_ptr<Elasticity::LocalEnergy<GridView,
-                                     FEBasis::LocalView::Tree::FiniteElement,
-                                     std::vector<Dune::FieldVector<adouble, dim> > > > elasticEnergy;
 
-    if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
-      elasticEnergy = std::make_shared<StVenantKirchhoffEnergy<GridView,
-                                                               FEBasis::LocalView::Tree::FiniteElement,
-                                                               adouble> >(materialParameters);
+    std::shared_ptr<Elasticity::LocalDensity<dim,adouble>> elasticDensity;
 
+    if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
+      elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,adouble>>(materialParameters);
     if (parameterSet.get<std::string>("energy") == "neohooke")
-      elasticEnergy = std::make_shared<NeoHookeEnergy<GridView,
-                                                      FEBasis::LocalView::Tree::FiniteElement,
-                                                      adouble> >(materialParameters);
-
+      elasticDensity = std::make_shared<Elasticity::NeoHookeDensity<dim,adouble>>(materialParameters);
     if (parameterSet.get<std::string>("energy") == "hencky")
-      elasticEnergy = std::make_shared<HenckyEnergy<GridView,
-                                                    FEBasis::LocalView::Tree::FiniteElement,
-                                                    adouble> >(materialParameters);
+      elasticDensity = std::make_shared<Elasticity::HenckyDensity<dim,adouble>>(materialParameters);
     if (parameterSet.get<std::string>("energy") == "exphencky")
-      elasticEnergy = std::make_shared<ExpHenckyEnergy<GridView,
-                                                       FEBasis::LocalView::Tree::FiniteElement,
-                                                       adouble> >(materialParameters);
+      elasticDensity = std::make_shared<Elasticity::ExpHenckyDensity<dim,adouble>>(materialParameters);
     if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
-      elasticEnergy = std::make_shared<MooneyRivlinEnergy<GridView,
-                                                               FEBasis::LocalView::Tree::FiniteElement,
-                                                               adouble> >(materialParameters);
+      elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,adouble>>(materialParameters);
 
-    if(!elasticEnergy)
+    if(!elasticDensity)
       DUNE_THROW(Exception, "Error: Selected energy not available!");
 
+    auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+                                          FEBasis::LocalView::Tree::FiniteElement,
+                                          adouble>>(elasticDensity);
+
     auto neumannEnergy = std::make_shared<NeumannEnergy<GridView,
                                                         FEBasis::LocalView::Tree::FiniteElement,
                                                         adouble> >(&neumannBoundary,neumannFunction.get());