diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 58fb6b3608894e0072b0ef14e86366d027c2f703..3c78ef0a5935d04bba949f10aef72af67c14c501 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -294,17 +294,14 @@ public:
     {
       QPcounter++;
       const auto& quadPos = quadPoint.position();
-      const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+      const auto geometryJacobianInverse = geometry.jacobianInverse(quadPos);
       const auto integrationElement = geometry.integrationElement(quadPos);
 
-      std::vector< Dune::FieldMatrix< double, 1, dim> > referenceGradients;
-      localFiniteElement.localBasis().evaluateJacobian(quadPos,
-                                                      referenceGradients);
-      // Compute the shape function gradients on the grid element
-      std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
+      std::vector<Dune::FieldMatrix<double,1,dim> > jacobians;
+      localFiniteElement.localBasis().evaluateJacobian(quadPos, jacobians);
 
-      for (size_t i=0; i<gradients.size(); i++)
-        jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+      for (size_t i=0; i< jacobians.size(); i++)
+        jacobians[i] = jacobians[i] * geometryJacobianInverse;
 
       // Compute the material phase that the current quadrature point is in
       const auto phase = localIndicatorFunction(quadPos);
@@ -315,10 +312,10 @@ public:
           size_t row = localView.tree().child(l).localIndex(j);
           // (scaled)  Deformation gradient of the test basis function
           MatrixRT defGradientV(0);
-          defGradientV[l][0] = gradients[j][0];                       // Y
-          defGradientV[l][1] = gradients[j][1];                       //X2
-  //         defGradientV[l][2] = (1.0/gamma)*gradients[j][2];           //X3
-          defGradientV[l][2] = gradients[j][2];           //X3
+          defGradientV[l][0] = jacobians[j][0][0];                       // Y
+          defGradientV[l][1] = jacobians[j][0][1];                       //X2
+  //         defGradientV[l][2] = (1.0/gamma)*jacobians[j][0][2];           //X3
+          defGradientV[l][2] = jacobians[j][0][2];           //X3
           
           defGradientV = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
           // "phi*phi"-part
@@ -327,10 +324,10 @@ public:
           {
               // (scaled) Deformation gradient of the ansatz basis function
               MatrixRT defGradientU(0);
-              defGradientU[k][0] = gradients[i][0];                       // Y
-              defGradientU[k][1] = gradients[i][1];                       //X2
-  //             defGradientU[k][2] = (1.0/gamma)*gradients[i][2];           //X3
-              defGradientU[k][2] = gradients[i][2];           //X3
+              defGradientU[k][0] = jacobians[i][0][0];                       // Y
+              defGradientU[k][1] = jacobians[i][0][1];                       //X2
+  //             defGradientU[k][2] = (1.0/gamma)*jacobians[i][0][2];           //X3
+              defGradientU[k][2] = jacobians[i][0][2];           //X3
               // printmatrix(std::cout, defGradientU , "defGradientU", "--");
               defGradientU = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientU));
 
@@ -562,15 +559,14 @@ public:
     for (const auto& quadPoint : quad)
     {
       const Dune::FieldVector<double,dim>& quadPos = quadPoint.position();
-      const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+      const auto geometryJacobianInverse = geometry.jacobianInverse(quadPos);
       const double integrationElement = geometry.integrationElement(quadPos);
 
-      std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
-      localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+      std::vector<Dune::FieldMatrix<double,1,dim> > jacobians;
+      localFiniteElement.localBasis().evaluateJacobian(quadPos, jacobians);
 
-      std::vector< Dune::FieldVector< double, dim> > gradients(referenceGradients.size());
-      for (size_t i=0; i< gradients.size(); i++)
-        jacobian.mv(referenceGradients[i][0], gradients[i]);
+      for (size_t i=0; i< jacobians.size(); i++)
+        jacobians[i] = jacobians[i] * geometryJacobianInverse;
       
       //TEST
   //     std::cout << "forceTerm(element.geometry().global(quadPos)):"  << std::endl;
@@ -595,10 +591,10 @@ public:
         {
           // Deformation gradient of the ansatz basis function
           MatrixRT defGradientV(0);
-          defGradientV[k][0] = gradients[i][0];                                 // Y
-          defGradientV[k][1] = gradients[i][1];                                 // X2
-  //         defGradientV[k][2] = (1.0/gamma_)*gradients[i][2];                     // 
-          defGradientV[k][2] = gradients[i][2];                     // X3
+          defGradientV[k][0] = jacobians[i][0][0];                                 // Y
+          defGradientV[k][1] = jacobians[i][0][1];                                 // X2
+  //         defGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2];                     //
+          defGradientV[k][2] = jacobians[i][0][2];                     // X3
           
           defGradientV = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));