From 5372b560bf765ac72c5950ce66bf369885868808 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Mon, 22 Apr 2024 22:11:56 +0200
Subject: [PATCH] Add Option to suppress microProblem output/print statements

---
 dune/microstructure/CorrectorComputer.hh      |  74 ++++--
 .../EffectiveQuantitiesComputer.hh            |  11 +-
 ...scretekirchhoffbendingenergyprestrained.hh |  43 ++--
 dune/microstructure/microproblem.hh           |  36 ++-
 dune/microstructure/prestrainedMaterial.hh    |   2 +-
 .../buckling_experiment.py                    |   2 +
 experiment/macrovariationtest.py              | 217 +++++++++++-------
 7 files changed, 248 insertions(+), 137 deletions(-)

diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 331ce234..0be324e8 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -133,14 +133,18 @@ public:
   {
       Dune::Timer StiffnessTimer;
       assembleCellStiffness(stiffnessMatrix_);
-      std::cout << "Stiffness assembly Timer: " << StiffnessTimer.elapsed() << std::endl;
+      if (parameterSet_.get<bool>("printMicroOutput ", false))
+        std::cout << "Stiffness assembly Timer: " << StiffnessTimer.elapsed() << std::endl;
+    
+
 
 
       Dune::Timer LoadVectorTimer;
       assembleCellLoad(load_alpha1_ ,x3G_1_);
       assembleCellLoad(load_alpha2_ ,x3G_2_);
       assembleCellLoad(load_alpha3_ ,x3G_3_);
-      std::cout << "Load vector assembly Timer: " << LoadVectorTimer.elapsed() << std::endl;
+      if (parameterSet_.get<bool>("printMicroOutput ", false))
+        std::cout << "Load vector assembly Timer: " << LoadVectorTimer.elapsed() << std::endl;
   };
 
 
@@ -151,7 +155,7 @@ public:
 
     void updateMaterial(std::shared_ptr<Material> mat)
     {
-        std::cout << "updateMaterial of CorrectorComputer" << std::endl;
+        // std::cout << "updateMaterial of CorrectorComputer" << std::endl;
         material_ = mat;
     }
 
@@ -253,7 +257,7 @@ public:
       for (int k = 0; k<3; k++)
         nb.add(row[k],row[k]);
     }
-    std::cout << "finished occupation pattern\n";
+    // std::cout << "finished occupation pattern\n";
   }
 
 
@@ -450,12 +454,16 @@ public:
     bool cacheElementMatrices = parameterSet_.get<bool>("cacheElementMatrices", true);
     // bool cacheElementMatrices = false;
     // if(parameterSet_.get<bool>("cacheElementMatrices", true))
-    if (cacheElementMatrices)
+    if (parameterSet_.get<bool>("printMicroOutput ", false))
     {
-      std::cout << "Use Caching of Element matrices" << std::endl;
+      if (cacheElementMatrices)
+      {
+        std::cout << "Use Caching of Element matrices" << std::endl;
+      }
     }
 
 
+
     std::map<std::vector<int>, ElementMatrixCT> elementMatrixCache;
 
     for (const auto& element : elements(basis_.gridView()))
@@ -616,7 +624,8 @@ public:
 
   void setOneBaseFunctionToZero()
   {
-    std::cout << "Setting one Basis-function to zero" << std::endl;
+    if (parameterSet_.get<bool>("printMicroOutput ", false))
+      std::cout << "Setting one Basis-function to zero" << std::endl;
 
   //   constexpr int dim = Basis::LocalView::Element::dimension;
     
@@ -767,11 +776,13 @@ public:
       x_2_ = load_alpha2_;
       x_3_ = load_alpha3_;
 
-      std::cout << "start corrector solver..." << std::endl;
+      if (parameterSet_.get<bool>("printMicroOutput ", false))
+        std::cout << "start corrector solver..." << std::endl;
       Dune::Timer SolverTimer;
       if (Solvertype==1)  // CG - SOLVER
       {
-          std::cout << "------------ CG - Solver ------------" << std::endl;
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+            std::cout << "------------ CG - Solver ------------" << std::endl;
           Dune::MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);
 
           // Sequential incomplete LU decomposition as the preconditioner
@@ -786,15 +797,19 @@ public:
                                   true    // Try to estimate condition_number
                                   );   // verbosity of the solver
           Dune::InverseOperatorResult statistics;
-          std::cout << "solve linear system for x_1.\n";
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+             std::cout << "solve linear system for x_1.\n";
           solver.apply(x_1_, load_alpha1_, statistics);
-          std::cout << "solve linear system for x_2.\n";
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+             std::cout << "solve linear system for x_2.\n";
           solver.apply(x_2_, load_alpha2_, statistics);
-          std::cout << "solve linear system for x_3.\n";
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+             std::cout << "solve linear system for x_3.\n";
           solver.apply(x_3_, load_alpha3_, statistics);
           // log_ << "Solver-type used: " <<" CG-Solver" << std::endl;
 
-          if(Solver_verbosity > 0)
+
+          if(parameterSet_.get<bool>("printMicroOutput ", false) && Solver_verbosity > 0)
           {
             std::cout << "statistics.converged " << statistics.converged << std::endl;
             std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
@@ -804,7 +819,8 @@ public:
       ////////////////////////////////////////////////////////////////////////////////////
       else if (Solvertype==2)  // GMRES - SOLVER
       {
-          std::cout << "------------ GMRES - Solver ------------" << std::endl;
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+             std::cout << "------------ GMRES - Solver ------------" << std::endl;
           // Turn the matrix into a linear operator
           Dune::MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_);
 
@@ -829,7 +845,7 @@ public:
           solver.apply(x_2_, load_alpha2_, statistics);
           solver.apply(x_3_, load_alpha3_, statistics);
           // log_ << "Solver-type used: " <<" GMRES-Solver" << std::endl;
-          if(Solver_verbosity > 0)
+          if(parameterSet_.get<bool>("printMicroOutput ", false) &&  Solver_verbosity > 0)
           {
             std::cout << "statistics.converged " << statistics.converged << std::endl;
             std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
@@ -839,7 +855,8 @@ public:
       ////////////////////////////////////////////////////////////////////////////////////
       else if ( Solvertype==3)// QR - SOLVER
       {
-          std::cout << "------------ QR - Solver ------------" << std::endl;
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+             std::cout << "------------ QR - Solver ------------" << std::endl;
           // log_ << "solveLinearSystems: We use QR solver.\n";
           //TODO install suitesparse
           Dune::SPQR<MatrixCT> sPQR(stiffnessMatrix_);
@@ -863,7 +880,7 @@ public:
           // std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
           // std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
           // log_ << "Solver-type used: " <<" QR-Solver" << std::endl;
-          if(Solver_verbosity > 0)
+          if(parameterSet_.get<bool>("printMicroOutput ", false) &&  Solver_verbosity > 0)
           {
             std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
             std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
@@ -874,7 +891,8 @@ public:
       ////////////////////////////////////////////////////////////////////////////////////
       else if (Solvertype==4)// UMFPACK - SOLVER
       {
-          std::cout << "------------ UMFPACK - Solver ------------" << std::endl;
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+              std::cout << "------------ UMFPACK - Solver ------------" << std::endl;
           // log_ << "solveLinearSystems: We use UMFPACK solver.\n";
 
           if(basis_.dimension() > 1e5)
@@ -907,7 +925,8 @@ public:
       ////////////////////////////////////////////////////////////////////////////////////
       else if (Solvertype==5)// CHOLDMOD - SOLVER
       {
-          std::cout << "------------ CHOLMOD - Solver ------------" << std::endl;
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+              std::cout << "------------ CHOLMOD - Solver ------------" << std::endl;
           // log_ << "solveLinearSystems: We use CHOLMOD solver.\n";
 
           if(basis_.dimension() > 1e5)
@@ -937,8 +956,8 @@ public:
   //         std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
           // log_ << "Solver-type used: " <<" CHOLMOD-Solver" << std::endl;
       }
-      std::cout <<  "--- Corrector computation finished ---" << std::endl;
-      std::cout << "Time required for solving:" << SolverTimer.elapsed() << std::endl;
+      if (parameterSet_.get<bool>("printMicroOutput ", false))
+        std::cout << "Corrector computation finished. Time required:" << SolverTimer.elapsed() << std::endl;
 
       ////////////////////////////////////////////////////////////////////////////////////
       // Extract phi_alpha  &  M_alpha coefficients
@@ -1004,14 +1023,16 @@ public:
       if(substract_integralMean)
       {
           Dune::Timer integralMeanTimer;
-          std::cout << " --- subtracting integralMean --- " << std::endl;
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+              std::cout << " --- subtracting integralMean --- " << std::endl;
           subtractIntegralMean(phi_1_);
           subtractIntegralMean(phi_2_);
           subtractIntegralMean(phi_3_);
           subtractIntegralMean(x_1_);
           subtractIntegralMean(x_2_);
           subtractIntegralMean(x_3_);
-          std::cout << "substract integral mean took " << integralMeanTimer.elapsed() << " seconds " << std::endl;
+          if (parameterSet_.get<bool>("printMicroOutput ", false))
+              std::cout << "substract integral mean took " << integralMeanTimer.elapsed() << " seconds " << std::endl;
           //////////////////////////////////////////
           // Check Integral-mean again:
           //////////////////////////////////////////
@@ -1058,7 +1079,9 @@ public:
   {
       std::string baseName = parameterSet_.get("baseName", "CellProblem-result");
       std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/" + baseName; 
-      std::cout << "vtkOutputName:" << vtkOutputName << std::endl;
+      if (parameterSet_.get<bool>("printMicroOutput ", false))
+        std::cout << "Write Correctors to VTK with Filename:" << vtkOutputName << std::endl;
+      
 
       int subsamplingRefinement = parameterSet_.get<int>("subsamplingRefinement", 2);
       Dune::SubsamplingVTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView(), Dune::refinementLevels(subsamplingRefinement));
@@ -1074,7 +1097,8 @@ public:
           Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_3_),
           Dune::VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
       vtkWriter.write(vtkOutputName  + "-level"+ std::to_string(level));
-      std::cout << "wrote Corrector-VTK data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;
+      if (parameterSet_.get<bool>("printMicroOutput ", false))
+          std::cout << "wrote Corrector-VTK data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;
   }
 
   // -----------------------------------------------------------------
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index 33b5c526..ddfad687 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -80,8 +80,12 @@ public:
   // --- Compute Effective Quantities
     void computeEffectiveQuantities()
     {
-        std::cout << "starting effective quantities computation..." << std::endl;
         Dune::Timer effectiveQuantitiesTimer;
+        Dune::ParameterTree parameterSet = correctorComputer_->getParameterSet();
+
+        if (parameterSet.get<bool>("printMicroOutput ", false))
+            std::cout << "starting effective quantities computation..." << std::endl;
+        
 
         auto MContainer = correctorComputer_->getMcontainer();
         auto MatrixBasisContainer = correctorComputer_->getMatrixBasiscontainer();
@@ -89,7 +93,7 @@ public:
 
         auto gamma = correctorComputer_->getGamma();
         auto basis = *correctorComputer_->getBasis();
-        Dune::ParameterTree parameterSet = correctorComputer_->getParameterSet();
+        
 
 		shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_->getCorr_phi1(), 
                                             correctorComputer_->getCorr_phi2(),
@@ -194,7 +198,8 @@ public:
         // log << "------------------------ " << std::endl;
 
         //   std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
-        std::cout << "Time required to solve for effective quantities: " << effectiveQuantitiesTimer.elapsed() << std::endl;
+        if (parameterSet.get<bool>("printMicroOutput ", false))
+            std::cout << "Time required to solve for effective quantities: " << effectiveQuantitiesTimer.elapsed() << std::endl;
         return ;
     }
 
diff --git a/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
index bbe5e804..c4f22217 100644
--- a/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
+++ b/dune/microstructure/energies/discretekirchhoffbendingenergyprestrained.hh
@@ -189,9 +189,25 @@ namespace Dune::GFE
          * @brief If the microstructure is constant on a macroscopic level,
          *        we read the effective quantitites from the ParameterTree 
          *        and set them up once here.
+         * 
+         *        If no effective quanitities are provided, we solve the 
+         *        micro problem once to obtain them.
+         * 
+         * 
+         *        If the microscture is varying on a macroscopic scale,
+         *        the microProblem is initialized with a default value
+         *        and gets updated for each each macro point in the quadrature loop.
          */
-        if(!(parameterSet_.get<bool>("macroscopically_varying_microstructure", 0)))
+        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;
 
           // If possible read the effective quantities from ParameterTree
@@ -212,10 +228,6 @@ namespace Dune::GFE
           printmatrix(std::cout, Qhom_, "Setup effective Quadratic form (Qhom) as  ", "--");
           printmatrix(std::cout, Beff_, "Setup effective prestrain (Beff) as  ", "--");
         }
-        else{
-           std::cout << "MACROSCOPICALLY VARYING MICROSTRUCTURE is used! " << std::endl;
-           microProblem_ = std::make_shared<MicroProblem>(microstructureClass_(initialMacroPoint),parameterSet_, module_); //TODO Replace 0  with Domain input type. using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-        }
           
     }
 
@@ -281,7 +293,9 @@ namespace Dune::GFE
 
 #if 1
       // TODO: More serious order selection
-      int quadOrder = 3;
+      // int quadOrder = 3;
+      int quadOrder = 2;
+      // int quadOrder = 1;  //too low. bad results.
       // int quadOrder = 6;
 #else
       if (localFiniteElement.localBasis().order()==1)
@@ -350,6 +364,8 @@ 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;
+      // std::cout << "Number of quadrature points: " << quadRule.size() << std::endl;
       for (auto&& quadPoint : quadRule)
       {
         const auto& quadPos = quadPoint.position();
@@ -360,20 +376,19 @@ namespace Dune::GFE
         //------------------------------------------------------------------------------------
         if(parameterSet_.get<bool>("macroscopically_varying_microstructure", 0))
         {
-          std::cout << "Compute effective quantities at current quadrature point." << quadPos << std::endl;
-
-          std::cout << "type_name<decltype(quadPos)>():" << type_name<decltype(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;
 
-          std::cout << "initialMacroPoint:" << initialMacroPoint << std::endl;
-
-          auto test = microstructureClass_(quadPos);
+          // auto test = microstructureClass_(element.geometry().global(quadPos));
 
             //USAGE:
-          microProblem_->updateMicrostructure(microstructureClass_(quadPos)); 
+          microProblem_->updateMicrostructure(microstructureClass_(element.geometry().global(quadPos))); 
           microProblem_->getEffectiveQuantities(Beff_,Qhom_);
+          // microProblem_->getEffectiveQuantities(Beff_,Qhom_,count);
 
-
+          count++;
 
         
 
diff --git a/dune/microstructure/microproblem.hh b/dune/microstructure/microproblem.hh
index ee677779..c527d61f 100644
--- a/dune/microstructure/microproblem.hh
+++ b/dune/microstructure/microproblem.hh
@@ -225,17 +225,21 @@ class MicroProblem
 
 
 
+
+        // template<class Btype,class Qtype>
+        // void getEffectiveQuantities(Btype& B, Qtype& Q, int count = 0)
         template<class Btype,class Qtype>
         void getEffectiveQuantities(Btype& B, Qtype& Q)
         {
-            std::cout << "getting effetive quantities.. " << std::endl;
+            Dune::Timer microProblemTimer;
+            // std::cout << "getting effective quantities.. " << std::endl;
 
             // --- Get scale ratio 
             // double gamma = parameterSet_.get<double>("gamma",1.0); 
             double gamma;
             microstructure_.get("gamma").toC<double>(gamma);
-            std::cout << "microstructure gamma:" << gamma << std::endl;
-            std::cout << "Scale ratio (gamma) set to : " << gamma << std::endl;
+            if (parameterSet_.get<bool>("printMicroOutput ", false))
+              std::cout << "Microstructure scale ratio (gamma) set to : " << gamma << std::endl;
 
 
             // /**
@@ -265,16 +269,21 @@ class MicroProblem
             // material.writeVTKMaterialFunctions(parameterSet_.get<int>("gridLevel", 0));
             // std::cout << "Part1 works." << std::endl;
 
+
             if (parameterSet_.get<bool>("write_materialFunctions", false))
                material_->writeVTKMaterialFunctions(parameterSet_.get<int>("gridLevel", 0));
 
-            std::cout << "Material setup done" << std::endl;
+            // exit(0);
+
+            // std::cout << "Material setup done" << std::endl;
 
             // // Compute Correctors
             // auto correctorComputer = CorrectorComputer(basis_, material, log, parameterSet_);
-            std::cout << "Starting Corrector assembly" << std::endl;
+            // #if 0
+
+            // std::cout << "Starting Corrector assembly" << std::endl;
             correctorComputer_->assemble();
-            std::cout << "Starting Corrector solve..." << std::endl;
+            // std::cout << "Starting Corrector solve..." << std::endl;
             correctorComputer_->solve();
 
 
@@ -292,8 +301,19 @@ class MicroProblem
             auto Beff = CoefficientsToSymMatrix(Beffvec);
             Q = effectiveQuantitiesComputer_.getQeff();
             B = CoefficientsToSymMatrix(Beffvec);
-            printmatrix(std::cout, Q, "Matrix Qeff", "--");
-            printmatrix(std::cout, B, "Beff", "--");
+
+
+            if (parameterSet_.get<bool>("printMicroOutput ", false))
+            {
+                printmatrix(std::cout, Q, "Matrix Qeff", "--");
+                printmatrix(std::cout, B, "Beff", "--");
+            }
+
+            std::cout << "Micro-Problem took " << microProblemTimer.elapsed() << " seconds." << std::endl;
+
+            // #endif
+
+            // exit(0);
         }
 
 };
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 3c513ed5..470b013d 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -119,7 +119,7 @@ public:
   
   void updateMicrostructure(const Python::Reference microstructure) 
   {
-      std::cout << "updateMicrostrucrue of prestrainedMaterial" << std::endl;
+      // std::cout << "updateMicrostrucrue of prestrainedMaterial" << std::endl;
       microstructure_ = microstructure;
   }
 
diff --git a/experiment/buckling_experiment/buckling_experiment.py b/experiment/buckling_experiment/buckling_experiment.py
index a785f5c4..1fb8ee00 100644
--- a/experiment/buckling_experiment/buckling_experiment.py
+++ b/experiment/buckling_experiment/buckling_experiment.py
@@ -118,6 +118,8 @@ parameterSet.prestrainFlag = True
 parameterSet.macroscopically_varying_microstructure = False
 parameterSet.read_effectiveQuantities_from_Parset = False # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
 
+parameterSet.printMicroOutput = False
+
 # parameterSet.read_effectiveQuantities_from_Parset = True
 effectivePrestrain= np.array([[1.0, 0.0],
                               [0.0, 1.0]]);
diff --git a/experiment/macrovariationtest.py b/experiment/macrovariationtest.py
index d10101e5..45aaf0b3 100644
--- a/experiment/macrovariationtest.py
+++ b/experiment/macrovariationtest.py
@@ -28,10 +28,10 @@ nY=8
 
 parameterSet.structuredGrid = 'simplex'
 parameterSet.lower = '0 0'
-parameterSet.upper = '4 1'
+parameterSet.upper = '1 1'
 parameterSet.elements = str(nX)+' '+  str(nY)
 
-parameterSet.numLevels = 1   #good results only for refinementLevel >=4
+parameterSet.numLevels = 2   #good results only for refinementLevel >=4
 
 #############################################
 #  Options
@@ -77,7 +77,7 @@ parameterSet.writeDOFvector = 0
 #  Dirichlet boundary indicator
 #############################################
 def dirichlet_indicator(x) :
-    if( (x[0] <= 0.01) or (x[0] >= 3.99)):
+    if( (x[0] <= 0.01) or (x[0] >= 0.99)):
         return True
     else:
         return False
@@ -177,8 +177,17 @@ class Microstructure:
         self.macroPhases = [] #store macro phase numbers.
         self.macroPhaseCount = 0
 
-        self.a = (1-macroPoint[0])*(2-np.sqrt(3))/2.0
-        self.b = macroPoint[0] * (2-np.sqrt(3))/2
+        self.a = (1.0-macroPoint[0])*(2.0-np.sqrt(3.0))/2.0
+        self.b = macroPoint[0] * (2.0-np.sqrt(3.0))/2.0
+
+        # self.a = (0.5)*(2.0-np.sqrt(3))
+        # self.b = 0.5 * (2.0-np.sqrt(3))
+
+        # self.a = 0
+        # self.b = (2.0-np.sqrt(3))/2
+
+        # print('self.a', self.a)
+        # print('self.b', self.b)
 
     # --- Indicator function for material phases
     # def indicatorFunction(self,x):
@@ -206,10 +215,14 @@ class Microstructure:
 
         #shift domain
         x = [y[0]+0.5,y[1]+0.5]
+        # indicator = (( (self.one_norm(x[0],x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0],x[1])) )\
+        #           | (( (self.one_norm(x[0]-1,x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0]-1,x[1])))\
+        #           | (x[0] < self.a/2.0) | (x[0] >  1-self.a/2.0) | (x[1] < self.a/2.0) | (x[1] >  1-self.a/2.0)  
+    
         indicator = (( (self.one_norm(x[0],x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0],x[1])) )\
-                  | (( (self.one_norm(x[0]-1,x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0]-1,x[1])))\
-                  | (x[0] < self.a/2.0) | (x[0] >  1-self.a/2.0) | (x[1] < self.a/2.0) | (x[1] >  1-self.a/2.0)  
-                
+                | (( (self.one_norm(x[0]-1,x[1])) < 1 + self.b/2.0) & (1 - self.b/2.0 < self.one_norm(x[0]-1,x[1])))\
+                | (x[0] < self.a/2.0) | (x[0] >  1-self.a/2.0) | (x[1] < self.a/2.0) | (x[1] >  1-self.a/2.0)  
+            
         # indicator = (x[0] < self.b/2.0 )
         # indicator = (np.abs(x[0]) + np.abs(x[1]) < 1 + self.b/2.0)
         # indicator = (self.one_norm(x[0],x[1]) < 1 + self.b/2.0)
@@ -232,7 +245,7 @@ class Microstructure:
         
     #--- Define prestrain function for each phase (also works with non-constant values)
     def prestrain_phase1(self,x):
-        return [[5, 0, 0], [0,5,0], [0,0,5]]
+        return [[0, 0, 0], [0,0,0], [0,0,0]]
         # rho = 0 
         # # rho = 1.0
         # return [[rho*1.0, 0, 0], [0,rho*1.0,0], [0,0,rho*1.0]]
@@ -276,18 +289,29 @@ class Microstructure:
 #  Microstructure
 ############################################
 parameterSet.prestrainFlag = True
-parameterSet.macroscopically_varying_microstructure = True
+parameterSet.macroscopically_varying_microstructure = False
+parameterSet.read_effectiveQuantities_from_Parset = True # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
+# parameterSet.macroscopically_varying_microstructure = True
+# parameterSet.read_effectiveQuantities_from_Parset = False # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
 
 
+parameterSet.printMicroOutput = False
 
-parameterSet.read_effectiveQuantities_from_Parset = False # Otherwise the Micro/Cell-Problem is solved once to obtain the quantities.
-effectivePrestrain= np.array([[1.0, 0.0],
+
+effectivePrestrain= np.array([[0.75, 0.0],
                               [0.0, 1.0]]);
 
-effectiveQuadraticForm = np.array([[1.0, 0.0, 0.0],
+effectiveQuadraticForm = np.array([[3.0, 0.0, 0.0],
                                    [0.0, 1.0, 0.0],
-                                   [0.0, 0.0, 1.0]]);
+                                   [0.0, 0.0, 0.8]]);
+
 
+# effectivePrestrain= np.array([[1.0, 0.0],
+#                               [0.0, 1.0]]);
+
+# effectiveQuadraticForm = np.array([[1.0, 0.0, 0.0],
+#                                    [0.0, 1.0, 0.0],
+#                                    [0.0, 0.0, 1.0]]);
 #############################################
 #  Initial iterate function
 #############################################
@@ -298,35 +322,57 @@ effectiveQuadraticForm = np.array([[1.0, 0.0, 0.0],
 #     return ((1,0),
 #             (0,1),
 #             (0,0))
-#Rotation:
-def R(beta):
-    return  [[math.cos(beta),0],
-            [0,1],
-            [-math.sin(beta),0]]
+# #Rotation:
+# def R(beta):
+#     return  [[math.cos(beta),0],
+#             [0,1],
+#             [-math.sin(beta),0]]
+
+
+# def f(x):
+#     a = 0.5
+#     if(x[0] <= 0.01):
+#         return [x[0], x[1], 0]
+#     elif(x[0] >= 3.99):
+#         return [x[0] - a, x[1], 0]
+#     else:
+#         return [x[0], x[1], 0]
+
+
+# beta = math.pi/4.0
+
+# def df(x):
+#     a = 0.5
+#     if(x[0] <= 0.01):
+#         return R(-1*beta)
+#     elif(x[0] >= 3.99):
+#         return R(beta)
+#     else:
+#         return ((1,0),
+#                 (0,1),
+#                 (0,0))
+    
+
+
+#RUMPF EXAMPLE
 
 
 def f(x):
-    a = 0.5
+    a = (3.0/16.0)
     if(x[0] <= 0.01):
-        return [x[0], x[1], 0]
-    elif(x[0] >= 3.99):
+        return [x[0]+a, x[1], 0]
+    elif(x[0] >= 0.99):
         return [x[0] - a, x[1], 0]
     else:
         return [x[0], x[1], 0]
 
 
-beta = math.pi/4.0
 
 def df(x):
-    a = 0.5
-    if(x[0] <= 0.01):
-        return R(-1*beta)
-    elif(x[0] >= 3.99):
-        return R(beta)
-    else:
-        return ((1,0),
-                (0,1),
-                (0,0))
+    return ((1,0),
+            (0,1),
+            (0,0))
+
 
 
 
@@ -355,10 +401,10 @@ fdf = (f, df)
 #############################################
 #  Force
 ############################################
-parameterSet.assemble_force_term = False
+parameterSet.assemble_force_term = True
 
 def force(x):
-    return [0, 0, 1e-5]
+    return [0, 0, 1e-2]
 
 
 #############################################
@@ -383,62 +429,62 @@ def force(x):
 # Microstructure used: Parametrized Laminate.
 
 # --- Choose scale ratio gamma:
-parameterSet.gamma = 1.0
-# --- Number of material phases
-parameterSet.Phases = 4 
-
-#--- Indicator function for material phases
-def indicatorFunction(x):
-    theta=0.25
-    factor=1
-    if (abs(x[0]) < (theta/2) and x[2] < 0 ):
-        return 1    #Phase1
-    elif (abs(x[0]) > (theta/2) and x[2] > 0 ):
-        return 2    #Phase2
-    elif (abs(x[0]) < (theta/2) and x[2] > 0 ):
-        return 3    #Phase3
-    else :
-        return 4    #Phase4
+# parameterSet.gamma = 1.0
+# # --- Number of material phases
+# parameterSet.Phases = 4 
+
+# #--- Indicator function for material phases
+# def indicatorFunction(x):
+#     theta=0.25
+#     factor=1
+#     if (abs(x[0]) < (theta/2) and x[2] < 0 ):
+#         return 1    #Phase1
+#     elif (abs(x[0]) > (theta/2) and x[2] > 0 ):
+#         return 2    #Phase2
+#     elif (abs(x[0]) < (theta/2) and x[2] > 0 ):
+#         return 3    #Phase3
+#     else :
+#         return 4    #Phase4
     
-########### Options for material phases: #################################
-#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
-#########################################################################
-## Notation - Parameter input :
-# isotropic (Lame parameters) : [mu , lambda]
-#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]   # see https://en.wikipedia.org/wiki/Poisson%27s_ratio with x=1,y=2,z=3
-# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
-# general_anisotropic         : full compliance matrix C
-######################################################################
+# ########### Options for material phases: #################################
+# #     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+# #########################################################################
+# ## Notation - Parameter input :
+# # isotropic (Lame parameters) : [mu , lambda]
+# #         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]   # see https://en.wikipedia.org/wiki/Poisson%27s_ratio with x=1,y=2,z=3
+# # transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# # general_anisotropic         : full compliance matrix C
+# ######################################################################
 
-#--- Define different material phases:
-#- PHASE 1
-parameterSet.phase1_type="isotropic"
-materialParameters_phase1 = [2.0, 0]   
+# #--- Define different material phases:
+# #- PHASE 1
+# parameterSet.phase1_type="isotropic"
+# materialParameters_phase1 = [2.0, 0]   
 
-#- PHASE 2
-parameterSet.phase2_type="isotropic"
-materialParameters_phase2 = [1.0, 0]   
+# #- PHASE 2
+# parameterSet.phase2_type="isotropic"
+# materialParameters_phase2 = [1.0, 0]   
 
-#- PHASE 3
-parameterSet.phase3_type="isotropic"
-materialParameters_phase3 = [2.0, 0]
+# #- PHASE 3
+# parameterSet.phase3_type="isotropic"
+# materialParameters_phase3 = [2.0, 0]
 
-#- PHASE 4
-parameterSet.phase4_type="isotropic"
-materialParameters_phase4 = [1.0, 0]
+# #- PHASE 4
+# parameterSet.phase4_type="isotropic"
+# materialParameters_phase4 = [1.0, 0]
 
-#--- Define prestrain function for each phase (also works with non-constant values)
-def prestrain_phase1(x):
-    return [[2, 0, 0], [0,2,0], [0,0,2]]
+# #--- Define prestrain function for each phase (also works with non-constant values)
+# def prestrain_phase1(x):
+#     return [[2, 0, 0], [0,2,0], [0,0,2]]
 
-def prestrain_phase2(x):
-    return [[1, 0, 0], [0,1,0], [0,0,1]]
+# def prestrain_phase2(x):
+#     return [[1, 0, 0], [0,1,0], [0,0,1]]
 
-def prestrain_phase3(x):
-    return [[0, 0, 0], [0,0,0], [0,0,0]]
+# def prestrain_phase3(x):
+#     return [[0, 0, 0], [0,0,0], [0,0,0]]
 
-def prestrain_phase4(x):
-    return [[0, 0, 0], [0,0,0], [0,0,0]]
+# def prestrain_phase4(x):
+#     return [[0, 0, 0], [0,0,0], [0,0,0]]
 
 
 
@@ -466,9 +512,9 @@ parameterSet.Solver_verbosity = 0  #(default = 2)  degree of information for sol
 #  Write/Output options      #(default=false)
 #############################################
 # --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
-parameterSet.write_materialFunctions = 1   # VTK indicator function for material/prestrain definition
+parameterSet.write_materialFunctions = 0   # VTK indicator function for material/prestrain definition
 #parameterSet.write_prestrainFunctions = 1  # VTK norm of B (currently not implemented)
-parameterSet.MaterialSubsamplingRefinement = 4
+parameterSet.MaterialSubsamplingRefinement = 1
 
 # --- (Additional debug output)
 parameterSet.print_debug = 0  #(default=false)
@@ -476,13 +522,12 @@ parameterSet.print_debug = 0  #(default=false)
 parameterSet.print_corrector_matrices = 0
 
 # --- Write Correctos to VTK-File:  
-parameterSet.write_VTK = 1
+parameterSet.write_VTK = 0
 
 # --- Use caching of element matrices:  
 parameterSet.cacheElementMatrices = 1
 
-# The grid can be refined several times for a higher resolution in the VTK-file.
-# parameterSet.subsamplingRefinement = 2  #TODO New Subsampling for cell problem
+
 
 # --- check orthogonality (75) from paper: 
 parameterSet.write_checkOrthogonality = 0
-- 
GitLab