diff --git a/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
index 0f8d9fce2d441f0c185886d9937f2b5ef9c9ee1c..c5fbc3a77e9b505a756783fbcb1018aa0b7cc363 100644
--- a/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
+++ b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
@@ -16,6 +16,8 @@
 #include <dune/microstructure/microproblem.hh>
 
 #include <cmath>
+#include <string> 
+#include <map>
 
 
   // template <class T>
@@ -75,7 +77,9 @@ namespace Dune::GFE
 
     using Domain = typename Basis::GridView ::template Codim<0>::Geometry::GlobalCoordinate;
 
-    using IndicatorType = std::function< int(const Domain&) >;
+    using IndicatorType = std::function<int(const Domain&)>;
+    using IndicatorGridViewFunctionType = GridViewFunction<int(const Domain&), GridView>;
+    using LocalIndicatorFunctionType = LocalFunction<int(const Domain&), typename GridViewEntitySet<GridView, 0>::Element >;
 
     LocalDiscreteKirchhoffFunction& localFunction_;
 
@@ -93,7 +97,10 @@ namespace Dune::GFE
 
 
     const Python::Callable globalMicrostructureClass_;  //we use this class to create microstructure objects passed to microproblem.
-    Python::Callable LocalMicrostructureClass_;
+    Python::Reference globalMicrostructure_;
+
+    // Python::Callable LocalMicrostructureClass_;
+    std::vector<Python::Callable> LocalMicrostructureClass_;
     // const Python::Callable microstructureClass_;  //we use this class to create microstructure objects passed to microproblem.
     // Dune::ParameterTree parameterSet_;
     // Python::Module module_;
@@ -102,11 +109,16 @@ namespace Dune::GFE
     // MicroProblem microProblem_; 
 
     // std::shared_ptr<Python::Callable> anotherMicroProblem_;  
-    std::shared_ptr<MicroProblem> microProblem_;  
+    // std::shared_ptr<MicroProblem> microProblem_;  
 
+    // std::vector<std::shared_ptr<MicroProblem>> microProblem_;  
+    std::map<int,std::shared_ptr<MicroProblem> > microProblem_;  
 
 
     IndicatorType macroPhaseIndicator_;
+    mutable LocalIndicatorFunctionType macroPhaseIndicatorLocalFunction_;
+
+    int numberOfMacroPhases_; 
 
 
     /**
@@ -177,6 +189,11 @@ namespace Dune::GFE
           std::cout << "There is no class 'GlobalMicrostructure' in the Python module!" << std::endl;
           std::cerr << "Error during  DiscreteKirchhoffBendingEnergyPrestrained setup" << e << std::endl;
         }
+
+      // TEST
+      // std::cout << "Trying to access inner class directly without creating an instance:" << std::endl;
+      // Python::Callable TestClass = module_.get("LocalMicrostructure_1");
+
       return globalMicrostructureClass;
     };
 
@@ -236,39 +253,95 @@ namespace Dune::GFE
 
 
 
-
-
         //create instance of globalMicrostructureClass_
-        Python::Reference globalMicrostructure = globalMicrostructureClass_();
+        // Python::Reference globalMicrostructure = globalMicrostructureClass_();
+        globalMicrostructure_ = globalMicrostructureClass_();
 
-        int numberOfMacroPhases; // = globalMicrostructure.get("macroPhases", 1)
-        globalMicrostructure.get("macroPhases").toC<int>(numberOfMacroPhases);
-
-        std::cout << "numberOfMacroPhases: " << numberOfMacroPhases << std::endl;
+         // = globalMicrostructure.get("macroPhases", 1)
+        globalMicrostructure_.get("macroPhases").toC<int>(numberOfMacroPhases_);
+        std::cout << "numberOfMacroPhases: " << numberOfMacroPhases_ << std::endl;
 
 
 
         // Get the macroPhaseIndicator function
         try {
-          macroPhaseIndicator_ = Python::make_function<int>(Python::Callable(globalMicrostructure.get("macroPhaseIndicator")));
+          macroPhaseIndicator_ = Python::make_function<int>(Python::Callable(globalMicrostructure_.get("macroPhaseIndicator")));
+          // auto macroPhaseIndicatorGVF = Dune::Functions::makeGridViewFunction(macroPhaseIndicator_, localFunction_.gridView());
+          macroPhaseIndicatorLocalFunction_ = localFunction(Dune::Functions::makeGridViewFunction(macroPhaseIndicator_, localFunction_.gridView()));
         } catch (Dune::Exception &e) {
-          std::cout << "Piecewise constant microstructure but no macroPhaseIndicator was provided. Assuming constant microstructure" << std::endl;
+          std::cout << "Piecewise constant microstructure but no macroPhaseIndicator was provided. Assuming macroscopically constant microstructure" << std::endl;
           // Using constant True function
           macroPhaseIndicator_ = Python::make_function<int>(Python::Callable(Python::evaluate("lambda x: 1")));
         }
+        //VTK-Write macroPhaseIndicator:
+        if(parameterSet_.get<bool>("VTKwriteMacroPhaseIndicator", 0))
+          VTKwriteMacroPhaseIndicator();
+        
+        std::cout << "Getting macroPhase Indicator works" << std::endl;
+
+
+        Qhom_.resize(numberOfMacroPhases_);
+        Beff_.resize(numberOfMacroPhases_);
+      
+
+        for(size_t i=0; i<numberOfMacroPhases_; i++)
+        {
+            std::cout << "------------------- " + std::to_string(i+1) + "-th macro phase setup " + "------------------- " << std::endl;
+            LocalMicrostructureClass_.push_back(globalMicrostructure_.get("LocalMicrostructure_"+std::to_string(i+1)));
 
+            // microProblem_.push_back(std::make_shared<MicroProblem>(LocalMicrostructureClass_[i](),parameterSet_, module_));
 
+            // also store these quantities in a vector instead?
+            //---Create localMicrostructure object
+            // auto LocalMicrostructure = microProblem_[i]->getMicrostructure();
+            auto LocalMicrostructure = LocalMicrostructureClass_[i]();
 
+            // std::string testString;
+            // LocalMicrostructure.get("phase" + std::to_string(1) + "_type").toC<std::string>(testString);
+            // std::cout << "Teststring:" << testString << std::endl;
+
+
+            // bool .macroPhase2_readEffectiveQuantities
+            bool constantMicrostructure = true;
+            globalMicrostructure_.get("macroPhase" + std::to_string(i+1) + "_constantMicrostructure").toC<bool>(constantMicrostructure);
+
+            if (constantMicrostructure)
+            {
+              std::cout << "[Constant Microstructure]" << std::endl;
+                          // Check whether to read effective quantities for current phase or to compute them.
+              bool readEffectiveQuantities = false; //default
+              globalMicrostructure_.get("macroPhase" + std::to_string(i+1) + "_readEffectiveQuantities").toC<bool>(readEffectiveQuantities);
+
+              std::cout << "readEffectiveQuantities:" << readEffectiveQuantities << std::endl;
+
+              if(readEffectiveQuantities)
+              { 
+                std::cout << "(Read effective quantities for phase " + std::to_string(i+1) + ")" << std::endl;
+                LocalMicrostructure.get("effectiveQuadraticForm").toC<Dune::FieldMatrix<double,3,3>>(Qhom_[i]);
+                LocalMicrostructure.get("effectivePrestrain").toC<Dune::FieldMatrix<double,2,2>>(Beff_[i]);
+              } else {
+                std::cout << "(Solving the Micro-Problem once to obtain effective quantities ...)" << std::endl;
+
+                microProblem_.insert({i,std::make_shared<MicroProblem>(LocalMicrostructureClass_[i](),parameterSet_, module_)});
+                (microProblem_[i])->getEffectiveQuantities(Beff_[i],Qhom_[i]);
+              }
+
+              printmatrix(std::cout, Qhom_[i], "Setup effective Quadratic form (Qhom) for Phase " +std::to_string(i+1) + "as ", "--");
+              printmatrix(std::cout, Beff_[i], "Setup effective prestrain (Beff) for Phase " +std::to_string(i+1) + "as ", "--");
+            }
+            else 
+              std::cout << "[Macroscopically varying Microstructure]" << std::endl;
+
+        } //end-macroPhaseLoop
 
         // TEST Call python inner-class
-        std::cout << "Try to access inner python-class: " << std::endl;
+        // std::cout << "Try to access inner python-class: " << std::endl;
         // Python::Reference innerClass; 
-        // globalMicrostructure.get("lg").toC<Python::Reference>(innerClass); //conversion to C of Python::Reference not implemented!
-
-        // Python::Callable LocalMicrostructureClass = globalMicrostructure.get(Python::Callable("LocalMicrostructure"));
-        // Python::Callable LocalMicrostructureClass = globalMicrostructure.get("LocalMicrostructure_1");
-        LocalMicrostructureClass_ = globalMicrostructure.get("LocalMicrostructure_1");
+        // globalMicrostructure_.get("lg").toC<Python::Reference>(innerClass); //conversion to C of Python::Reference not implemented!
 
+        // Python::Callable LocalMicrostructureClass = globalMicrostructure_.get(Python::Callable("LocalMicrostructure"));
+        // Python::Callable LocalMicrostructureClass = globalMicrostructure_.get("LocalMicrostructure_1");
+        // LocalMicrostructureClass_[0] = globalMicrostructure_.get("LocalMicrostructure_1");
 
         // create instance of inner class: 
         // auto LocalMicrostructure = LocalMicrostructureClass();
@@ -276,86 +349,85 @@ namespace Dune::GFE
 
 
 
-        // int numberOfMacroPhases = parameterSet_.get<int>("numberOfMacroPhases", 1);
+        // int numberOfMacroPhases_ = parameterSet_.get<int>("numberOfMacroPhases_", 1);
+
 
-        Qhom_.resize(numberOfMacroPhases);
-        Beff_.resize(numberOfMacroPhases);
 
         //  microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(),parameterSet_, module_);
-        microProblem_ = std::make_shared<MicroProblem>(LocalMicrostructureClass_(),parameterSet_, module_);
+        // microProblem_ = std::make_shared<MicroProblem>(LocalMicrostructureClass_[0](),parameterSet_, module_);
 
-        auto LocalMicrostructure = microProblem_->getMicrostructure();
+        // auto LocalMicrostructure = microProblem_->getMicrostructure();
 
 
         // access members of inner class:
         // auto test = LocalMicrostructure.get("name");
-        std::string test;
-        LocalMicrostructure.get("name").toC<std::string>(test);
-        std::cout << "access member of inner class: " << test << std::endl;
 
-        
 
 
-        //VTK-Write macroPhaseIndicator:
-        if(parameterSet_.get<bool>("VTKwriteMacroPhaseIndicator", 0))
-          VTKwriteMacroPhaseIndicator();
+        // std::string test;
+        // LocalMicrostructure.get("name").toC<std::string>(test);
+        // std::cout << "access member of inner class: " << test << std::endl;
+
         
 
 
 
-        if(parameterSet_.get<bool>("macroscopically_varying_microstructure", 0))
-        {
-           std::cout << "MACROSCOPICALLY VARYING MICROSTRUCTURE is used! " << std::endl;
-           /*
-            *  This 'default' microstructureClass_() is just used for an initial setup.
-           */
-          //  microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(initialMacroPoint),parameterSet_, module_);
-          //  microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(),parameterSet_, module_);
-        }
-        else{
-          std::cout << "CONSTANT MICROSTRUCTURE is used with effective quantities..." << std::endl;
 
 
-          // Get Number of MacroPhases
 
-          // If possible read the effective quantities from ParameterTree
-          if(parameterSet_.get<bool>("read_effectiveQuantities_from_Parset", 1))
-          {
+        // if(parameterSet_.get<bool>("macroscopically_varying_microstructure", 0))
+        // {
+        //    std::cout << "MACROSCOPICALLY VARYING MICROSTRUCTURE is used! " << std::endl;
+        //    /*
+        //     *  This 'default' microstructureClass_() is just used for an initial setup.
+        //    */
+        //   //  microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(initialMacroPoint),parameterSet_, module_);
+        //   //  microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(),parameterSet_, module_);
+        // }
+        // else{
+        //   std::cout << "CONSTANT MICROSTRUCTURE is used with effective quantities..." << std::endl;
+
 
+        //   // Get Number of MacroPhases
 
+        //   // If possible read the effective quantities from ParameterTree
+        //   if(parameterSet_.get<bool>("read_effectiveQuantities_from_Parset", 1))
+        //   {
 
-            std::cout << "READ EFFECTIVE QUANTITIES FROM PARSET..." << std::endl;
-            // module_.get("effectiveQuadraticForm").toC<Dune::FieldMatrix<double,3,3>>(Qhom_[0]);
-            // module_.get("effectivePrestrain").toC<Dune::FieldMatrix<double,2,2>>(Beff_[0]);
 
 
-            std::cout << "Obtain effective quantitites from Microstructure-Class:" << std::endl;
+        //     std::cout << "READ EFFECTIVE QUANTITIES FROM PARSET..." << std::endl;
+        //     // module_.get("effectiveQuadraticForm").toC<Dune::FieldMatrix<double,3,3>>(Qhom_[0]);
+        //     // module_.get("effectivePrestrain").toC<Dune::FieldMatrix<double,2,2>>(Beff_[0]);
+
+
+        //     std::cout << "Obtain effective quantitites from Microstructure-Class:" << std::endl;
 
        
-            // (microProblem_->microstructure_).get("phase" + std::to_string(1) + "_type").toC<std::string>(testString);
+        //     // (microProblem_->microstructure_).get("phase" + std::to_string(1) + "_type").toC<std::string>(testString);
 
 
 
-            std::string testString;
-            LocalMicrostructure.get("phase" + std::to_string(1) + "_type").toC<std::string>(testString);
-            std::cout << "Teststring:" << testString << std::endl;
+        //     std::string testString;
+        //     LocalMicrostructure.get("phase" + std::to_string(1) + "_type").toC<std::string>(testString);
+        //     std::cout << "Teststring:" << testString << std::endl;
 
 
-            LocalMicrostructure.get("effectiveQuadraticForm").toC<Dune::FieldMatrix<double,3,3>>(Qhom_[0]);
-            LocalMicrostructure.get("effectivePrestrain").toC<Dune::FieldMatrix<double,2,2>>(Beff_[0]);
+        //     LocalMicrostructure.get("effectiveQuadraticForm").toC<Dune::FieldMatrix<double,3,3>>(Qhom_[0]);
+        //     LocalMicrostructure.get("effectivePrestrain").toC<Dune::FieldMatrix<double,2,2>>(Beff_[0]);
 
-          }
-          else{
-            std::cout << "Solving the Micro/Cell-Problem (once) to obtain effective quantities ..." << std::endl;
+        //   }
+        //   else{
+        //     std::cout << "Solving the Micro/Cell-Problem (once) to obtain effective quantities ..." << std::endl;
 
-            // microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(),parameterSet_, module_);
-            // microProblem_.updateMicrostructure(microstructureClass_()); // pass a constant microstructure object.
-            microProblem_->getEffectiveQuantities(Beff_[0],Qhom_[0]);
-          }
+        //     // microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(),parameterSet_, module_);
+        //     // microProblem_.updateMicrostructure(microstructureClass_()); // pass a constant microstructure object.
+        //     microProblem_->getEffectiveQuantities(Beff_[0],Qhom_[0]);
+        //   }
 
-          printmatrix(std::cout, Qhom_[0], "Setup effective Quadratic form (Qhom) as  ", "--");
-          printmatrix(std::cout, Beff_[0], "Setup effective prestrain (Beff) as  ", "--");
-        }
+        //   printmatrix(std::cout, Qhom_[0], "Setup effective Quadratic form (Qhom) as  ", "--");
+        //   printmatrix(std::cout, Beff_[0], "Setup effective prestrain (Beff) as  ", "--");
+        // }
           
     }
 
@@ -385,17 +457,17 @@ namespace Dune::GFE
 
 
     template<class Itype>
-    void resetEffectiveQuadraticForm(const Dune::FieldMatrix<Itype,3,3>& A)const
+    void resetEffectiveQuadraticForm(const Dune::FieldMatrix<Itype,3,3>& A, const int macroPhase) const
     {
         std::cout << "reset quadratic form ..." << std::endl;
-        convertFieldMatrix<RT>(A,Qhom_[0]);
+        convertFieldMatrix<RT>(A,Qhom_[macroPhase]);
     }
 
     template<class Itype>
-    void resetEffectivePrestrain(const Dune::FieldMatrix<Itype,2,2>& B) const
+    void resetEffectivePrestrain(const Dune::FieldMatrix<Itype,2,2>& B, const int macroPhase) const
     {
         std::cout << "reset effective prestrain ..." << std::endl;
-        convertFieldMatrix<RT>(B,Beff_[0]);
+        convertFieldMatrix<RT>(B,Beff_[macroPhase]);
     }
 
 
@@ -418,7 +490,7 @@ namespace Dune::GFE
 
     /** \brief Assemble the energy for a single element */
     virtual RT energy (const typename Basis::LocalView& localView,         
-               const  std::vector<TargetSpace>& localSolution) const       
+               const  std::vector<TargetSpace>& localSolution) const
     {
       RT energy = 0;
 
@@ -464,6 +536,11 @@ namespace Dune::GFE
       // bind to current element.
       localFunction_.bind(element);
       localForce_.bind(element);
+
+      macroPhaseIndicatorLocalFunction_.bind(element);
+
+
+
       /**
        * @brief We need to update the Coefficients of localFunction_
        *        on the current element.
@@ -515,7 +592,7 @@ namespace Dune::GFE
       // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
       // Trapezoidal-rule + evaluation on edge centers (only for testing purposes):
       // std::vector<Dune::QuadraturePoint<double,2>> quadRule = { {{0.0,0.0}, 1.0/6.0}, {{0.5,0.0}, 1.0/6.0}, {{1.0,0.0}, 1.0/6.0}, {{0.0,0.5}, 1.0/6.0}, {{0.5,0.5}, 1.0/6.0}, {{0.0,1.0}, 1.0/6.0} };
-      int count = 0;
+      // int count = 0;
       // std::cout << "Number of quadrature points: " << quadRule.size() << std::endl;
       for (auto&& quadPoint : quadRule)
       {
@@ -523,45 +600,77 @@ namespace Dune::GFE
         const auto integrationElement = geometry.integrationElement(quadPos);
 
 
+        // exit(0);
 
-        //------------------------------------------------------------------------------------
-        if(parameterSet_.get<bool>("macroscopically_varying_microstructure", 0))
+        // Determine macro phase and type 
+        int currentPhase = macroPhaseIndicatorLocalFunction_(quadPos);
+
+        bool constantMicrostructure = true;
+        globalMicrostructure_.get("macroPhase" + std::to_string(currentPhase) + "_constantMicrostructure").toC<bool>(constantMicrostructure);
+
+        if(!constantMicrostructure)
         {
-          std::cout << "Compute effective quantities at current quadrature point." << element.geometry().global(quadPos) << std::endl;
+           std::cout << "Compute effective quantities at current quadrature point: " << element.geometry().global(quadPos) << std::endl;
 
-          // std::cout << "type_name<decltype(quadPos)>():" << type_name<decltype(quadPos)>() << std::endl;
-          // std::cout << "initialMacroPoint:" << initialMacroPoint << std::endl;
+          auto microProblem = microProblem_.find(currentPhase);
 
-          // auto test = microstructureClass_(element.geometry().global(quadPos));
+          ((*microProblem).second)->updateMicrostructure(LocalMicrostructureClass_[(currentPhase-1)](element.geometry().global(quadPos))); 
+          ((*microProblem).second)->getEffectiveQuantities(Beff_[(currentPhase-1)],Qhom_[(currentPhase-1)]);
+        
+        }
+
+
+        // //------------------------------------------------------------------------------------
+        // if(parameterSet_.get<bool>("macroscopically_varying_microstructure", 0))
+        // {
+        //   std::cout << "Compute effective quantities at current quadrature point." << element.geometry().global(quadPos) << std::endl;
+
+        //   // std::cout << "type_name<decltype(quadPos)>():" << type_name<decltype(quadPos)>() << std::endl;
+        //   // std::cout << "initialMacroPoint:" << initialMacroPoint << std::endl;
+
+        //   // auto test = microstructureClass_(element.geometry().global(quadPos));
+
+        //     //USAGE:
+
+        //   //TODO 
+        //   //create instance of globalMicrostructure_Class_
+        //   // Python::Reference globalMicrostructure_ = globalMicrostructure_Class_();
+        //   // //get LocalMicrostructureClass
+        //   // Python::Callable LocalMicrostructureClass = globalMicrostructure_.get("LocalMicrostructure_1");
 
-            //USAGE:
+          
+
+        //   auto microProblem = microProblem_.find(0);
+
+        //   ((*microProblem).second)->updateMicrostructure(LocalMicrostructureClass_[0](element.geometry().global(quadPos))); 
+        //   ((*microProblem).second)->getEffectiveQuantities(Beff_[0],Qhom_[0]);
+
+        //   // (*microProblem_[0])->updateMicrostructure(LocalMicrostructureClass_[0](element.geometry().global(quadPos))); 
+        //   // (*microProblem_[0])->getEffectiveQuantities(Beff_[0],Qhom_[0]);
 
-          //TODO 
-          //create instance of globalMicrostructureClass_
-          // Python::Reference globalMicrostructure = globalMicrostructureClass_();
-          // //get LocalMicrostructureClass
-          // Python::Callable LocalMicrostructureClass = globalMicrostructure.get("LocalMicrostructure_1");
+          
+        //   //easier to simply create new microProblem object ...
 
 
-          microProblem_->updateMicrostructure(LocalMicrostructureClass_(element.geometry().global(quadPos))); 
-          microProblem_->getEffectiveQuantities(Beff_[0],Qhom_[0]);
-          // microProblem_->updateMicrostructure(microstructureClass_(element.geometry().global(quadPos))); 
-          // microProblem_->getEffectiveQuantities(Beff_,Qhom_,count);
+        //   // microProblem_->updateMicrostructure(LocalMicrostructureClass_[0](element.geometry().global(quadPos))); 
+        //   // microProblem_->getEffectiveQuantities(Beff_[0],Qhom_[0]);
+        //   // microProblem_->updateMicrostructure(microstructureClass_(element.geometry().global(quadPos))); 
+        //   // microProblem_->getEffectiveQuantities(Beff_,Qhom_,count);
 
-          count++;
+        //   // count++;
 
         
 
-          /**
-           * @brief  TODO .. 
-           * 
-           *        1.information from each quadrature point has to enter.
-           *        2. Caching of Effective Quantities
-           * 
-           */
+        //   /**
+        //    * @brief  TODO .. 
+        //    * 
+        //    *        1.information from each quadrature point has to enter.
+        //    *        2. Caching of Effective Quantities
+        //    * 
+        //    */
 
-          // exit(0);
-        }
+        //   // exit(0);
+        // }
 
 
 
@@ -669,7 +778,7 @@ namespace Dune::GFE
 
         // apply Qhom to each slice of the Hessian... 
         for(int k=0; k<3; k++)
-          term1 += applyQhom(Qhom_[0],X[k], X[k]); 
+          term1 += applyQhom(Qhom_[currentPhase-1],X[k], X[k]); 
       
 
         term1 *= weight;
@@ -681,14 +790,14 @@ namespace Dune::GFE
          * @brief Compute the second (mixed) term
          */
         FieldVector<RT,3> secondFFCoefficients = matrixToSymCoefficients(secondFF);
-        FieldVector<RT,3> BeffCoefficients = matrixToSymCoefficients(Beff_[0]);
-        RT term2 = 2.0 * weight* applyQhom(Qhom_[0],BeffCoefficients,secondFFCoefficients) ;
+        FieldVector<RT,3> BeffCoefficients = matrixToSymCoefficients(Beff_[currentPhase-1]);
+        RT term2 = 2.0 * weight* applyQhom(Qhom_[currentPhase-1],BeffCoefficients,secondFFCoefficients) ;
 
 
         /**
          * @brief Compute the third term
          */
-        RT term3 = weight * applyQhom(Qhom_[0],BeffCoefficients,BeffCoefficients);
+        RT term3 = weight * applyQhom(Qhom_[currentPhase-1],BeffCoefficients,BeffCoefficients);
 
 
         harmonicEnergy +=  term1 - term2 + term3;     
diff --git a/dune/microstructure/microproblem.hh b/dune/microstructure/microproblem.hh
index 1803445f5330942c676649ca7e174dfa8db89d58..3a4bbfdbc73ba88e85614fc92a650f2e01f0f9db 100644
--- a/dune/microstructure/microproblem.hh
+++ b/dune/microstructure/microproblem.hh
@@ -189,7 +189,7 @@ class MicroProblem
             int microGridLevel = parameterSet_.get<int>("microGridLevel", 0);
             std::array<int, dim> nElements = {(int)std::pow(2,microGridLevel) ,(int)std::pow(2,microGridLevel) ,(int)std::pow(2,microGridLevel)};
             // std::cout << "Number of Grid-Elements in each direction: " << nElements << std::endl;
-            std::cout << "Number of Grid-Elements in each direction: " << (int)std::pow(2,microGridLevel) << std::endl;
+            std::cout << "Number of Micro Grid-Elements in each direction: " << (int)std::pow(2,microGridLevel) << std::endl;
 
             return CellGridType(lower,upper,nElements);
         }
diff --git a/experiment/buckling_experiment/buckling_experiment.py b/experiment/buckling_experiment/buckling_experiment.py
index ca57eb0429ee722ba0436a3fafeb0c9806b13aeb..9b7ec0dd56a3ab4cf412ef8b854208794444ba76 100644
--- a/experiment/buckling_experiment/buckling_experiment.py
+++ b/experiment/buckling_experiment/buckling_experiment.py
@@ -256,48 +256,133 @@ parameterSet.printMicroOutput = False
 
 
 # Microstructure used:  Isotropic matrix material (phase 2) with prestrained fibers (phase 1) in the top layer aligned with the e2-direction.
-class Microstructure:
-    def __init__(self):
-        # gamma = 1.0
-        # self.macroPoint = macroPoint
-        self.gamma = 1.0    #in the future this might change depending on macroPoint.
-
-        self.phases = 2     #in the future this might change depending on macroPoint.
-        #--- Define different material phases:
-        #- PHASE 1
-        self.phase1_type="isotropic"
-        self.materialParameters_phase1 = [200, 1.0]   
-        #- PHASE 2
-        self.phase2_type="isotropic"
-        self.materialParameters_phase2 = [100, 1.0]    
+
+
+
+class GlobalMicrostructure:
+    """ Class that represents the global microstructure.
+    
+        The global microstructure can be defined individually 
+        for different (finitely many) macroscopic phases. 
+        Each macroscopic phase corresponds to a 'LocalMicrostructure_'
+        sub-class that defines the necessary data for the local 
+        micro-problem. 
+
+        Currently, there are three possibilities for the LocalMicrostructure:
+            (1) Read the effective quantities (Qhom,Beff) from this parset.
+            (Constant local microstructure)
+            (2) Solve the micro-problem once to obtain the effective quantities. 
+            (Constant local microstructure)
+            (3) Solve for micro-problem for each macroscopic quadrature point.
+            (Macroscopocally varying local microstructure)) 
+    """
+    def __init__(self,macroPoint=[0,0]):
+        self.macroPhases = 1
+
+        # define first local microstructure options.
+        self.macroPhase1_constantMicrostructure = True
+        self.macroPhase1_readEffectiveQuantities = False;
+
+
+
+    def macroPhaseIndicator(self,y): #y :macroscopic point
+            """ Indicatorfunction that determines the domains 
+                i.e. macro-phases of different local microstructures.
+            """
+            return 1;
+                
+
+        
+    # Represents the local microstructure in Phase 1
+    class LocalMicrostructure_1:  
+        def __init__(self,macroPoint=[0,0]):
+            # gamma = 1.0
+            # self.macroPoint = macroPoint
+            self.gamma = 1.0    #in the future this might change depending on macroPoint.
+
+            self.phases = 2     #in the future this might change depending on macroPoint.
+            #--- Define different material phases:
+            #- PHASE 1
+            self.phase1_type="isotropic"
+            self.materialParameters_phase1 = [200, 1.0]   
+            #- PHASE 2
+            self.phase2_type="isotropic"
+            self.materialParameters_phase2 = [100, 1.0]    
+
+            # self.effectivePrestrain= np.array([[-0.725, 0.0],
+            #                                    [0.0, -1.0]]);
+
+            # self.effectiveQuadraticForm = np.array([[19.3, 0.8, 0.0],
+            #                                         [0.8, 20.3, 0.0],
+            #                                         [0.0, 0.0, 19.3]]);
+
+
+        #--- Indicator function for material phases
+        def indicatorFunction(self,x):
+            # if (abs(x[0]) < (theta/2) and x[2] >= 0 ):
+            if (abs(x[0]) < (1.0/4.0) and x[2] <= 0 ):
+                return 1    #Phase1   
+            else :
+                return 2    #Phase2
+            
+        #--- Define prestrain function for each phase (also works with non-constant values)
+        def prestrain_phase1(self,x):
+            # return [[2, 0, 0], [0,2,0], [0,0,2]]
+            # rho = 1.0
+            rho = 1.0
+            return [[rho*1.0, 0, 0], [0,rho*1.0,0], [0,0,rho*1.0]]
+
+        def prestrain_phase2(self,x):
+            return [[0, 0, 0], [0,0,0], [0,0,0]]
+            
+
+
+
+
+
+
+# class Microstructure:
+#     def __init__(self):
+#         # gamma = 1.0
+#         # self.macroPoint = macroPoint
+#         self.gamma = 1.0    #in the future this might change depending on macroPoint.
+
+#         self.phases = 2     #in the future this might change depending on macroPoint.
+#         #--- Define different material phases:
+#         #- PHASE 1
+#         self.phase1_type="isotropic"
+#         self.materialParameters_phase1 = [200, 1.0]   
+#         #- PHASE 2
+#         self.phase2_type="isotropic"
+#         self.materialParameters_phase2 = [100, 1.0]    
 
         
-        # self.effectivePrestrain= np.array([[-0.725, 0.0],
-        #                                    [0.0, -1.0]]);
+#         # self.effectivePrestrain= np.array([[-0.725, 0.0],
+#         #                                    [0.0, -1.0]]);
 
-        # self.effectiveQuadraticForm = np.array([[19.3, 0.8, 0.0],
-        #                                         [0.8, 20.3, 0.0],
-        #                                         [0.0, 0.0, 19.3]]);
+#         # self.effectiveQuadraticForm = np.array([[19.3, 0.8, 0.0],
+#         #                                         [0.8, 20.3, 0.0],
+#         #                                         [0.0, 0.0, 19.3]]);
         
 
 
-    #--- Indicator function for material phases
-    def indicatorFunction(self,x):
-        # if (abs(x[0]) < (theta/2) and x[2] >= 0 ):
-        if (abs(x[0]) < (1.0/4.0) and x[2] <= 0 ):
-            return 1    #Phase1   
-        else :
-            return 2    #Phase2
+#     #--- Indicator function for material phases
+#     def indicatorFunction(self,x):
+#         # if (abs(x[0]) < (theta/2) and x[2] >= 0 ):
+#         if (abs(x[0]) < (1.0/4.0) and x[2] <= 0 ):
+#             return 1    #Phase1   
+#         else :
+#             return 2    #Phase2
         
-    #--- Define prestrain function for each phase (also works with non-constant values)
-    def prestrain_phase1(self,x):
-        # return [[2, 0, 0], [0,2,0], [0,0,2]]
-        # rho = 1.0
-        rho = 1.0
-        return [[rho*1.0, 0, 0], [0,rho*1.0,0], [0,0,rho*1.0]]
-
-    def prestrain_phase2(self,x):
-        return [[0, 0, 0], [0,0,0], [0,0,0]]
+#     #--- Define prestrain function for each phase (also works with non-constant values)
+#     def prestrain_phase1(self,x):
+#         # return [[2, 0, 0], [0,2,0], [0,0,2]]
+#         # rho = 1.0
+#         rho = 1.0
+#         return [[rho*1.0, 0, 0], [0,rho*1.0,0], [0,0,rho*1.0]]
+
+#     def prestrain_phase2(self,x):
+#         return [[0, 0, 0], [0,0,0], [0,0,0]]
 
 
 
diff --git a/experiment/self-folding-box/self-folding-box.py b/experiment/self-folding-box/self-folding-box.py
index 339fb4e328b1bd52da82da9751ee1d649cae8d42..4cf65174a1beb5b9f3d0557bc283a16ed1d0d0cd 100644
--- a/experiment/self-folding-box/self-folding-box.py
+++ b/experiment/self-folding-box/self-folding-box.py
@@ -160,44 +160,68 @@ parameterSet.read_effectiveQuantities_from_Parset = True # Otherwise the Micro/C
 parameterSet.printMicroOutput = False
 
 #New
-parameterSet.piecewiseConstantMicrostructure = True
+# parameterSet.piecewiseConstantMicrostructure = True
 parameterSet.VTKwriteMacroPhaseIndicator = True
 parameterSet.MacroPhaseSubsamplingRefinement = 3
 
 
 
-# Represents the global microstructure
+ 
+
 class GlobalMicrostructure:
+    """ Class that represents the global microstructure.
+    
+        The global microstructure can be defined individually 
+        for different (finitely many) macroscopic phases. 
+        Each macroscopic phase corresponds to a 'LocalMicrostructure_'
+        sub-class that defines the necessary data for the local 
+        micro-problem. 
+
+        Currently, there are three possibilities for the LocalMicrostructure:
+            (1) Read the effective quantities (Qhom,Beff) from this parset.
+            (Constant local microstructure)
+            (2) Solve the micro-problem once to obtain the effective quantities. 
+            (Constant local microstructure)
+            (3) Solve for micro-problem for each macroscopic quadrature point.
+            (Macroscopocally varying local microstructure)) 
+    """
     def __init__(self,macroPoint=[0,0]):
         self.macroPhases = 2
 
+        # define first local microstructure options.
+        self.macroPhase1_constantMicrostructure = True
         self.macroPhase1_readEffectiveQuantities = False;
-    
 
-        # self.LocalMicrostructure = self.LocalMicrostructure()
+        # define second local microstructure options.
+        self.macroPhase2_constantMicrostructure = True
+        self.macroPhase2_readEffectiveQuantities = True;
+    
 
-    def macroPhaseIndicator(self,x):
+    def macroPhaseIndicator(self,y): #y :macroscopic point
+        """ Indicatorfunction that determines the domains 
+            i.e. macro-phases of different local microstructures.
+        """
         a = 1.0/4.0
 
-        self.macroPoint = x #just for visualization purposes
+        # self.macroPoint = x #just for visualization purposes
 
-        if(self.macroPoint[1] < 1.0 and self.macroPoint[1] > 1.0-a):
+        if(y[1] < 1.0 and y[1] > 1.0-a):
             return 1
-        elif(self.macroPoint[1] < 2.0 and self.macroPoint[1] > 2.0-a):
+        elif(y[1] < 2.0 and y[1] > 2.0-a):
             return 1
-        elif(self.macroPoint[1] < 3.0 + a and self.macroPoint[1] > 3.0):
+        elif(y[1] < 3.0 + a and y[1] > 3.0):
             return 1
-        elif(self.macroPoint[0] < 0.0 and self.macroPoint[0] > -a):
+        elif(y[0] < 0.0 and y[0] > -a):
             return 1
-        elif(self.macroPoint[0] < 1.0 + a and self.macroPoint[0] > 1.0):
+        elif(y[0] < 1.0 + a and y[0] > 1.0):
             return 1 
         else:
-            return 0
+            return 2
         
 
         
     # Represents the local microstructure in Phase 1
-    class LocalMicrostructure_1:  #inner class for access testing
+    class LocalMicrostructure_1:  
         def __init__(self,macroPoint=[0,0]):
             self.name = 'Light Green'
             self.code = '024avc'
@@ -205,11 +229,10 @@ class GlobalMicrostructure:
             self.macroPoint = macroPoint
             self.gamma = 1.0    #in the future this might change depending on macroPoint.
             self.phases = 2     #in the future this might change depending on macroPoint.
-            #--- Define different material phases:
-            #- PHASE 1
+            #--- Define different local material phases:
             self.phase1_type="isotropic"
             self.materialParameters_phase1 = [200, 1.0]   
-            #- PHASE 2
+
             self.phase2_type="isotropic"
             self.materialParameters_phase2 = [100, 1.0]    
 
@@ -237,13 +260,24 @@ class GlobalMicrostructure:
 
         def prestrain_phase2(self,x):
             return [[0, 0, 0], [0,0,0], [0,0,0]]
+        
 
+    # Represents the local microstructure in Phase 2
+    class LocalMicrostructure_2:  
+        def __init__(self,macroPoint=[0,0]):
+            self.name = 'Light Green'
+            self.code = '024avc'
 
+            self.macroPoint = macroPoint
+            self.gamma = 1.0 
 
+            self.effectivePrestrain= np.array([[1.0, 0.0],
+                                            [0.0, 1.0]])
 
-
-
-
+            self.effectiveQuadraticForm = np.array([[1.0, 0.0, 0.0],
+                                                    [0.0, 1.0, 0.0],
+                                                    [0.0, 0.0, 1.0]])
+        
 
 
 # # Represents the local microstructure