diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index f0d989cb933e0bfa3e05a3d647caf4d46de4bf27..a7c6f8cdf9591e737d4575dba88aeee4d484583b 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -60,8 +60,8 @@ protected:
   VectorCT phi_1_, phi_2_, phi_3_;      // Corrector phi_i coefficient vectors 
   Dune::FieldVector<double,3> m_1_, m_2_, m_3_;  // Corrector m_i coefficient vectors
 
-  MatrixRT M1_, M2_, M3_;  // (assembled) corrector matrices M_i
-  const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
+  // (assembled) corrector matrices M_i
+  std::array<MatrixRT, 3 > mContainer;
   const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
 
   // ---- Basis for R_sym^{2x2}
@@ -1094,30 +1094,27 @@ public:
       }
       // assemble M_alpha's (actually iota(M_alpha) )
 
-      // MatrixRT M_1(0), M_2(0), M_3(0);
-
-      M1_ = 0;
-      M2_ = 0;
-      M3_ = 0;
+      for (auto& matrix : mContainer)
+        matrix = 0;
 
       for(size_t i=0; i<3; i++)
       {
-          M1_ += m_1_[i]*matrixBasis_[i];
-          M2_ += m_2_[i]*matrixBasis_[i];
-          M3_ += m_3_[i]*matrixBasis_[i];
+          mContainer[0] += m_1_[i]*matrixBasis_[i];
+          mContainer[1] += m_2_[i]*matrixBasis_[i];
+          mContainer[2] += m_3_[i]*matrixBasis_[i];
       }
 
       std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
-      printmatrix(std::cout, M1_, "Corrector-Matrix M_1", "--");
-      printmatrix(std::cout, M2_, "Corrector-Matrix M_2", "--");
-      printmatrix(std::cout, M3_, "Corrector-Matrix M_3", "--");
+      printmatrix(std::cout, mContainer[0], "Corrector-Matrix M_1", "--");
+      printmatrix(std::cout, mContainer[1], "Corrector-Matrix M_2", "--");
+      printmatrix(std::cout, mContainer[2], "Corrector-Matrix M_3", "--");
       log_ << "---------- OUTPUT ----------" << std::endl;
       log_ << " --------------------" << std::endl;
-      log_ << "Corrector-Matrix M_1: \n" << M1_ << std::endl;
+      log_ << "Corrector-Matrix M_1: \n" << mContainer[0] << std::endl;
       log_ << " --------------------" << std::endl;
-      log_ << "Corrector-Matrix M_2: \n" << M2_ << std::endl;
+      log_ << "Corrector-Matrix M_2: \n" << mContainer[1] << std::endl;
       log_ << " --------------------" << std::endl;
-      log_ << "Corrector-Matrix M_3: \n" << M3_ << std::endl;
+      log_ << "Corrector-Matrix M_3: \n" << mContainer[2] << std::endl;
       log_ << " --------------------" << std::endl;
 
 
@@ -1204,9 +1201,9 @@ public:
       computeL2Norm();
       computeL2SymGrad();
 
-      std::cout<< "Frobenius-Norm of M1_: " << M1_.frobenius_norm() << std::endl;
-      std::cout<< "Frobenius-Norm of M2_: " << M2_.frobenius_norm() << std::endl;
-      std::cout<< "Frobenius-Norm of M3_: " << M3_.frobenius_norm() << std::endl;
+      std::cout<< "Frobenius-Norm of M1_: " << mContainer[0].frobenius_norm() << std::endl;
+      std::cout<< "Frobenius-Norm of M2_: " << mContainer[1].frobenius_norm() << std::endl;
+      std::cout<< "Frobenius-Norm of M3_: " << mContainer[2].frobenius_norm() << std::endl;
   }
 
   void computeL2Norm()
@@ -1379,7 +1376,7 @@ public:
           const double integrationElement = geometry.integrationElement(quadPos);
           
           
-          auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + *mContainer[b];
+          auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + mContainer[b];
 
           auto strain1 = DerPhi1(quadPos);
 
@@ -1391,7 +1388,7 @@ public:
           auto G = matrixFieldG(quadPos);
       //       auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
           
-          auto tmp = G + *mContainer[a] + strain1;
+          auto tmp = G + mContainer[a] + strain1;
 
           double energyDensity= scalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
           // double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index e7deec3720486f76f10fd09032ef1dbb5c538656..b6f3e8aff123858541694bd5eaa089bedcc79ac7 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -162,7 +162,7 @@ public:
                     const auto& quadPos = quadPoint.position();
                     const double integrationElement = geometry.integrationElement(quadPos);
                     
-                    auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a];
+                    auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + MContainer[a];
                     
                     
                     auto G1 = matrixFieldG1(quadPos);