diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index cdc6dd1a1ee9b791804f4bcdf8bd0a801302d365..f0b97b97cb7d117d46636bf89c0ac4f5137dc178 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -309,6 +309,9 @@ public:
       for (size_t i=0; i<gradients.size(); i++)
         jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
 
+      // Compute the material phase that the current quadrature point is in
+      const auto phase = localIndicatorFunction(quadPos);
+
       for (size_t l=0; l< dimworld; l++)
       for (size_t j=0; j < nSf; j++ )
       {
@@ -378,7 +381,7 @@ public:
               // printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
 
 
-              double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
+              double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,phase),sym(defGradientV));
               // 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));
@@ -419,7 +422,7 @@ public:
 
 
 
-              double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
+              double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(defGradientV));
               // 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;
@@ -468,7 +471,7 @@ public:
 
 
 
-          double energyDensityGG = scalarProduct(material_.applyElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
+          double energyDensityGG = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(basisContainer[n]));
           // double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
           // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
           // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
@@ -1473,4 +1476,4 @@ public:
 
 }; // end class
 
-#endif
\ No newline at end of file
+#endif