diff --git a/dune/elasticity/materials/CMakeLists.txt b/dune/elasticity/materials/CMakeLists.txt
index fced7f1f79078f0ed1812d72112c3e7794d7587e..6c150109bb7b407d363df5cafd768d16f7862bd2 100644
--- a/dune/elasticity/materials/CMakeLists.txt
+++ b/dune/elasticity/materials/CMakeLists.txt
@@ -15,6 +15,7 @@ install(FILES
     neohookeenergy.hh
     neumannenergy.hh
     ogdenmaterial.hh
+    quasiconvexitytestdensities.hh
     stvenantkirchhoffdensity.hh
     stvenantkirchhoffenergy.hh
     sumenergy.hh
diff --git a/dune/elasticity/materials/coshenergy.hh b/dune/elasticity/materials/coshenergy.hh
deleted file mode 100644
index f3d1070b00bccc74660531888830bb82e7c58ed1..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/coshenergy.hh
+++ /dev/null
@@ -1,105 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_COSHENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_COSHENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class CosHEnergy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  CosHEnergy(const Dune::ParameterTree& parameters)
-  {
-    eps_ = parameters.get<double>("eps");
-    K0_ = parameters.get<double>("K0");
-  }
-
-  /** \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;
-
-  double eps_;
-  double K0_;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-CosHEnergy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto K = normSquared / (2*det);
-
-    using std::cosh;
-    using std::pow;
-
-    auto h2 = cosh(pow( (K-K0_)/eps_, 4)) - 1;
-
-    auto energyDensity = h2 + 0;//100*f2;
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_COSHENERGY_HH
diff --git a/dune/elasticity/materials/dacorognaenergy.hh b/dune/elasticity/materials/dacorognaenergy.hh
deleted file mode 100644
index a087471d4e76adb992d04d6973609e552203ba2c..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/dacorognaenergy.hh
+++ /dev/null
@@ -1,97 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune::Elasticity {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class DacorognaEnergy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  DacorognaEnergy(const Dune::ParameterTree& parameters)
-  {
-    gamma_ = parameters.template get<double>("gamma");
-  }
-
-  /** \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;
-
-  double gamma_;
-
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-DacorognaEnergy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto energyDensity = normSquared * ( normSquared - gamma_ * det);
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune::Elasticity
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_HH
diff --git a/dune/elasticity/materials/expacoshenergy.hh b/dune/elasticity/materials/expacoshenergy.hh
deleted file mode 100644
index b4baf69ab3f2540388ea60b2db5d1be05ddfdc58..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/expacoshenergy.hh
+++ /dev/null
@@ -1,106 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class ExpACosHEnergy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  ExpACosHEnergy(const Dune::ParameterTree& parameters)
-  {
-    c1_ = parameters.get<double>("c1");
-    c2_ = parameters.get<double>("c2");
-  }
-
-  /** \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;
-
-  double c1_;
-  double c2_;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-ExpACosHEnergy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto K = normSquared / (2*det);
-
-    using std::acosh;
-
-    auto h = exp(c1_ * acosh(K) * acosh(K));
-
-    auto f = - c2_ * log(det);
-
-    auto energyDensity = h + f;
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
diff --git a/dune/elasticity/materials/fungenergy.hh b/dune/elasticity/materials/fungenergy.hh
deleted file mode 100644
index 0d79c60b742cb495307151cfe8e49924d4d5f8f7..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/fungenergy.hh
+++ /dev/null
@@ -1,121 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_FUNGENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_FUNGENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune::Elasticity {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class FungEnergy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  FungEnergy(const Dune::ParameterTree& parameters)
-  {
-    c_ = parameters.template get<double>("c");
-    A_[0] = parameters.template get<double>("A1");
-    A_[1] = parameters.template get<double>("A2");
-    A_[2] = parameters.template get<double>("A3");
-    A_[3] = parameters.template get<double>("A4");
-    A_[4] = parameters.template get<double>("A5");
-    A_[5] = parameters.template get<double>("A6");
-  }
-
-  /** \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;
-
-  double c_;
-
-  std::array<double, 6> A_;
-
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-FungEnergy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    FieldMatrix<field_type,gridDim,gridDim> greenLagrangeStrain(0);
-    for (int i=0; i<gridDim; i++)
-      for (int j=0; j<gridDim; j++)
-      {
-        greenLagrangeStrain[i][j] = -(i==j);
-        for (int k=0; k<gridDim; k++)
-          greenLagrangeStrain[i][j] += derivative[k][i] * derivative[k][j];
-      }
-
-    //auto psi = A1 * E11*E11 + A2 * E22 * E22 + 2*A3*E11*E22 + A4 * E12 * E12 + 2*A5 * E12 * E11 + 2*A6 * E12 * E22;
-    auto  psi = A_[0] * greenLagrangeStrain[0][0] * greenLagrangeStrain[0][0]
-              + A_[1] * greenLagrangeStrain[1][1] * greenLagrangeStrain[1][1]
-              + 2*A_[2]*greenLagrangeStrain[0][0] * greenLagrangeStrain[1][1]
-              + A_[3] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][1]
-              + 2*A_[4] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][0]
-              + 2*A_[5] * greenLagrangeStrain[0][1] * greenLagrangeStrain[1][1];
-
-    using std::exp;
-
-    auto energyDensity = 0.5 * c_ * (exp(psi) - 1);
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune::Elasticity
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_FUNGENERGY_HH
diff --git a/dune/elasticity/materials/neffenergy.hh b/dune/elasticity/materials/neffenergy.hh
deleted file mode 100644
index f1de30e75a35d28019feaa8addec3e2e17a0597e..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/neffenergy.hh
+++ /dev/null
@@ -1,96 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_NEFFENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_NEFFENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffEnergy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  NeffEnergy(const Dune::ParameterTree& parameters)
-  {
-    scaling_ = parameters.get<double>("scaling");
-  }
-
-  /** \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;
-
-private:
-  double scaling_;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeffEnergy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    energy += qp.weight() * integrationElement *
-           ((scaling_/(4*std::sqrt(3))) * (normSquared*normSquared/(det*det) -4) + std::log(det) - det);
-  }
-
-  return energy;
-}
-
-}  // namespace Dune
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFENERGY_HH
diff --git a/dune/elasticity/materials/nefflogenergy.hh b/dune/elasticity/materials/nefflogenergy.hh
deleted file mode 100644
index 2e9acc0633e15d8b01aea37aac8e5b088f638197..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/nefflogenergy.hh
+++ /dev/null
@@ -1,94 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffLogEnergy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  NeffLogEnergy(const Dune::ParameterTree& parameters)
-  {}
-
-  /** \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;
-
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeffLogEnergy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto energyDensity = std::log( normSquared / (2*det) + std::sqrt(normSquared*normSquared/(4*det*det) -1) );
-    energyDensity = energyDensity * energyDensity;
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY_HH
diff --git a/dune/elasticity/materials/nefflogenergy2.hh b/dune/elasticity/materials/nefflogenergy2.hh
deleted file mode 100644
index e7449b8ec2073d07e28fbe46bd82ecc73b1b7a75..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/nefflogenergy2.hh
+++ /dev/null
@@ -1,104 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY2_HH
-#define DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY2_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffLogEnergy2 final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  NeffLogEnergy2(const Dune::ParameterTree& parameters,
-                 double w0=1, double w1=1)
-  {
-    c0_ = parameters.template get<double>("c0");
-    weights_ = {w0,w1};
-  }
-
-  /** \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;
-
-  double c0_;
-
-  // This is needed to compute the two parts of the energy separately
-  std::array<double,2> weights_;
-
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeffLogEnergy2<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto K = 0.5 * normSquared / det;
-
-    auto energyDensity = weights_[0] * c0_ * K * std::log(K) + weights_[1] * std::log(det);
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY2_HH
diff --git a/dune/elasticity/materials/neffreducedenergy.hh b/dune/elasticity/materials/neffreducedenergy.hh
deleted file mode 100644
index 85cfd0afc4460b0ce3196127b2b397ab39f05fe1..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/neffreducedenergy.hh
+++ /dev/null
@@ -1,93 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_NEFFREDUCEDENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_NEFFREDUCEDENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffReducedEnergy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  NeffReducedEnergy(const Dune::ParameterTree& parameters)
-  {}
-
-  /** \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;
-
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeffReducedEnergy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto energyDensity = normSquared / det;
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFREDUCEDENERGY_HH
diff --git a/dune/elasticity/materials/neffw2energy.hh b/dune/elasticity/materials/neffw2energy.hh
deleted file mode 100644
index 52525775bee9231f5ab0c3a81d93cba58564151a..0000000000000000000000000000000000000000
--- a/dune/elasticity/materials/neffw2energy.hh
+++ /dev/null
@@ -1,102 +0,0 @@
-#ifndef DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune::Elasticity {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffW2Energy final
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // 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:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  NeffW2Energy(const Dune::ParameterTree& parameters)
-  {
-    c_ = 1.0;
-    if (parameters.hasKey("c"))
-      c_ = parameters.template get<double>("c");
-  }
-
-  /** \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;
-  double c_;
-
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeffW2Energy<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());
-
-  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 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());
-
-    // 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]);
-
-    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]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto K = normSquared / (2*det);
-
-    using std::acosh;
-
-    auto energyDensity = K + acosh(K) - c_*log(det);
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune::Elasticity
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_HH
diff --git a/dune/elasticity/materials/quasiconvexitytestdensities.hh b/dune/elasticity/materials/quasiconvexitytestdensities.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7839e47f8682b48e809a38bb0b0f19312b15bb25
--- /dev/null
+++ b/dune/elasticity/materials/quasiconvexitytestdensities.hh
@@ -0,0 +1,402 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_ELASTICITY_MATERIAL_QUASICONVEXITYTESTDENSITIES_HH
+#define DUNE_ELASTICITY_MATERIAL_QUASICONVEXITYTESTDENSITIES_HH
+
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune:: Elasticity
+{
+
+template<int dim, class field_type, class ctype=double>
+struct BurkholderDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  BurkholderDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    p_ = parameters.template get<double>("p");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto lambdaMaxSquared = 0.5 * (normSquared + std::sqrt(normSquared*normSquared - 4*det*det));
+    auto lambdaMax = std::sqrt(lambdaMaxSquared);
+
+    return -0.5 * p_ * det * std::pow(lambdaMax,p_-2) - 0.5*(2-p_) * std::pow(lambdaMax,p_);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double p_;
+};
+
+template<int dim, class field_type, class ctype=double>
+struct CosHDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  CosHDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    eps_ = parameters.get<double>("eps");
+    K0_ = parameters.get<double>("K0");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*det);
+
+    using std::cosh;
+    using std::pow;
+
+    auto h2 = cosh(pow( (K-K0_)/eps_, 4)) - 1;
+
+    return h2 + 0;//100*f2;
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double eps_;
+  double K0_;
+};
+
+template<int dim, class field_type, class ctype=double>
+struct DacorognaDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  DacorognaDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    gamma_ = parameters.template get<double>("gamma");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    return normSquared * ( normSquared - gamma_ * det);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double gamma_;
+};
+
+template<int dim, class field_type, class ctype=double>
+struct ExpACosHDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  ExpACosHDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    c1_ = parameters.get<double>("c1");
+    c2_ = parameters.get<double>("c2");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*det);
+
+    using std::acosh;
+    return exp(c1_ * acosh(K) * acosh(K)) - c2_ * log(det);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double c1_;
+  double c2_;
+};
+
+template<int dim, class field_type, class ctype=double>
+struct FungDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  FungDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    c_ = parameters.template get<double>("c");
+    A_[0] = parameters.template get<double>("A1");
+    A_[1] = parameters.template get<double>("A2");
+    A_[2] = parameters.template get<double>("A3");
+    A_[3] = parameters.template get<double>("A4");
+    A_[4] = parameters.template get<double>("A5");
+    A_[5] = parameters.template get<double>("A6");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+
+    FieldMatrix<field_type,dim,dim> greenLagrangeStrain(0);
+    for (int i=0; i<dim; i++)
+      for (int j=0; j<dim; j++)
+      {
+        greenLagrangeStrain[i][j] = -(i==j);
+        for (int k=0; k<dim; k++)
+          greenLagrangeStrain[i][j] += F[k][i] * F[k][j];
+      }
+
+    //auto psi = A1 * E11*E11 + A2 * E22 * E22 + 2*A3*E11*E22 + A4 * E12 * E12 + 2*A5 * E12 * E11 + 2*A6 * E12 * E22;
+    auto  psi = A_[0] * greenLagrangeStrain[0][0] * greenLagrangeStrain[0][0]
+              + A_[1] * greenLagrangeStrain[1][1] * greenLagrangeStrain[1][1]
+              + 2*A_[2]*greenLagrangeStrain[0][0] * greenLagrangeStrain[1][1]
+              + A_[3] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][1]
+              + 2*A_[4] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][0]
+              + 2*A_[5] * greenLagrangeStrain[0][1] * greenLagrangeStrain[1][1];
+
+    using std::exp;
+
+    return 0.5 * c_ * (exp(psi) - 1);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double c_;
+  std::array<double, 6> A_;
+};
+
+template<int dim, class field_type, class ctype=double>
+struct KPowerDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  KPowerDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    p_ = parameters.template get<double>("p");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto detF = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*detF);
+
+    return std::pow(K, p_);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double p_;
+};
+
+
+template<int dim, class field_type, class ctype=double>
+struct MagicDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  MagicDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    c_ = 1.0;
+    if (parameters.hasKey("c"))
+      c_ = parameters.template get<double>("c");
+  }
+
+  /** \brief Implement space-dependent coefficients */
+  double c(const Dune::FieldVector<ctype,dim>& x) const
+  {
+    return c_;
+#if 0  // Three circles
+    bool inCircle0 = (x-FieldVector<double,dim>{-0.5,0}).two_norm() <0.2;
+    bool inCircle1 = (x-FieldVector<double,dim>{0.3535,0.3535}).two_norm() <0.2;
+    bool inCircle2 = (x-FieldVector<double,dim>{0.3535,-0.3535}).two_norm() <0.2;
+
+    return (inCircle0 or inCircle1 or inCircle2) ? c_ : 1;
+#endif
+//    return (x.two_norm() < 0.9) ? 1 : c_;
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*det);
+
+    using std::acosh;
+    using std::sqrt;
+
+    return K + sqrt(K*K-1) - acosh(K) + c(x)*log(det);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double c_;
+};
+
+
+template<int dim, class field_type, class ctype=double>
+struct NeffDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  NeffDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    scaling_ = parameters.get<double>("scaling");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    return (scaling_/(4*std::sqrt(3))) * (normSquared*normSquared/(det*det) -4) + std::log(det) - det;
+  }
+
+private:
+  FieldMatrix<ctype,dim,dim> F0_;
+  double scaling_;
+};
+
+
+template<int dim, class field_type, class ctype=double>
+struct NeffLogDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  NeffLogDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {}
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto energyDensity = std::log( normSquared / (2*det) + std::sqrt(normSquared*normSquared/(4*det*det) -1) );
+    return energyDensity * energyDensity;
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+};
+
+
+template<int dim, class field_type, class ctype=double>
+struct NeffLog2Density final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  NeffLog2Density(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters,
+                  double w0=1, double w1=1)
+  : F0_(F0)
+  {
+    c0_ = parameters.template get<double>("c0");
+    weights_ = {w0,w1};
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = 0.5 * normSquared / det;
+
+    return weights_[0] * c0_ * K * std::log(K) + weights_[1] * std::log(det);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double c0_;
+
+  // This is needed to compute the two parts of the energy separately
+  std::array<double,2> weights_;
+};
+
+
+template<int dim, class field_type, class ctype=double>
+struct NeffReducedDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  NeffReducedDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {}
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    return normSquared / det;
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+};
+
+
+template<int dim, class field_type, class ctype=double>
+struct NeffW2Density final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  NeffW2Density(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    c_ = 1.0;
+    if (parameters.hasKey("c"))
+      c_ = parameters.template get<double>("c");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*det);
+
+    using std::log;
+    using std::acosh;
+
+    return K + acosh(K) - c_*log(det);
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double c_;
+};
+
+
+template<int dim, class field_type, class ctype=double>
+struct VossDensity final
+: public Elasticity::LocalDensity<dim,field_type,ctype>
+{
+  VossDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
+  {
+    c_ = parameters.template get<double>("c");
+  }
+
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
+  {
+    auto F = derivative + F0_;
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto KBold = normSquared / (2*det);
+
+    using std::acosh;
+
+    auto K = KBold + std::sqrt(KBold*KBold - 1);
+
+    auto h = (1.0/(2*K)) * (1 + K*K + (1+K)*std::sqrt(1 + 6*K + K*K)) - std::log(K)
+           + std::log(1 + 4*K + K*K + (1+K)*std::sqrt(1 + 6*K + K*K));
+
+    return h + c_ * std::log( det );
+  }
+
+  FieldMatrix<ctype,dim,dim> F0_;
+  double c_;
+};
+
+}   // namespace Dune::Elasticity
+
+#endif
diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 0ca4a1a42630f814f28aa660e11c9764356fca28..42df972b910993afc1e48ed5f41637c00d78244f 100644
--- a/src/quasiconvexity-test.cc
+++ b/src/quasiconvexity-test.cc
@@ -37,16 +37,8 @@
 #include <dune/elasticity/common/trustregionsolver.hh>
 #include <dune/elasticity/assemblers/localadolcstiffness.hh>
 #include <dune/elasticity/assemblers/feassembler.hh>
-//#include <dune/elasticity/materials/neffenergy.hh>
-//#include <dune/elasticity/materials/coshenergy.hh>
-//#include <dune/elasticity/materials/dacorognaenergy.hh>
-//#include <dune/elasticity/materials/expacoshenergy.hh>
-//#include <dune/elasticity/materials/fungenergy.hh>
 #include <dune/elasticity/materials/localintegralenergy.hh>
-//#include <dune/elasticity/materials/nefflogenergy.hh>
-//#include <dune/elasticity/materials/nefflogenergy2.hh>
-//#include <dune/elasticity/materials/neffreducedenergy.hh>
-//#include <dune/elasticity/materials/neffw2energy.hh>
+#include <dune/elasticity/materials/quasiconvexitytestdensities.hh>
 
 // grid dimension
 const int dim = 2;
@@ -55,164 +47,6 @@ const int order = 1;
 
 using namespace Dune;
 
-template<int dim, class field_type, class ctype=double>
-struct BurkholderDensity final
-: public Elasticity::LocalDensity<dim,field_type,ctype>
-{
-  BurkholderDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
-  : F0_(F0)
-  {
-    p_ = parameters.template get<double>("p");
-  }
-
-  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
-  {
-    auto F = derivative + F0_;
-    auto det = F.determinant();
-    auto normSquared = F.frobenius_norm2();
-
-    auto lambdaMaxSquared = 0.5 * (normSquared + std::sqrt(normSquared*normSquared - 4*det*det));
-    auto lambdaMax = std::sqrt(lambdaMaxSquared);
-
-    return -0.5 * p_ * det * std::pow(lambdaMax,p_-2) - 0.5*(2-p_) * std::pow(lambdaMax,p_);
-  }
-
-  FieldMatrix<ctype,dim,dim> F0_;
-  double p_;
-};
-
-
-template<int dim, class field_type, class ctype=double>
-struct ExpACosHDensity final
-: public Elasticity::LocalDensity<dim,field_type,ctype>
-{
-  ExpACosHDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
-  : F0_(F0)
-  {
-    c1_ = parameters.get<double>("c1");
-    c2_ = parameters.get<double>("c2");
-  }
-
-  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
-  {
-    auto F = derivative + F0_;
-    auto det = F.determinant();
-    auto normSquared = F.frobenius_norm2();
-
-    auto K = normSquared / (2*det);
-
-    using std::acosh;
-    return exp(c1_ * acosh(K) * acosh(K)) - c2_ * log(det);
-  }
-
-  FieldMatrix<ctype,dim,dim> F0_;
-  double c1_;
-  double c2_;
-};
-
-
-template<int dim, class field_type, class ctype=double>
-struct KPowerDensity final
-: public Elasticity::LocalDensity<dim,field_type,ctype>
-{
-  KPowerDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
-  : F0_(F0)
-  {
-    p_ = parameters.template get<double>("p");
-  }
-
-  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
-  {
-    auto F = derivative + F0_;
-    auto detF = F.determinant();
-    auto normSquared = F.frobenius_norm2();
-
-    auto K = normSquared / (2*detF);
-
-    return std::pow(K, p_);
-  }
-
-  FieldMatrix<ctype,dim,dim> F0_;
-  double p_;
-};
-
-
-template<int dim, class field_type, class ctype=double>
-struct MagicDensity final
-: public Elasticity::LocalDensity<dim,field_type,ctype>
-{
-  MagicDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
-  : F0_(F0)
-  {
-    c_ = 1.0;
-    if (parameters.hasKey("c"))
-      c_ = parameters.template get<double>("c");
-  }
-
-  /** \brief Implement space-dependent coefficients */
-  double c(const Dune::FieldVector<ctype,dim>& x) const
-  {
-    return c_;
-#if 0  // Three circles
-    bool inCircle0 = (x-FieldVector<double,dim>{-0.5,0}).two_norm() <0.2;
-    bool inCircle1 = (x-FieldVector<double,dim>{0.3535,0.3535}).two_norm() <0.2;
-    bool inCircle2 = (x-FieldVector<double,dim>{0.3535,-0.3535}).two_norm() <0.2;
-
-    return (inCircle0 or inCircle1 or inCircle2) ? c_ : 1;
-#endif
-//    return (x.two_norm() < 0.9) ? 1 : c_;
-  }
-
-  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
-  {
-    auto F = derivative + F0_;
-    auto det = F.determinant();
-    auto normSquared = F.frobenius_norm2();
-
-    auto K = normSquared / (2*det);
-
-    using std::acosh;
-
-    return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det);
-  }
-
-  FieldMatrix<ctype,dim,dim> F0_;
-  double c_;
-};
-
-
-template<int dim, class field_type, class ctype=double>
-struct VossDensity final
-: public Elasticity::LocalDensity<dim,field_type,ctype>
-{
-  VossDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
-  : F0_(F0)
-  {
-    c_ = parameters.template get<double>("c");
-  }
-
-  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
-  {
-    auto F = derivative + F0_;
-    auto det = F.determinant();
-    auto normSquared = F.frobenius_norm2();
-
-    auto KBold = normSquared / (2*det);
-
-    using std::acosh;
-
-    auto K = KBold + std::sqrt(KBold*KBold - 1);
-
-    auto h = (1.0/(2*K)) * (1 + K*K + (1+K)*std::sqrt(1 + 6*K + K*K)) - std::log(K)
-           + std::log(1 + 4*K + K*K + (1+K)*std::sqrt(1 + 6*K + K*K));
-
-    return h + c_ * std::log( det );
-  }
-
-  FieldMatrix<ctype,dim,dim> F0_;
-  double c_;
-};
-
 template <class Basis, class Vector>
 void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldMatrix<double,2,2>& F0, std::string filename)
 {
@@ -625,92 +459,53 @@ int main (int argc, char *argv[]) try
 
     // Assembler using ADOL-C
     std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-    std::shared_ptr<Elasticity::LocalEnergy<FEBasis::LocalView,adouble> > elasticEnergy;
+    std::shared_ptr<Elasticity::LocalDensity<dim,adouble> > density;
 
-    if (parameterSet.get<std::string>("energy") == "kpower")
-    {
-      auto density = std::make_shared<KPowerDensity<dim,adouble> >(F0, materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
-                                                                       adouble> >(density);
-    }
+    if (parameterSet.get<std::string>("energy") == "burkholder")
+      density = std::make_shared<Elasticity::BurkholderDensity<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "cosh")
+      density = std::make_shared<Elasticity::CosHDensity<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "dacorogna")
+      density = std::make_shared<Elasticity::DacorognaDensity<dim,adouble> >(F0,materialParameters);
 
-//     if (parameterSet.get<std::string>("energy") == "neff")
-//       elasticEnergy = std::make_shared<NeffEnergy<GridView,
-//                                                       FEBasis::LocalView::Tree::FiniteElement,
-//                                                       adouble> >(materialParameters);
-//
-//     if (parameterSet.get<std::string>("energy") == "nefflog")
-//       elasticEnergy = std::make_shared<NeffLogEnergy<GridView,
-//                                                       FEBasis::LocalView::Tree::FiniteElement,
-//                                                       adouble> >(materialParameters);
-//
-//     if (parameterSet.get<std::string>("energy") == "nefflog2")
-//       elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
-//                                                       FEBasis::LocalView::Tree::FiniteElement,
-//                                                       adouble> >(materialParameters);
-//
-//     if (parameterSet.get<std::string>("energy") == "neffreduced")
-//       elasticEnergy = std::make_shared<NeffReducedEnergy<GridView,
-//                                                          FEBasis::LocalView::Tree::FiniteElement,
-//                                                          adouble> >(materialParameters);
-//
-//     if (parameterSet.get<std::string>("energy") == "cosh")
-//       elasticEnergy = std::make_shared<CosHEnergy<GridView,
-//                                                          FEBasis::LocalView::Tree::FiniteElement,
-//                                                          adouble> >(materialParameters);
-//
-#if 0
-    if (parameterSet.get<std::string>("energy") == "expacosh")
-      elasticEnergy = std::make_shared<ExpACosHEnergy<GridView,
-                                                         FEBasis::LocalView::Tree::FiniteElement,
-                                                         adouble> >(materialParameters);
-#else
     if (parameterSet.get<std::string>("energy") == "expacosh")
-    {
-      auto density = std::make_shared<ExpACosHDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
-                                                                       adouble> >(density);
-    }
-#endif
-//     if (parameterSet.get<std::string>("energy") == "fung")
-//       elasticEnergy = std::make_shared<Elasticity::FungEnergy<GridView,
-//                                                   FEBasis::LocalView::Tree::FiniteElement,
-//                                                   adouble> >(materialParameters);
+      density = std::make_shared<Elasticity::ExpACosHDensity<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "fung")
+      density = std::make_shared<Elasticity::FungDensity<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "kpower")
+      density = std::make_shared<Elasticity::KPowerDensity<dim,adouble> >(F0, materialParameters);
 
     if (parameterSet.get<std::string>("energy") == "magic")
-    {
-      auto density = std::make_shared<MagicDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
-                                                                       adouble> >(density);
-    }
+      density = std::make_shared<Elasticity::MagicDensity<dim,adouble> >(F0,materialParameters);
 
-    if (parameterSet.get<std::string>("energy") == "burkholder")
-    {
-      auto density = std::make_shared<BurkholderDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
-                                                                       adouble> >(density);
-    }
+    if (parameterSet.get<std::string>("energy") == "neff")
+      density = std::make_shared<Elasticity::NeffDensity<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "nefflog")
+      density = std::make_shared<Elasticity::NeffLogDensity<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "nefflog2")
+      density = std::make_shared<Elasticity::NeffLog2Density<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "neffreduced")
+      density = std::make_shared<Elasticity::NeffReducedDensity<dim,adouble> >(F0,materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "w2")
+      density = std::make_shared<Elasticity::NeffW2Density<dim,adouble> >(F0,materialParameters);
 
-//     if (parameterSet.get<std::string>("energy") == "w2")
-//       elasticEnergy = std::make_shared<Elasticity::NeffW2Energy<GridView,
-//                                                                 FEBasis::LocalView::Tree::FiniteElement,
-//                                                                 adouble> >(materialParameters);
-//
-//     if (parameterSet.get<std::string>("energy") == "dacorogna")
-//       elasticEnergy = std::make_shared<Elasticity::DacorognaEnergy<GridView,
-//                                                                 FEBasis::LocalView::Tree::FiniteElement,
-//                                                                 adouble> >(materialParameters);
-//
     if (parameterSet.get<std::string>("energy") == "voss")
-    {
-      auto density = std::make_shared<VossDensity<dim,adouble> >(F0,materialParameters);
-      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
-                                                                       adouble> >(density);
-    }
+      auto density = std::make_shared<Elasticity::VossDensity<dim,adouble> >(F0,materialParameters);
 
-    if(!elasticEnergy)
+    if(!density)
       DUNE_THROW(Exception, "Error: Selected energy not available!");
 
+    auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
+                                                                          adouble> >(density);
+
     auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy);
 
     Elasticity::FEAssembler<FEBasis,SolutionType> assembler(feBasis, localADOLCStiffness);