diff --git a/src/quasiconvexity-test-micromorphic.cc b/src/quasiconvexity-test-micromorphic.cc
index 10a9b1f38f77528d7f2f352b17813395083a91e2..5f41b951f2d3813b7c5f8fd292fbefbcf86759d0 100644
--- a/src/quasiconvexity-test-micromorphic.cc
+++ b/src/quasiconvexity-test-micromorphic.cc
@@ -76,7 +76,7 @@ void writeDisplacement(const Basis& basis, const Vector& deformation, const Fiel
     // store values of shape functions
     std::vector<Dune::FieldVector<double,dim> > shapeFunctionValues(localFiniteElement.size());
 
-    auto p = element.geometry().center();
+    auto p = referenceElement(element).position(0,0);  // Element center in local coordinates
 
     // get gradients of shape functions
     localFiniteElement.localBasis().evaluateFunction(p, shapeFunctionValues);
diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 94905a9e788e55f873f26a62af04aa268643b0d8..173bf59f7d607773d4f85caca7fd64deb0edb9c8 100644
--- a/src/quasiconvexity-test.cc
+++ b/src/quasiconvexity-test.cc
@@ -139,7 +139,7 @@ void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldM
     std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients(localFiniteElement.size());
     std::vector<Dune::FieldVector<double,dim> > gradients(localFiniteElement.size());
 
-    auto p = element.geometry().center();
+    auto p = referenceElement(element).position(0,0);  // Element center in local coordinates
 
     const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(p);