diff --git a/dune/elasticity/materials/exphenckydensity.hh b/dune/elasticity/materials/exphenckydensity.hh
index b20a440fb061768f3c2cc6216dddeb23b0eedb6a..9258e57a9a507283f51809fa0d0db67ba904d7eb 100644
--- a/dune/elasticity/materials/exphenckydensity.hh
+++ b/dune/elasticity/materials/exphenckydensity.hh
@@ -8,7 +8,7 @@
 
 namespace Dune::Elasticity {
 
-template<int dim, class field_type = double>
+template<int dim, class field_type = double, class ctype = double>
 class ExpHenckyDensity final
   : public LocalDensity<dim,field_type>
 {
@@ -40,7 +40,7 @@ public:
   * \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
+  field_type operator() (const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
   {
     /////////////////////////////////////////////////////////
     // compute C = F^TF
diff --git a/dune/elasticity/materials/henckydensity.hh b/dune/elasticity/materials/henckydensity.hh
index 81a69562b711dd280c96a5f6e88ee0f46d5d1acd..a4dc90edd563568309873a4e0bb58407f836ab4a 100644
--- a/dune/elasticity/materials/henckydensity.hh
+++ b/dune/elasticity/materials/henckydensity.hh
@@ -8,7 +8,7 @@
 
 namespace Dune::Elasticity {
 
-template<int dim, class field_type = double>
+template<int dim, class field_type = double, class ctype = double>
 class HenckyDensity final
   : public LocalDensity<dim,field_type>
 {
@@ -30,7 +30,7 @@ public:
   * \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
+  field_type operator() (const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
   {
     /////////////////////////////////////////////////////////
     // compute C = F^TF
diff --git a/dune/elasticity/materials/localdensity.hh b/dune/elasticity/materials/localdensity.hh
index c091b4d7f84961c6e93c4c6870ae7c8d1a328f7e..80156d942ad7d1247bf1fd874d17afc7d001e9c1 100644
--- a/dune/elasticity/materials/localdensity.hh
+++ b/dune/elasticity/materials/localdensity.hh
@@ -5,8 +5,12 @@
 
 namespace Dune::Elasticity {
 
-/** \brief A base class for energy densities to be evaluated in an integral energy */
-template<int dim, class field_type = double>
+/** \brief A base class for energy densities to be evaluated in an integral energy
+ *
+ * \tparam field_type type of the gradient entries
+ * \tparam ctype type of the coordinates, i.e., the x entries
+ */
+template<int dim, class field_type = double, class ctype = double>
 class LocalDensity
 {
 
@@ -17,7 +21,7 @@ public:
    * \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;
+  virtual field_type operator() (const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const = 0;
 
 };
 
diff --git a/dune/elasticity/materials/localintegralenergy.hh b/dune/elasticity/materials/localintegralenergy.hh
index 515d4322a4bb1fb3822ba060b1cfd5ee0816a25e..1ba6759d57845360700910e271e38a984bbbc32a 100644
--- a/dune/elasticity/materials/localintegralenergy.hh
+++ b/dune/elasticity/materials/localintegralenergy.hh
@@ -27,7 +27,7 @@ public:
 
   /** \brief Constructor with a local energy density
     */
-  LocalIntegralEnergy(const std::shared_ptr<LocalDensity<dim,field_type>>& ld)
+  LocalIntegralEnergy(const std::shared_ptr<LocalDensity<dim,field_type,DT>>& ld)
   : localDensity_(ld)
   {}
 
@@ -41,7 +41,7 @@ public:
                      const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
 
 protected:
-  const std::shared_ptr<LocalDensity<dim,field_type>> localDensity_ = nullptr;
+  const std::shared_ptr<LocalDensity<dim,field_type,DT>> localDensity_ = nullptr;
 
 };
 
diff --git a/dune/elasticity/materials/mooneyrivlindensity.hh b/dune/elasticity/materials/mooneyrivlindensity.hh
index 4dcacdcde00188d0af641449690e09371149f99d..8069301faf3fbf547f8457b61f2be02d01aabdff 100644
--- a/dune/elasticity/materials/mooneyrivlindensity.hh
+++ b/dune/elasticity/materials/mooneyrivlindensity.hh
@@ -8,7 +8,7 @@
 
 namespace Dune::Elasticity {
 
-template<int dim, class field_type = double>
+template<int dim, class field_type = double, class ctype = double>
 class MooneyRivlinDensity final
     : public LocalDensity<dim,field_type>
 {
@@ -43,7 +43,7 @@ public:
   * \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
+  field_type operator() (const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
   {
     /////////////////////////////////////////////////////////
     // compute C = F^TF
diff --git a/dune/elasticity/materials/neohookedensity.hh b/dune/elasticity/materials/neohookedensity.hh
index 47aa6d61ab2f71f6ba7df83e4ec9f64264522baa..1541e3c25554e72542b62a2e1760df614ee070c2 100644
--- a/dune/elasticity/materials/neohookedensity.hh
+++ b/dune/elasticity/materials/neohookedensity.hh
@@ -8,7 +8,7 @@
 
 namespace Dune::Elasticity {
 
-template<int dim, class field_type = double>
+template<int dim, class field_type = double, class ctype = double>
 class NeoHookeDensity final
   : public LocalDensity<dim,field_type>
 {
@@ -30,7 +30,7 @@ public:
   * \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
+  field_type operator() (const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
   {
     /////////////////////////////////////////////////////////
     // compute C = F^TF
diff --git a/dune/elasticity/materials/stvenantkirchhoffdensity.hh b/dune/elasticity/materials/stvenantkirchhoffdensity.hh
index 7d7bfbea9d43af7d3791154ad3dc9750f1fedd12..31a535aa178e28b6820cef9621656367f2e872df 100644
--- a/dune/elasticity/materials/stvenantkirchhoffdensity.hh
+++ b/dune/elasticity/materials/stvenantkirchhoffdensity.hh
@@ -7,7 +7,7 @@
 
 namespace Dune::Elasticity {
 
-template<int dim, class field_type = double>
+template<int dim, class field_type = double, class ctype = double>
 class StVenantKirchhoffDensity final
   : public LocalDensity<dim,field_type>
 {
@@ -28,7 +28,7 @@ public:
   * \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
+  field_type operator() (const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
   {
     //////////////////////////////////////////////////////////
     // compute strain E = 1/2 *( F^T F - I)