diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 3c78ef0a5935d04bba949f10aef72af67c14c501..8a6f9c93887119b81f3aed98a7826b48737f75f5 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -303,6 +303,20 @@ public:
       for (size_t i=0; i< jacobians.size(); i++)
         jacobians[i] = jacobians[i] * geometryJacobianInverse;
 
+      // Compute the scaled deformation gradient for each shape function
+      std::vector<std::array<MatrixRT,dimworld> > deformationGradient(nSf);
+
+      for (size_t i=0; i < nSf; i++)
+      {
+        for (size_t k=0; k< dimworld; k++)
+        {
+          MatrixRT defGradientV(0);
+          defGradientV[k] = jacobians[i][0];
+
+          deformationGradient[i][k] = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
+        }
+      }
+
       // Compute the material phase that the current quadrature point is in
       const auto phase = localIndicatorFunction(quadPos);
 
@@ -310,29 +324,11 @@ public:
       for (size_t j=0; j < nSf; j++ )
       {
           size_t row = localView.tree().child(l).localIndex(j);
-          // (scaled)  Deformation gradient of the test basis function
-          MatrixRT defGradientV(0);
-          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
           for (size_t k=0; k < dimworld; k++)
           for (size_t i=0; i < nSf; i++)
           {
-              // (scaled) Deformation gradient of the ansatz basis function
-              MatrixRT defGradientU(0);
-              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));
-
-
-
               // auto etmp = material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)); 
               // auto etmp = elasticityTensor(defGradientU,element.geometry().global(quadPos)); 
               // auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos); 
@@ -375,7 +371,7 @@ public:
               // printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
 
 
-              double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,phase),sym(defGradientV));
+              double energyDensity= scalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),sym(deformationGradient[j][l]));
               // double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
               // double energyDensity= scalarProduct(elasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
               // double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
@@ -416,7 +412,7 @@ public:
 
 
 
-              double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(defGradientV));
+              double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(deformationGradient[j][l]));
               // double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
               // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
               // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) <<std::endl;