diff --git a/dune/elasticity/materials/neffmagicenergy.hh b/dune/elasticity/materials/neffmagicenergy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ec6f4e0a2af726e5c8440cdffa0c027bf5c96807
--- /dev/null
+++ b/dune/elasticity/materials/neffmagicenergy.hh
@@ -0,0 +1,96 @@
+#ifndef DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_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 NeffMagicEnergy
+  : 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
+   */
+  NeffMagicEnergy(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
+NeffMagicEnergy<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 + std::sqrt(K*K-1) - acosh(K) + log(det);
+
+    energy += qp.weight() * integrationElement * energyDensity;
+  }
+
+  return energy;
+}
+
+}  // namespace Dune::Elasticity
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 36c82af9e3836910ca179116a1a991c48722a55b..91459e70a18b196353b18419bf8a7b376a866e24 100644
--- a/src/quasiconvexity-test.cc
+++ b/src/quasiconvexity-test.cc
@@ -39,6 +39,7 @@
 #include <dune/elasticity/materials/expacoshenergy.hh>
 #include <dune/elasticity/materials/nefflogenergy.hh>
 #include <dune/elasticity/materials/nefflogenergy2.hh>
+#include <dune/elasticity/materials/neffmagicenergy.hh>
 #include <dune/elasticity/materials/neffreducedenergy.hh>
 
 // grid dimension
@@ -271,6 +272,11 @@ int main (int argc, char *argv[]) try
                                                          FEBasis::LocalView::Tree::FiniteElement,
                                                          adouble> >(materialParameters);
 
+    if (parameterSet.get<std::string>("energy") == "magic")
+      elasticEnergy = std::make_shared<Elasticity::NeffMagicEnergy<GridView,
+                                                                   FEBasis::LocalView::Tree::FiniteElement,
+                                                                   adouble> >(materialParameters);
+
     if(!elasticEnergy)
       DUNE_THROW(Exception, "Error: Selected energy not available!");