diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index c11727750567d796725c90ac33c43f817a0437ba..eebca42d6bbb6775fbd64204c3a99a049e90c42a 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -244,9 +244,11 @@ public:
   /** \brief Compute the stiffness matrix for one element
    *
    * \param quadRule The quadrature rule to be used
+   * \param phaseAtQuadPoint The material phase (i.e., the type of material) at each quadrature point
    */
   void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
                                      const Dune::QuadratureRule<double,dim>& quadRule,
+                                     const std::vector<int>& phaseAtQuadPoint,
                                     ElementMatrixCT& elementMatrix
                                     )
   {
@@ -262,14 +264,6 @@ public:
 
 
     // auto elasticityTensor = material_.getElasticityTensor();
-    // auto indicatorFunction = material_.getIndicatorFunction();
-    auto localIndicatorFunction = material_.getLocalIndicatorFunction();
-    localIndicatorFunction.bind(element);
-    // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
-    // auto indicatorFunction = localFunction(indicatorFunctionGVF);
-    // indicatorFunction.bind(element);
-    // Func2int indicatorFunction = *material_.getIndicatorFunction();
-
 
     // LocalBasis-Offset
     const int localPhiOffset = localView.size();
@@ -308,7 +302,7 @@ public:
       }
 
       // Compute the material phase that the current quadrature point is in
-      const auto phase = localIndicatorFunction(quadPos);
+      const auto phase = phaseAtQuadPoint[QPcounter-1];
 
       for (size_t l=0; l< dimworld; l++)
       for (size_t j=0; j < nSf; j++ )
@@ -658,9 +652,18 @@ public:
       int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
       const auto& quadRule = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
 
+      // Determine the material phase at each quadrature point
+      auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+      localIndicatorFunction.bind(element);
+
+      std::vector<int> phaseAtQuadPoint(quadRule.size());
+
+      for (std::size_t i=0; i<quadRule.size(); ++i)
+        phaseAtQuadPoint[i] = localIndicatorFunction(quadRule[i].position());
+
       ElementMatrixCT elementMatrix;
       // computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
-      computeElementStiffnessMatrix(localView, quadRule, elementMatrix);
+      computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
       
       
   //     std::cout << "local assembly time:" << Time.elapsed() << std::endl;