src/Cell-Problem.cc

+#include <config.h>
+#include <array>
+#include <vector>
+#include <fstream>
+
+#include <iostream>
+#include <dune/common/indices.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/math.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/spqr.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/periodicbasis.hh>
+
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/microstructure/prestrain_material_geometry.hh>
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/vtk_filler.hh>    //TEST
+
+
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+
+// #include <dune/fufem/discretizationerror.hh>
+// #include <boost/multiprecision/cpp_dec_float.hpp>
+// #include <boost/any.hpp>
+
+#include <any>
+#include <variant>
+#include <string>
+#include <iomanip>
+
+using namespace Dune;
+using namespace MatrixOperations;
+
+
+//////////////////////////////////////////////////////////////////////
+// Helper functions for Table-Output
+//////////////////////////////////////////////////////////////////////
+/*! Center-aligns string within a field of width w. Pads with blank spaces
+    to enforce alignment. */
+std::string center(const std::string s, const int w) {
+    std::stringstream ss, spaces;
+    int padding = w - s.size();                 // count excess room to pad
+    for(int i=0; i<padding/2; ++i)
+        spaces << " ";
+    ss << spaces.str() << s << spaces.str();    // format with padding
+    if(padding>0 && padding%2!=0)               // if odd #, add 1 space
+        ss << " ";
+    return ss.str();
+}
+
+/* Convert double to string with specified number of places after the decimal
+   and left padding. */
+template<class type>
+std::string prd(const type x, const int decDigits, const int width) {
+    std::stringstream ss;
+//     ss << std::fixed << std::right;
+    ss << std::scientific << std::right;                     // Use scientific Output!
+    ss.fill(' ');        // fill space around displayed #
+    ss.width(width);     // set  width around displayed #
+    ss.precision(decDigits); // set # places after decimal
+    ss << x;
+    return ss.str();
+}
+
+
+
+
+
+template<class Basis>
+auto arbitraryComponentwiseIndices(const Basis& basis,
+                                   const int elementNumber,
+                                   const int leafIdx
+                                   )
+{
+  // (Local Approach -- works for non Lagrangian-Basis too)
+  // Input  : elementNumber & localIdx
+  // Output : determine all Component-Indices that correspond to that basis-function
+  auto localView = basis.localView();
+
+  FieldVector<int,3> row;
+  int elementCounter = 0;
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    if(elementCounter == elementNumber)             // get arbitraryIndex(global) for each Component        ..alternativ:gridView.indexSet
+    {
+      localView.bind(element);
+
+      for (int k = 0; k < 3; k++)
+      {
+        auto rowLocal = localView.tree().child(k).localIndex(leafIdx);    //Input: LeafIndex! TODO bräuchte hier (Inverse )  Local-to-Leaf Idx Map
+        row[k] = localView.index(rowLocal);
+        //         std::cout << "rowLocal:" << rowLocal << std::endl;
+        //         std::cout << "row[k]:" << row[k] << std::endl;
+      }
+    }
+    elementCounter++;
+  }
+  return row;
+}
+
+
+
+
+template<class LocalView, class Matrix, class localFunction1, class localFunction2>
+void computeElementStiffnessMatrix(const LocalView& localView,
+                                   Matrix& elementMatrix,
+                                   const localFunction1& mu,
+                                   const localFunction2& lambda,
+                                   const double gamma
+                                   )
+{
+  // Local StiffnessMatrix of the form:
+  // | phi*phi    m*phi |
+  // | phi *m     m*m   |
+  using Element = typename LocalView::Element;
+  const Element element = localView.element();
+  auto geometry = element.geometry();
+  constexpr int dim = Element::dimension;
+  constexpr int dimWorld = dim;
+  using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+  elementMatrix.setSize(localView.size()+3, localView.size()+3);         //extend by dim ´R_sym^{2x2}
+  elementMatrix = 0;
+
+  // LocalBasis-Offset
+  const int localPhiOffset = localView.size();
+
+  const auto& localFiniteElement = localView.tree().child(0).finiteElement();     
+  const auto nSf = localFiniteElement.localBasis().size();
+//   std::cout << "localView.size(): " << localView.size() << std::endl;
+//   std::cout << "nSf : " << nSf << std::endl;
+
+  ///////////////////////////////////////////////
+  // Basis for R_sym^{2x2}            // wird nicht als Funktion benötigt da konstant...
+  //////////////////////////////////////////////
+  MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+  MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+  MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+  std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+
+//   int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;  // TEST
+  int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+  const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+  
+  
+  
+//   std::cout << "Print QuadratureOrder:" << orderQR << std::endl;  //TEST
+
+  for (const auto& quadPoint : quad)
+  {
+    const auto& quadPos = quadPoint.position();
+    const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+    const auto integrationElement = geometry.integrationElement(quadPos);
+
+    std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
+    localFiniteElement.localBasis().evaluateJacobian(quadPos,
+                                                     referenceGradients);
+    // Compute the shape function gradients on the grid element
+    std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+
+    for (size_t i=0; i<gradients.size(); i++)
+      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+    for (size_t k=0; k < dimWorld; k++)
+      for (size_t i=0; i < nSf; i++)
+      {
+        // "phi*phi"-part
+        for (size_t l=0; l< dimWorld; l++)
+          for (size_t j=0; j < nSf; j++ )
+          {
+            // (scaled) Deformation gradient of the ansatz basis function
+            MatrixRT defGradientU(0);
+            defGradientU[k][0] = gradients[i][0];                       // Y
+            defGradientU[k][1] = gradients[i][1];                       //X2
+            defGradientU[k][2] = (1.0/gamma)*gradients[i][2];           //X3
+            // printmatrix(std::cout, defGradientU , "defGradientU", "--");
+            // (scaled)  Deformation gradient of the test basis function
+            MatrixRT defGradientV(0);
+            defGradientV[l][0] = gradients[j][0];                       // Y
+            defGradientV[l][1] = gradients[j][1];                       //X2
+            defGradientV[l][2] = (1.0/gamma)*gradients[j][2];           //X3
+              
+            double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+//             double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);  // also works..
+
+            size_t col = localView.tree().child(k).localIndex(i);                       
+            size_t row = localView.tree().child(l).localIndex(j);
+            elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
+            
+            // "m*phi"  & "phi*m" - part
+            for( size_t m=0; m<3; m++)
+            {
+                double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+
+                auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+                elementMatrix[row][localPhiOffset+m] += value;
+                elementMatrix[localPhiOffset+m][row] += value;
+            }
+          }
+      }
+    // "m*m"-part
+    for(size_t m=0; m<3; m++)
+      for(size_t n=0; n<3; n++)
+      {
+        double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
+        elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
+      }
+  }
+}
+
+
+
+// Get the occupation pattern of the stiffness matrix
+template<class Basis, class ParameterSet>
+void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet& parameterSet)
+{
+  //  OccupationPattern:
+  // | phi*phi    m*phi |
+  // | phi *m     m*m   |
+
+  auto localView = basis.localView();
+  const int phiOffset = basis.dimension();
+
+  nb.resize(basis.size()+3, basis.size()+3);
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    for (size_t i=0; i<localView.size(); i++)
+    {
+      // The global index of the i-th vertex of the element
+      auto row = localView.index(i);
+      for (size_t j=0; j<localView.size(); j++ )
+      {
+        // The global index of the j-th vertex of the element
+        auto col = localView.index(j);
+        nb.add(row[0],col[0]);                   // nun würde auch nb.add(row,col) gehen..
+      }
+    }
+    for (size_t i=0; i<localView.size(); i++)
+    {
+      auto row = localView.index(i);
+
+      for (size_t j=0; j<3; j++ )
+      {
+        nb.add(row,phiOffset+j);               // m*phi part of StiffnessMatrix
+        nb.add(phiOffset+j,row);               // phi*m part of StiffnessMatrix
+      }
+    }
+    for (size_t i=0; i<3; i++ )
+      for (size_t j=0; j<3; j++ )
+      {
+        nb.add(phiOffset+i,phiOffset+j);       // m*m part of StiffnessMatrix
+      }
+  }
+  //////////////////////////////////////////////////////////////////
+  // setOneBaseFunctionToZero
+  //////////////////////////////////////////////////////////////////
+  if(parameterSet.template get<bool>("set_oneBasisFunction_Zero ", true)){
+    FieldVector<int,3> row;
+    unsigned int arbitraryLeafIndex =  parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
+    unsigned int arbitraryElementNumber =  parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
+
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();    // macht keinen Unterschied ob hier k oder 0..
+    const auto nSf = localFiniteElement.localBasis().size();
+    assert(arbitraryLeafIndex < nSf );
+    assert(arbitraryElementNumber  < basis.gridView().size(0));   // "arbitraryElementNumber larger than total Number of Elements"
+
+    //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+    row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
+
+    for (int k = 0; k<3; k++)
+      nb.add(row[k],row[k]);
+  }
+}
+
+
+
+// Compute the source term for a single element
+// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
+template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void computeElementLoadVector( const LocalView& localView,
+                               LocalFunction1& mu,
+                               LocalFunction2& lambda,
+                               const double gamma,
+                               Vector& elementRhs,
+                               const Force& forceTerm
+                               )
+{
+  // | f*phi|
+  // | ---  |
+  // | f*m  |
+  using Element = typename LocalView::Element;
+  const auto element = localView.element();
+  const auto geometry = element.geometry();
+  constexpr int dim = Element::dimension;
+  constexpr int dimWorld = dim;
+
+  using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+  // Set of shape functions for a single element
+  const auto& localFiniteElement= localView.tree().child(0).finiteElement();
+  const auto nSf = localFiniteElement.localBasis().size();
+
+  elementRhs.resize(localView.size() +3);
+  elementRhs = 0;
+
+  // LocalBasis-Offset
+  const int localPhiOffset = localView.size();
+
+  ///////////////////////////////////////////////
+  // Basis for R_sym^{2x2}
+  //////////////////////////////////////////////
+  MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+  MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+  MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+  std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+
+//   int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;  // TEST
+  int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);   // für simplex
+  const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+  for (const auto& quadPoint : quad)
+  {
+    const FieldVector<double,dim>& quadPos = quadPoint.position();
+    const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+    const double integrationElement = geometry.integrationElement(quadPos);
+
+    std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+
+    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+    std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());  
+    for (size_t i=0; i< gradients.size(); i++)
+      jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+    // "f*phi"-part
+    for (size_t i=0; i < nSf; i++)
+      for (size_t k=0; k < dimWorld; k++)
+      {
+        // Deformation gradient of the ansatz basis function
+        MatrixRT defGradientV(0);
+        defGradientV[k][0] = gradients[i][0];                                 // Y
+        defGradientV[k][1] = gradients[i][1];                                 // X2
+        defGradientV[k][2] = (1.0/gamma)*gradients[i][2];                     // X3
+
+        double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientV, forceTerm(quadPos));
+        size_t row = localView.tree().child(k).localIndex(i);
+        elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
+      }
+
+    // "f*m"-part
+    for (size_t m=0; m<3; m++)
+    {
+      double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], forceTerm(quadPos));
+      elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
+    }
+  }
+}
+
+
+
+template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2, class ParameterSet>
+void assembleCellStiffness(const Basis& basis,
+                           LocalFunction1& muLocal,
+                           LocalFunction2& lambdaLocal,
+                           const double gamma,
+                           Matrix& matrix,
+                           ParameterSet& parameterSet
+                           )
+{
+  std::cout << "assemble Stiffness-Matrix begins." << std::endl;
+
+  MatrixIndexSet occupationPattern;
+  getOccupationPattern(basis, occupationPattern, parameterSet);
+  occupationPattern.exportIdx(matrix);
+  matrix = 0;
+
+  auto localView = basis.localView();
+  const int phiOffset = basis.dimension();
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    muLocal.bind(element);
+    lambdaLocal.bind(element);
+    const int localPhiOffset = localView.size();
+    
+    //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+    Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+    computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
+    //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+    //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
+    //std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
+
+    //////////////////////////////////////////////////////////////////////////////
+    // GLOBAL STIFFNES ASSEMBLY
+    //////////////////////////////////////////////////////////////////////////////
+    for (size_t i=0; i<localPhiOffset; i++)
+    for (size_t j=0; j<localPhiOffset; j++ )
+    {
+        auto row = localView.index(i);
+        auto col = localView.index(j);
+        matrix[row][col] += elementMatrix[i][j];
+    }
+    for (size_t i=0; i<localPhiOffset; i++)
+    for(size_t m=0; m<3; m++)
+    {
+        auto row = localView.index(i);
+        matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m];
+        matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
+    }
+    for (size_t m=0; m<3; m++ )
+    for (size_t n=0; n<3; n++ )
+        matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
+
+    //  printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
+  }
+}
+
+
+template<class Basis, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void assembleCellLoad(const Basis& basis,
+                      LocalFunction1& muLocal,
+                      LocalFunction2& lambdaLocal,
+                      const double gamma,
+                      Vector& b,
+                      const Force& forceTerm
+                      )
+{
+  //     std::cout << "assemble load vector." << std::endl;
+  b.resize(basis.size()+3);
+  b = 0;
+
+  auto localView = basis.localView();
+  const int phiOffset = basis.dimension();
+
+  // Transform G_alpha's to GridViewFunctions/LocalFunctions 
+  auto loadGVF  = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
+  auto loadFunctional = localFunction(loadGVF);      
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    muLocal.bind(element);
+    lambdaLocal.bind(element);
+    loadFunctional.bind(element);
+
+    const int localPhiOffset = localView.size();
+    //         std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+
+    BlockVector<FieldVector<double,1> > elementRhs;
+    computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
+//     printvector(std::cout, elementRhs, "elementRhs", "--");
+//     printvector(std::cout, elementRhs, "elementRhs", "--");
+    //////////////////////////////////////////////////////////////////////////////
+    // GLOBAL LOAD ASSEMBLY
+    //////////////////////////////////////////////////////////////////////////////
+    for (size_t p=0; p<localPhiOffset; p++)
+    {
+      auto row = localView.index(p);
+      b[row] += elementRhs[p];
+    }
+    for (size_t m=0; m<3; m++ )
+      b[phiOffset+m] += elementRhs[localPhiOffset+m];
+      //printvector(std::cout, b, "b", "--");
+  }
+}
+
+
+
+template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
+auto energy(const Basis& basis,
+            LocalFunction1& mu,
+            LocalFunction2& lambda,
+            const MatrixFunction& matrixFieldFuncA,
+            const MatrixFunction& matrixFieldFuncB)
+{
+
+//   TEST HIGHER PRECISION
+//   using float_50 = boost::multiprecision::cpp_dec_float_50;
+//   float_50 higherPrecEnergy = 0.0;
+  double energy = 0.0;
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  constexpr int dimWorld = dim;
+
+  auto localView = basis.localView();
+
+  auto matrixFieldAGVF  = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
+  auto matrixFieldA = localFunction(matrixFieldAGVF);
+  auto matrixFieldBGVF  = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
+  auto matrixFieldB = localFunction(matrixFieldBGVF);
+
+  using GridView = typename Basis::GridView;
+  using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+  using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+//   TEST
+//   FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
+//   printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
+
+  for (const auto& e : elements(basis.gridView()))
+  {
+    localView.bind(e);
+    matrixFieldA.bind(e);
+    matrixFieldB.bind(e);
+    mu.bind(e);
+    lambda.bind(e);
+
+    double elementEnergy = 0.0;
+    //double elementEnergy_HP = 0.0;
+
+    auto geometry = e.geometry();
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+//     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 );  // TEST
+    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+    const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+
+    for (const auto& quadPoint : quad) 
+    {
+      const auto& quadPos = quadPoint.position();
+      const double integrationElement = geometry.integrationElement(quadPos);
+
+      auto strain1 = matrixFieldA(quadPos);
+      auto strain2 = matrixFieldB(quadPos);
+
+      double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
+
+      elementEnergy += energyDensity * quadPoint.weight() * integrationElement;          
+      //elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
+    }
+    energy += elementEnergy;
+    //higherPrecEnergy += elementEnergy_HP;
+  }
+//   TEST
+//   std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
+  return energy;
+}
+
+
+
+template<class Basis, class Matrix, class Vector, class ParameterSet>
+void setOneBaseFunctionToZero(const Basis& basis,
+                              Matrix& stiffnessMatrix,
+                              Vector& load_alpha1,
+                              Vector& load_alpha2,
+                              Vector& load_alpha3,
+                              ParameterSet& parameterSet
+                              )
+{
+  std::cout << "Setting one Basis-function to zero" << std::endl;
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  
+  unsigned int arbitraryLeafIndex =  parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
+  unsigned int arbitraryElementNumber =  parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
+  //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+  FieldVector<int,3> row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
+
+  for (int k = 0; k < dim; k++)
+  {
+    load_alpha1[row[k]] = 0.0;
+    load_alpha2[row[k]] = 0.0;
+    load_alpha3[row[k]] = 0.0;
+    auto cIt    = stiffnessMatrix[row[k]].begin();
+    auto cEndIt = stiffnessMatrix[row[k]].end();
+    for (; cIt!=cEndIt; ++cIt)
+      *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
+  }
+}
+
+
+
+template<class Basis>
+auto childToIndexMap(const Basis& basis,
+                     const int k
+                     )
+{
+  // Input  : child/component
+  // Output : determine all Indices that belong to that component
+  auto localView = basis.localView();
+
+  std::vector<int> r = { };
+  //     for (int n: r)
+  //         std::cout << n << ","<< std::endl;
+
+  // Determine global indizes for each component k = 1,2,3.. in order to subtract correct (component of) integral Mean
+  // (global) Indices that correspond to component k = 1,2,3
+  for(const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    const auto& localFiniteElement = localView.tree().child(k).finiteElement();        
+    const auto nSf = localFiniteElement.localBasis().size();                           
+
+    for(size_t j=0; j<nSf; j++)
+    {
+      auto Localidx = localView.tree().child(k).localIndex(j);  // local  indices
+      r.push_back(localView.index(Localidx));                   // global indices
+    }
+  }
+  // Delete duplicate elements
+  // first remove consecutive (adjacent) duplicates
+  auto last = std::unique(r.begin(), r.end());
+  r.erase(last, r.end());
+  // sort followed by unique, to remove all duplicates
+  std::sort(r.begin(), r.end());
+  last = std::unique(r.begin(), r.end());
+  r.erase(last, r.end());
+
+  return r;
+}
+
+
+template<class Basis, class Vector>
+auto integralMean(const Basis& basis,
+                  Vector& coeffVector
+                  )
+{
+  // computes integral mean of given LocalFunction
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+
+  auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
+  auto localfun = localFunction(GVFunction);
+
+  auto localView = basis.localView();
+
+  FieldVector<double,3> r = {0.0, 0.0, 0.0};
+  double area = 0.0;
+
+  // Compute Area integral & integral of FE-function
+  for(const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    localfun.bind(element);
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+    
+//     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; //TEST
+    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+    const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+
+    for(const auto& quadPoint : quad)
+    {
+      const auto& quadPos = quadPoint.position();
+      const double integrationElement = element.geometry().integrationElement(quadPos);
+      area += quadPoint.weight() * integrationElement;
+      
+      r += localfun(quadPos) * quadPoint.weight() * integrationElement;
+    }
+  }
+  //     std::cout << "Domain-Area: " << area << std::endl;
+  return (1.0/area) * r ;
+}
+
+
+template<class Basis, class Vector>
+auto subtractIntegralMean(const Basis& basis,
+                          Vector& coeffVector
+                          )
+{
+  // Substract correct Integral mean from each associated component function
+  auto IM = integralMean(basis, coeffVector);
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+
+  for(size_t k=0; k<dim; k++)
+  {
+    //std::cout << "Integral-Mean: " << IM[k] << std::endl;
+    auto idx = childToIndexMap(basis,k);
+
+    for ( int i : idx)
+      coeffVector[i] -= IM[k];
+  }
+}
+
+
+
+//////////////////////////////////////////////////
+//   Infrastructure for handling periodicity
+//////////////////////////////////////////////////
+
+// Check whether two points are equal on R/Z x R/Z x R
+auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
+                {
+                    return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+                            and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
+                            and (FloatCmp::eq(x[2],y[2]))
+                        );
+                };
+
+
+                
+////////////////////////////////////////////////////////////// L2-ERROR /////////////////////////////////////////////////////////////////
+template<class Basis, class Vector, class MatrixFunction>
+double computeL2SymError(const Basis& basis,
+                         Vector& coeffVector,
+                         const double gamma,
+                         const MatrixFunction& matrixFieldFunc)
+{
+  double error = 0.0;
+
+
+  auto localView = basis.localView();
+  
+
+  constexpr int dim = Basis::LocalView::Element::dimension;             //TODO TEST 
+  constexpr int dimWorld = 3;                                           // Hier auch möglich? 
+  
+  auto matrixFieldGVF  = Dune::Functions::makeGridViewFunction(matrixFieldFunc, basis.gridView());
+  auto matrixField = localFunction(matrixFieldGVF);
+  using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    matrixField.bind(element);
+    auto geometry = element.geometry();
+
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();             
+    const auto nSf = localFiniteElement.localBasis().size();
+
+//     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; //TEST
+    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );  
+    const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+    for (const auto& quadPoint : quad)
+    {
+        const auto& quadPos = quadPoint.position();
+        const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+        const auto integrationElement = geometry.integrationElement(quadPos);
+
+        std::vector< FieldMatrix<double, 1, dim>> referenceGradients;
+        localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+
+        // Compute the shape function gradients on the grid element
+        std::vector<FieldVector<double,dim>> gradients(referenceGradients.size());
+
+        for (size_t i=0; i<gradients.size(); i++)
+          jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+        MatrixRT tmp(0);
+        double sum = 0.0;
+
+        for (size_t k=0; k < dimWorld; k++)
+        for (size_t i=0; i < nSf; i++)
+        {
+            size_t localIdx = localView.tree().child(k).localIndex(i);  // hier i:leafIdx
+            size_t globalIdx = localView.index(localIdx);
+
+            // (scaled) Deformation gradient of the ansatz basis function
+            MatrixRT defGradientU(0);
+            defGradientU[k][0] = coeffVector[globalIdx]*gradients[i][0];                       // Y  //hier i:leafIdx
+            defGradientU[k][1] = coeffVector[globalIdx]*gradients[i][1];                       //X2
+            defGradientU[k][2] = coeffVector[globalIdx]*(1.0/gamma)*gradients[i][2];           //X3
+
+            tmp += sym(defGradientU);
+        }
+//         printmatrix(std::cout, matrixField(quadPos), "matrixField(quadPos)", "--");
+//         printmatrix(std::cout, tmp, "tmp", "--");                                    // TEST Symphi_1 hat andere Struktur als analytic?
+//         tmp = tmp - matrixField(quadPos);
+//         printmatrix(std::cout, tmp - matrixField(quadPos), "Difference", "--");
+        for (int ii=0; ii<dimWorld; ii++)
+        for (int jj=0; jj<dimWorld; jj++)
+        {
+            sum +=  std::pow(tmp[ii][jj]-matrixField(quadPos)[ii][jj],2);
+        }
+//         std::cout << "sum:" << sum << std::endl;
+        error += sum * quadPoint.weight() * integrationElement;
+//         std::cout << "error:" << error << std::endl;
+    }
+  }
+  return sqrt(error);
+}
+////////////////////////////////////////////////////////////// L2-NORM /////////////////////////////////////////////////////////////////
+template<class Basis, class Vector>
+double computeL2Norm(const Basis& basis,
+                     Vector& coeffVector)
+{
+  // IMPLEMENTATION with makeDiscreteGlobalBasisFunction
+  double error = 0.0;
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  constexpr int dimWorld = dim;
+  auto localView = basis.localView();
+  auto GVFunc = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
+  auto localfun = localFunction(GVFunc);
+
+  for(const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    localfun.bind(element);
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+    
+//     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; //TEST
+    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+    const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+
+    for(const auto& quadPoint : quad)
+    {
+      const auto& quadPos = quadPoint.position();
+      const double integrationElement = element.geometry().integrationElement(quadPos);
+      error += localfun(quadPos)*localfun(quadPos) * quadPoint.weight() * integrationElement;
+    }
+  }
+  return sqrt(error);
+}
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[])
+{
+  MPIHelper::instance(argc, argv);
+
+  ParameterTree parameterSet;
+  if (argc < 2)
+    ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
+  else
+  {
+    ParameterTreeParser::readINITree(argv[1], parameterSet);
+    ParameterTreeParser::readOptions(argc, argv, parameterSet);
+  }
+
+  // Output setter
+//   std::string outputPath = parameterSet.get("outputPath", "../../outputs/output.txt");
+//   std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt");
+  std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs");
+//   std::string MatlabPath = parameterSet.get("MatlabPath", "/home/klaus/Desktop/DUNE/dune-microstructure/Matlab-Programs");
+//     std::string outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt";
+  std::fstream log;
+  log.open(outputPath + "/output.txt" ,std::ios::out);
+
+  std::cout << "outputPath:" << outputPath << std::endl;
+  
+//   parameterSet.report(log); // short Alternativ
+  
+  
+  constexpr int dim = 3;
+  constexpr int dimWorld = 3;
+
+  ///////////////////////////////////
+  // Get Parameters/Data
+  ///////////////////////////////////
+  double gamma = parameterSet.get<double>("gamma",1.0);   // ratio dimension reduction to homogenization
+  double alpha = parameterSet.get<double>("alpha", 2.0);
+  double theta = parameterSet.get<double>("theta",1.0/4.0); 
+  ///////////////////////////////////
+  // Get Material Parameters
+  ///////////////////////////////////
+  std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
+  log << "material_prestrain used: "<< imp  << std::endl;
+  double beta = parameterSet.get<double>("beta",2.0); 
+  double mu1 = parameterSet.get<double>("mu1",10.0);;
+  double mu2 = beta*mu1;
+  double lambda1 = parameterSet.get<double>("lambda1",0.0);;
+  double lambda2 = beta*lambda1;
+  
+  // Plate geometry settings
+  double width = parameterSet.get<double>("width", 1.0);   //geometry cell, cross section
+  //     double len  = parameterSet.get<double>("len", 10.0); //length
+  //     double height  = parameterSet.get<double>("h", 0.1); //rod thickness
+  //     double eps  = parameterSet.get<double>("eps", 0.1); //size of perioticity cell
+
+  ///////////////////////////////////
+  // Get Prestrain/Parameters
+  ///////////////////////////////////
+//   unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", "analytical_Example");  //OLD 
+//   std::string prestraintype = parameterSet.get<std::string>("prestrainType", "analytical_Example");
+//   double rho1 = parameterSet.get<double>("rho1", 1.0);
+//   double rho2 = alpha*rho1;
+//   auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
+//   auto B_Term = prestrainImp.getPrestrain(prestraintype);
+  
+  
+  
+  auto prestrainImp = PrestrainImp<dim>(); //NEW 
+  auto B_Term = prestrainImp.getPrestrain(parameterSet);
+
+  log << "----- Input Parameters -----: " << std::endl;
+//   log << "prestrain imp: " <<  prestraintype << "\nrho1 = " << rho1 << "\nrho2 = " << rho2  << std::endl;
+  log << "alpha: " << alpha << std::endl;
+  log << "gamma: " << gamma << std::endl;
+  log << "theta: " << theta << std::endl;
+  log << "beta: " << beta << std::endl;
+  log << "material parameters: " << std::endl;
+  log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+  log << "lambda1: " << lambda1 <<"\nlambda2: " << lambda2 << std::endl;
+  log << "----------------------------: " << std::endl;
+
+  ///////////////////////////////////
+  // Generate the grid
+  ///////////////////////////////////
+
+  //Corrector Problem Domain:
+  unsigned int cellDomain = parameterSet.get<unsigned int>("cellDomain", 1);
+  // (shifted) Domain (-1/2,1/2)^3
+  FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
+  FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
+  if (cellDomain == 2)
+  {
+    // Domain : [0,1)^2 x (-1/2, 1/2)
+    FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
+    FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
+  }
+
+  std::array<int,2> numLevels = parameterSet.get<std::array<int,2>>("numLevels", {1,3});
+
+  int levelCounter = 0;
+
+  ///////////////////////////////////
+  // Create Data Storage
+  ///////////////////////////////////
+  // Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4  L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
+  std::vector<std::variant<std::string, size_t , double>> Storage_Error;
+  
+  // Storage:: #1 level #2 |q1_a-q1_c| #3 |q2_a-q2_c| #4 |q3_a-q3_c| #5 |b1_a-b1_c| #6 |b2_a-b2_c| #7 |b3_a-b3_c|           
+   std::vector<std::variant<std::string, size_t , double>> Storage_Quantities;         
+
+
+//   for(const size_t &level : numLevels)                             // explixite Angabe der levels.. {2,4}
+  for(size_t level = numLevels[0] ; level <= numLevels[1]; level++)    // levels von bis.. [2,4]
+  {
+    std::cout << " ----------------------------------" << std::endl;
+    std::cout << "Level: " << level << std::endl;
+    std::cout << " ----------------------------------" << std::endl;
+
+    Storage_Error.push_back(level);
+    Storage_Quantities.push_back(level);
+    std::array<int, dim> nElements = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
+
+    std::cout << "Number of Elements in each direction: " << nElements << std::endl;
+    log << "Number of Elements in each direction: " << nElements << std::endl;
+
+    using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+    CellGridType grid_CE(lower,upper,nElements);
+    using GridView = CellGridType::LeafGridView;
+    const GridView gridView_CE = grid_CE.leafGridView();
+
+    using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+    using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+    using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+    using VectorCT = BlockVector<FieldVector<double,1> >;
+    using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
+
+    ///////////////////////////////////
+    //  Create Lambda-Functions for material Parameters depending on microstructure
+    ///////////////////////////////////
+    
+    auto materialImp = IsotropicMaterialImp<dim>();
+    auto muTerm = materialImp.getMu(parameterSet);
+    auto lambdaTerm = materialImp.getLambda(parameterSet);
+
+/*  
+    auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+                    if (abs(z[0]) >= (theta/2.0))                                                   
+                        return mu1;
+                    else
+                        return mu2;
+                    };
+
+    auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
+                    if (abs(z[0]) >= (theta/2.0))
+                        return lambda1;
+                    else
+                        return lambda2;
+                    };*/
+
+    auto muGridF  = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
+    auto muLocal = localFunction(muGridF);
+    auto lambdaGridF  = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
+    auto lambdaLocal = localFunction(lambdaGridF);
+
+    ///////////////////////////////////
+    // --- Choose Solver ---
+    // 1 : CG-Solver
+    // 2 : GMRES
+    // 3 : QR
+    ///////////////////////////////////
+    unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+    
+    // Print Options 
+    bool print_debug = parameterSet.get<bool>("print_debug", false);
+    
+    
+    bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
+    bool write_prestrainFunctions  = parameterSet.get<bool>("write_prestrainFunctions", false);
+    
+    
+    bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
+    bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
+    bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
+    bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
+    bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
+
+    /////////////////////////////////////////////////////////
+    // Choose a finite element space for Cell Problem
+    /////////////////////////////////////////////////////////
+    using namespace Functions::BasisFactory;
+    Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+
+    // Don't do the following in real life: It has quadratic run-time in the number of vertices.
+    for (const auto& v1 : vertices(gridView_CE))
+        for (const auto& v2 : vertices(gridView_CE))
+            if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
+            {
+                periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
+            }
+
+    // First order periodic Lagrange-Basis
+    auto Basis_CE = makeBasis(
+        gridView_CE,
+        power<dim>(                                                                             // eig dimworld?!?! 
+        Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
+        flatLexicographic()));     
+
+    log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
+    const int phiOffset = Basis_CE.size();
+
+    /////////////////////////////////////////////////////////
+    // Data structures: Stiffness matrix and right hand side vectors
+    /////////////////////////////////////////////////////////
+    VectorCT load_alpha1, load_alpha2, load_alpha3;
+    MatrixCT stiffnessMatrix_CE;
+
+    bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
+    bool set_oneBasisFunction_Zero = false;
+    bool substract_integralMean = false;
+    if(set_IntegralZero)
+    {
+        set_oneBasisFunction_Zero = true;
+        substract_integralMean = true;
+    }
+
+    /////////////////////////////////////////////////////////
+    // Define "rhs"-functions
+    /////////////////////////////////////////////////////////
+    Func2Tensor x3G_1 = [] (const Domain& x) {
+                            return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};    //TODO könnte hier sign übergeben?
+                        };
+
+    Func2Tensor x3G_2 = [] (const Domain& x) {
+                            return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
+                        };
+
+    Func2Tensor x3G_3 = [] (const Domain& x) {                                                                               
+                            return MatrixRT{{0.0, 1.0/sqrt(2.0)*x[2], 0.0}, {1.0/sqrt(2.0)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+                        };
+                        
+    Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {return -1.0*x3G_1(x);};
+    Func2Tensor x3G_2neg = [x3G_2] (const Domain& x) {return -1.0*x3G_2(x);};
+    Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {return -1.0*x3G_3(x);};
+    
+    
+    //TEST 
+//     auto QuadraticForm = [] (const double mu, const double lambda, const MatrixRT& M) {
+// 
+//                                 return lambda*std::pow(trace(M),2) + 2*mu*pow(norm(sym(M)),2);
+//                         };
+
+
+    // TEST - energy method ///
+    // different indicatorFunction in muTerm has impact here !!
+    //     double GGterm = 0.0;
+    //     GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3G_1, x3G_1  );   // <L i(x3G_alpha) , i(x3G_beta) >
+    //     std::cout << "GGTerm:" << GGterm << std::endl;
+    //     std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+
+
+    ///////////////////////////////////////////////
+    // Basis for R_sym^{2x2}
+    //////////////////////////////////////////////
+    MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+    MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+    MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+    std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+    //  printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+    //  printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+    //  printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
+    //  log <<  basisContainer[0] << std::endl;
+
+
+    //////////////////////////////////////////////////////////////////////
+    // Determine global indices of components for arbitrary (local) index
+    //////////////////////////////////////////////////////////////////////
+    int arbitraryLeafIndex =  parameterSet.get<unsigned int>("arbitraryLeafIndex", 0);            // localIdx..assert Number < #ShapeFcts 
+    int arbitraryElementNumber =  parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
+
+    FieldVector<int,3> row;
+     row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
+    if(print_debug)
+      printvector(std::cout, row, "row" , "--");
+
+    /////////////////////////////////////////////////////////
+    // Assemble the system
+    /////////////////////////////////////////////////////////
+    assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma,  stiffnessMatrix_CE, parameterSet);
+    assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1neg);
+    assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2neg);
+    assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg);
+    //   printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
+    //   printvector(std::cout, load_alpha1, "load_alpha1", "--");
+
+
+
+    // set one basis-function to zero
+    if(set_oneBasisFunction_Zero)
+    {
+        setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3, parameterSet);
+    //     printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
+    //     printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
+    }
+
+    ////////////////////////////////////////////////////
+    // Compute solution
+    ////////////////////////////////////////////////////
+    VectorCT x_1 = load_alpha1;
+    VectorCT x_2 = load_alpha2;
+    VectorCT x_3 = load_alpha3;
+
+//     auto load_alpha1BS = load_alpha1;
+    //   printvector(std::cout, load_alpha1, "load_alpha1 before SOLVER", "--" );
+    //   printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );
+
+    if (Solvertype == 1)  // CG - SOLVER
+    {
+        std::cout << "------------ CG - Solver ------------" << std::endl;
+        MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
+
+
+
+        // Sequential incomplete LU decomposition as the preconditioner
+        SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_CE,1.0);
+        int iter = parameterSet.get<double>("iterations_CG", 999);
+        // Preconditioned conjugate-gradient solver
+        CGSolver<VectorCT> solver(op,
+                                ilu0, //NULL,
+                                1e-8, // desired residual reduction factorlack
+                                iter,   // maximum number of iterations
+                                2);   // verbosity of the solver
+        InverseOperatorResult statistics;
+        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";
+        solver.apply(x_2, load_alpha2, statistics);
+        std::cout << "solve linear system for x_3.\n";
+        solver.apply(x_3, load_alpha3, statistics);
+        log << "Solver-type used: " <<" CG-Solver" << std::endl;
+    }
+    ////////////////////////////////////////////////////////////////////////////////////
+
+    else if (Solvertype ==2)  // GMRES - SOLVER
+    {
+
+        std::cout << "------------ GMRES - Solver ------------" << std::endl;
+        // Turn the matrix into a linear operator
+        MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
+
+        // Fancy (but only) way to not have a preconditioner at all
+        Richardson<VectorCT,VectorCT> preconditioner(1.0);
+
+        // Construct the iterative solver
+        RestartedGMResSolver<VectorCT> solver(
+            stiffnessOperator, // Operator to invert
+            preconditioner,    // Preconditioner
+            1e-10,             // Desired residual reduction factor
+            500,               // Number of iterations between restarts,
+                            // here: no restarting
+            500,               // Maximum number of iterations
+            2);                // Verbosity of the solver
+
+        // Object storing some statistics about the solving process
+        InverseOperatorResult statistics;
+
+        // solve for different Correctors (alpha = 1,2,3)
+        solver.apply(x_1, load_alpha1, statistics);
+        solver.apply(x_2, load_alpha2, statistics);
+        solver.apply(x_3, load_alpha3, statistics);
+        log << "Solver-type used: " <<" GMRES-Solver" << std::endl;
+    }
+    ////////////////////////////////////////////////////////////////////////////////////
+    else if ( Solvertype == 3)// QR - SOLVER
+    {
+
+        std::cout << "------------ QR - Solver ------------" << std::endl;
+        log << "solveLinearSystems: We use QR solver.\n";
+        //TODO install suitesparse
+        SPQR<MatrixCT> sPQR(stiffnessMatrix_CE);
+        sPQR.setVerbosity(1);
+        InverseOperatorResult statisticsQR;
+
+        sPQR.apply(x_1, load_alpha1, statisticsQR);
+        sPQR.apply(x_2, load_alpha2, statisticsQR);
+        sPQR.apply(x_3, load_alpha3, statisticsQR);
+        log << "Solver-type used: " <<" QR-Solver" << std::endl;
+    }
+    //     printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
+    //     printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
+    //     printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
+
+    ////////////////////////////////////////////////////////////////////////////////////
+    // Extract phi_alpha  &  M_alpha coefficients
+    ////////////////////////////////////////////////////////////////////////////////////
+    VectorCT phi_1, phi_2, phi_3;
+    phi_1.resize(Basis_CE.size());
+    phi_1 = 0;
+    phi_2.resize(Basis_CE.size());
+    phi_2 = 0;
+    phi_3.resize(Basis_CE.size());
+    phi_3 = 0;
+
+    for(size_t i=0; i<Basis_CE.size(); i++)
+    {
+        phi_1[i] = x_1[i];
+        phi_2[i] = x_2[i];
+        phi_3[i] = x_3[i];
+    }
+
+    FieldVector<double,3> m_1, m_2, m_3;
+
+    for(size_t i=0; i<3; i++)
+    {
+        m_1[i] = x_1[phiOffset+i];
+        m_2[i] = x_2[phiOffset+i];
+        m_3[i] = x_3[phiOffset+i];
+    }
+    // assemble M_alpha's (actually iota(M_alpha) )
+
+    MatrixRT M_1(0), M_2(0), M_3(0);
+
+    for(size_t i=0; i<3; i++)
+    {
+        M_1 += m_1[i]*basisContainer[i];
+        M_2 += m_2[i]*basisContainer[i];
+        M_3 += m_3[i]*basisContainer[i];
+    }
+
+    std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
+    printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
+    printmatrix(std::cout, M_2, "Corrector-Matrix M_2", "--");
+    printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
+    log << "---------- OUTPUT ----------" << std::endl;
+    log << " --------------------" << std::endl;
+    log << "Corrector-Matrix M_1: \n" << M_1 << std::endl;
+    log << " --------------------" << std::endl;
+    log << "Corrector-Matrix M_2: \n" << M_2 << std::endl;
+    log << " --------------------" << std::endl;
+    log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
+    log << " --------------------" << std::endl;
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Substract Integral-mean
+    ////////////////////////////////////////////////////////////////////////////                                                                         
+    using SolutionRange = FieldVector<double,dim>;
+
+    if(write_IntegralMean)
+    {
+        std::cout << "check integralMean phi_1: " << std::endl;
+        auto A = integralMean(Basis_CE,phi_1);
+        for(size_t i=0; i<3; i++)
+        {
+            std::cout << "Integral-Mean : " << A[i] << std::endl;
+        }
+    }
+    if(substract_integralMean)
+    {
+        std::cout << " --- substracting integralMean --- " << std::endl;
+        subtractIntegralMean(Basis_CE,  phi_1);
+        subtractIntegralMean(Basis_CE,  phi_2);
+        subtractIntegralMean(Basis_CE,  phi_3);
+        subtractIntegralMean(Basis_CE,  x_1);
+        subtractIntegralMean(Basis_CE,  x_2);
+        subtractIntegralMean(Basis_CE,  x_3);
+        //////////////////////////////////////////
+        // CHECK INTEGRAL-MEAN:
+        //////////////////////////////////////////
+        if(write_IntegralMean)
+        {
+            auto A = integralMean(Basis_CE, phi_1);
+            for(size_t i=0; i<3; i++)
+            {
+            std::cout << "Integral-Mean (CHECK)  : " << A[i] << std::endl;
+            }
+        }
+    }
+
+    /////////////////////////////////////////////////////////
+    // Write Solution (Corrector Coefficients) in Logs
+    /////////////////////////////////////////////////////////
+//     log << "\nSolution of Corrector problems:\n";
+    if(write_corrector_phi1)
+    {
+        log << "\nSolution of Corrector problems:\n";
+        log << "\n Corrector_phi1:\n";
+        log << x_1 << std::endl;
+    }
+    if(write_corrector_phi2)
+    {
+        log << "-----------------------------------------------------";
+        log << "\n Corrector_phi2:\n";
+        log << x_2 << std::endl;
+    }
+    if(write_corrector_phi3)
+    {
+        log << "-----------------------------------------------------";
+        log << "\n Corrector_phi3:\n";
+        log << x_3 << std::endl;
+    }
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Make a discrete function from the FE basis and the coefficient vector
+    ////////////////////////////////////////////////////////////////////////////                                                                                     // ERROR
+    auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+        Basis_CE,
+        phi_1);
+
+    auto correctorFunction_2 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+        Basis_CE,
+        phi_2);
+
+    auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+        Basis_CE,
+        phi_3);
+
+    /////////////////////////////////////////////////////////
+    // Create Containers for Basis-Matrices, Correctors and Coefficients
+    /////////////////////////////////////////////////////////
+    const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
+    const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
+    const std::array<VectorCT, 3> loadContainer = {load_alpha1, load_alpha2, load_alpha3};
+    const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
+
+
+    /////////////////////////////////////////////////////////
+    // Compute effective quantities: Elastic law & Prestrain
+    /////////////////////////////////////////////////////////
+    MatrixRT Q(0);
+    VectorCT tmp1,tmp2;
+    FieldVector<double,3> B_hat ;
+
+    // Compute effective elastic law Q
+    for(size_t a = 0; a < 3; a++)
+        for(size_t b=0; b < 3; b++)
+        {
+            assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]);   // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
+
+            double GGterm = 0.0;
+            GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b]  );   // <L i(x3G_alpha) , i(x3G_beta) >
+
+            // TEST
+            // std::setprecision(std::numeric_limits<float>::digits10);
+
+            Q[a][b] =  (coeffContainer[a]*tmp1) + GGterm;                         // seems symmetric...check positiv definitness?
+        
+            if (print_debug)
+            {
+                std::cout << "analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+                std::cout << "GGTerm:" << GGterm << std::endl;
+                std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
+            }
+        }
+    printmatrix(std::cout, Q, "Matrix Q", "--");
+    log << "Effective Matrix Q: " << std::endl;
+    log << Q << std::endl;
+
+    // compute B_hat
+    for(size_t a = 0; a<3; a++)
+    {
+        assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B_Term);  // <L i(M_alpha) + sym(grad phi_alpha), B >
+        auto GBterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
+        B_hat[a] = (coeffContainer[a]*tmp2) + GBterm;
+        
+        if (print_debug)
+        {
+            std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl;  //see orthotropic.tex
+        }
+    }
+
+//     log << B_hat << std::endl;
+//     log << "Prestrain B_hat: " << std::endl;
+//     log << B_hat << std::endl;
+
+    std::cout << "check this Contribution: " << (coeffContainer[2]*tmp2) << std::endl;  //see orthotropic.tex
+
+    ///////////////////////////////
+    // Compute effective Prestrain B_eff
+    //////////////////////////////
+    FieldVector<double, 3> Beff;
+    Q.solve(Beff,B_hat);
+    
+    log << "--- Prestrain Output --- " << std::endl;
+    log << "B_hat: " << B_hat << std::endl;
+    log << "B_eff: " << Beff <<  " (Effective Prestrain)" << std::endl;
+    log << "------------------------ " << std::endl;
+//     log << Beff << std::endl;
+//     log << "Effective Prestrain B_eff: " << std::endl;
+//     log << Beff << std::endl;
+//     printvector(std::cout, Beff, "Beff", "--");
+
+    auto q1 = Q[0][0];
+    auto q2 = Q[1][1];
+    auto q3 = Q[2][2];
+    
+    std::cout<< "q1 : " << q1 << std::endl;
+    std::cout<< "q2 : " << q2 << std::endl;
+    std::cout.precision(10);
+    std::cout << "q3 : " << std::fixed << q3 << std::endl;
+//     std::cout<< "q3 : " << q3 << std::endl;
+    printvector(std::cout, B_hat, "B_hat", "--");
+    printvector(std::cout, Beff, "Beff", "--");
+    
+    std::cout << "Theta: " << theta << std::endl;
+    std::cout << "Gamma: " << gamma << std::endl;
+    std::cout << "Number of Elements in each direction: " << nElements << std::endl;
+    log << "q1=" << q1 << std::endl;
+    log << "q2=" << q2 << std::endl;
+    log << "q3=" << q3 << std::endl;
+    log << "q12=" << Q[0][1] << std::endl;
+    
+    log << "b1=" << Beff[0] << std::endl;
+    log << "b2=" << Beff[1] << std::endl;
+    log << "b3=" << Beff[2] << std::endl;
+    log << "b1_hat: " << B_hat[0] << std::endl;
+    log << "b2_hat: " << B_hat[1] << std::endl;
+    log << "b3_hat: " << B_hat[2] << std::endl;
+    log << "mu_gamma=" << q3 << std::endl;           // added for Python-Script
+ 
+    log << std::fixed << std::setprecision(6) << "q_onetwo=" << Q[0][1] << std::endl;
+//     log << "q_onetwo=" << Q[0][1] << std::endl;           // added for Python-Script
+
+    //////////////////////////////////////////////////////////////
+    // Define Analytic Solutions
+    //////////////////////////////////////////////////////////////
+    // double b1 = (mu1*p1/4)*(beta/(theta+(1-theta)*beta))*(1-theta*(1+alpha));
+    // double b2 = (mu1*p1/8)*(1-theta*(1+beta*alpha));
+    double b1_hat_ana = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
+    double b2_hat_ana = -(theta/4.0)*mu2;
+    double b3_hat_ana = 0.0;
+   
+    double b1_eff_ana = (-3.0/2.0)*theta;
+    double b2_eff_ana = (-3.0/2.0)*(theta*mu2)/(mu1*(1-theta)+mu2*theta);
+    double b3_eff_ana = 0.0;
+    
+//     double q1_ana = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+//     double q2_ana = ((1.0-theta)*mu1+theta*mu2)/6.0;
+    double q1_ana = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);  // 1/6 * harmonicMean
+    double q2_ana = mu1*((1-theta)+theta*beta)        *(1.0/6.0);     // 1/6 * arithmeticMean
+    
+    
+    
+    
+    
+    std::cout << "Test B_hat & B_eff" << std::endl;
+     double p1 = parameterSet.get<double>("rho1", 1.0);
+     double alpha = parameterSet.get<double>("alpha", 2.0);
+     double p2 = alpha*p1;
+    
+    std::cout << ((3.0*p1)/2.0)*beta*(1-(theta*(1+alpha)))   << std::endl;  // TODO ERROR in paper ?? 
+
+    std::cout << "----- print analytic solutions -----" << std::endl;
+    std::cout << "b1_hat_ana : " << b1_hat_ana << std::endl;
+    std::cout << "b2_hat_ana : " << b2_hat_ana << std::endl;
+    std::cout << "b3_hat_ana : " << b3_hat_ana << std::endl;
+    std::cout << "b1_eff_ana : " << b1_eff_ana << std::endl;
+    std::cout << "b2_eff_ana : " << b2_eff_ana << std::endl;
+    std::cout << "b3_eff_ana : " << b3_eff_ana << std::endl;
+    
+    std::cout << "q1_ana : "     << q1_ana << std::endl;
+    std::cout << "q2_ana : "     << q2_ana << std::endl;
+    std::cout << "q3 should be between q1 and q2"  << std::endl;
+    log << "----- print analytic solutions -----" << std::endl;
+    log << "b1_hat_ana : " << b1_hat_ana << std::endl;
+    log << "b2_hat_ana : " << b2_hat_ana << std::endl;
+    log << "b3_hat_ana : " << b3_hat_ana << std::endl;
+    log << "b1_eff_ana : " << b1_eff_ana << std::endl;
+    log << "b2_eff_ana : " << b2_eff_ana << std::endl;
+    log << "b3_eff_ana : " << b3_eff_ana << std::endl;
+    log << "q1_ana : "     << q1_ana << std::endl;
+    log << "q2_ana : "     << q2_ana << std::endl;
+    log << "q3 should be between q1 and q2"  << std::endl;
+    
+    if (imp == "analytical_Example")   // print Errors only for analytical_Example
+    {
+        Storage_Quantities.push_back(std::abs(q1_ana - q1));
+        Storage_Quantities.push_back(std::abs(q2_ana - q2));
+        Storage_Quantities.push_back(q3);
+        Storage_Quantities.push_back(std::abs(b1_eff_ana - Beff[0]));
+        Storage_Quantities.push_back(std::abs(b2_eff_ana - Beff[1]));
+        Storage_Quantities.push_back(std::abs(b3_eff_ana - Beff[2]));
+    }
+    else
+    {
+        Storage_Quantities.push_back(q1);
+        Storage_Quantities.push_back(q2);
+        Storage_Quantities.push_back(q3);
+        Storage_Quantities.push_back(Beff[0]);
+        Storage_Quantities.push_back(Beff[1]);
+        Storage_Quantities.push_back(Beff[2]);
+    }
+
+
+    auto symPhi_1_analytic = [mu1, mu2, theta, muTerm] (const Domain& x) {
+                        return MatrixRT{{  (((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x))) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+                    };
+
+    auto zeroFunction = [] (const Domain& x) {
+                    return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+                    };
+
+    if(write_L2Error)
+    {
+//         std::cout << " ----- L2ErrorSym ----" << std::endl;
+        auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+//         std::cout << "L2SymError: " << L2SymError << std::endl;
+//         std::cout << " -----------------" << std::endl;
+
+//         std::cout << " ----- L2NORMSym ----" << std::endl;
+        auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction);                           
+//         std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl;           
+        VectorCT zeroVec;
+        zeroVec.resize(Basis_CE.size());
+        zeroVec = 0;
+        auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma, symPhi_1_analytic);
+//         std::cout << "L2-Norm(SymAnalytic): " << L2Norm_SymAnalytic << std::endl;
+//         std::cout << " -----------------" << std::endl;
+
+//         std::cout << " ----- L2NORM ----" << std::endl;
+        auto L2Norm = computeL2Norm(Basis_CE,phi_1);
+//         std::cout << "L2Norm(phi_1): "  << L2Norm << std::endl;
+//         std::cout << " -----------------" << std::endl;
+        
+        
+        
+//         log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
+//         log << "L2-Norm(Symphi_1): "    << L2Norm_Symphi<< std::endl;
+//         log << "L2-Norm(SymAnalytic): " << L2Norm_SymAnalytic << std::endl;
+
+        double EOC = 0.0;
+        Storage_Error.push_back(L2SymError);
+
+        //Compute Experimental order of convergence (EOC)
+        if(levelCounter > 0)
+        {
+            // Storage_Error:: #1 level #2 L2SymError #3 L2SymErrorOrder #4  L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
+//             std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter-1)*6)+1 << std::endl;           // Besser std::map ???
+//             std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl;             // für Storage_Error[idx] muss idx zur compile time feststehen?!
+            auto ErrorOld = std::get<double>(Storage_Error[((levelCounter-1)*6)+1]);
+            auto ErrorNew = std::get<double>(Storage_Error[(levelCounter*6)+1]);
+//
+            EOC = std::log(ErrorNew/ErrorOld)/(-1*std::log(2));
+            //   std::cout << "Storage_Error[0] :" << std::get<1>(Storage_Error[0]) << std::endl;
+            //   std::cout << "output variant :" << std::get<std::string>(Storage_Error[1]) << std::endl;
+            //   std::cout << "output variant :" << std::get<0>(Storage_Error[1]) << std::endl;
+        }
+//         Storage_Error.push_back(level);
+        Storage_Error.push_back(EOC);
+        Storage_Error.push_back(L2Norm_Symphi);
+        Storage_Error.push_back(L2Norm_SymAnalytic);
+        Storage_Error.push_back(L2Norm);
+    }
+  //////////////////////////////////////////////////////////////////////////////////////////////
+  // Write Data to Matlab / Optimization-Code
+  //////////////////////////////////////////////////////////////////////////////////////////////
+//   writeMatrixToMatlab(Q, "../../Matlab-Programs/QMatrix.txt");
+  writeMatrixToMatlab(Q, outputPath + "/QMatrix.txt");
+  
+  // write effective Prestrain in Matrix for Output
+  FieldMatrix<double,1,3> BeffMat;
+  BeffMat[0] = Beff;
+//   writeMatrixToMatlab(BeffMat, "../../Matlab-Programs/BMatrix.txt");
+  writeMatrixToMatlab(BeffMat, outputPath + "/BMatrix.txt");
+  
+  //////////////////////////////////////////////////////////////////////////////////////////////
+  // Write result to VTK file
+  //////////////////////////////////////////////////////////////////////////////////////////////
+  std::string vtkOutputName = outputPath + "/CellProblem-result";
+  
+  VTKWriter<GridView> vtkWriter(gridView_CE);
+
+  vtkWriter.addVertexData(
+    correctorFunction_1,
+    VTK::FieldInfo("Corrector phi_1 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));     
+  vtkWriter.addVertexData(
+    correctorFunction_2,
+    VTK::FieldInfo("Corrector phi_2 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
+  vtkWriter.addVertexData(
+    correctorFunction_3,
+    VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
+  vtkWriter.write( vtkOutputName  + "-level"+ std::to_string(level));
+//     vtkWriter.pwrite( vtkOutputName  + "-level"+ std::to_string(level), outputPath, "");   // TEST Write to folder "/outputs" 
+//   vtkWriter.pwrite( vtkOutputName  + "-level"+ std::to_string(level), outputPath, "", VTK::OutputType::ascii, 0, 0 );
+  std::cout << "wrote data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;      
+
+  
+  
+   if (write_materialFunctions)
+   {
+        using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+        VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
+        using GridViewVTK = VTKGridType::LeafGridView;
+        const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+        
+        auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+        auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+        std::vector<double> mu_CoeffP0;
+        Functions::interpolate(scalarP0FeBasis, mu_CoeffP0, muTerm);
+        auto mu_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, mu_CoeffP0);
+        
+        std::vector<double> mu_CoeffP1;
+        Functions::interpolate(scalarP1FeBasis, mu_CoeffP1, muTerm);
+        auto mu_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, mu_CoeffP1);
+        
+        
+        std::vector<double> lambda_CoeffP0;
+        Functions::interpolate(scalarP0FeBasis, lambda_CoeffP0, lambdaTerm);
+        auto lambda_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, lambda_CoeffP0);
+        
+        std::vector<double> lambda_CoeffP1;
+        Functions::interpolate(scalarP1FeBasis, lambda_CoeffP1, lambdaTerm);
+        auto lambda_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, lambda_CoeffP1);
+        
+        VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+        
+        MaterialVtkWriter.addVertexData(
+            mu_DGBF_P1,
+            VTK::FieldInfo("mu_P1", VTK::FieldInfo::Type::scalar, 1));    
+        MaterialVtkWriter.addCellData(
+            mu_DGBF_P0,
+            VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));    
+        MaterialVtkWriter.addVertexData(
+            lambda_DGBF_P1,
+            VTK::FieldInfo("lambda_P1", VTK::FieldInfo::Type::scalar, 1));    
+        MaterialVtkWriter.addCellData(
+            lambda_DGBF_P0,
+            VTK::FieldInfo("lambda_P0", VTK::FieldInfo::Type::scalar, 1));    
+        
+        MaterialVtkWriter.write(outputPath + "/MaterialFunctions-level"+ std::to_string(level) );
+        std::cout << "wrote data to file:" + outputPath +"/MaterialFunctions-level" + std::to_string(level) << std::endl;  
+        
+   }
+  
+   if (write_prestrainFunctions)
+   {
+    using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+    VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
+    using GridViewVTK = VTKGridType::LeafGridView;
+    const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+   
+    FTKfillerContainer<dim> VTKFiller;
+    VTKFiller.vtkPrestrainNorm(gridView_VTK, B_Term, "PrestrainBNorm");
+    
+    // WORKS Too 
+    VTKFiller.vtkProblemCell(gridView_VTK, B_Term, muLocal,"VTKProblemCell");;
+   }
+  
+  
+  levelCounter++; 
+  } // Level-Loop End
+  
+    
+  
+
+    /////////////////////////////////////////
+    // Print Storage
+    /////////////////////////////////////////
+    int tableWidth = 12;
+    if (imp == "analytical_Example")   // print Errors only for analytical_Example
+    {
+        std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+        std::cout << center("Levels",tableWidth)       << " | "
+                << center("L2SymError",tableWidth)     << " | "
+                << center("Order",tableWidth)     << " | "
+                << center("L2SymNorm",tableWidth)     << " | "
+                << center("L2SymNorm_ana",tableWidth)     << " | "
+                << center("L2Norm",tableWidth)  << " | " << "\n";
+        std::cout << std::string(tableWidth*6 + 3*6, '-') << "\n";
+        log     << std::string(tableWidth*6 + 3*6, '-') << "\n";
+        log     << center("Levels",tableWidth)       << " | "
+                << center("L2SymError",tableWidth)     << " | "
+                << center("Order",tableWidth)     << " | "
+                << center("L2SymNorm",tableWidth)     << " | "
+                << center("L2SNorm_ana",tableWidth)     << " | "
+                << center("L2Norm",tableWidth)  << " | " << "\n";
+        log     << std::string(tableWidth*6 + 3*6, '-') << "\n";
+        
+        int StorageCount = 0;
+
+        for(auto& v: Storage_Error) 
+        {
+            std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth)      << " | ";}, v);      // Anzahl-Nachkommastellen
+            std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth)      << " & ";}, v);
+            StorageCount++;
+            if(StorageCount % 6 == 0 )
+            {
+                std::cout << std::endl;
+                log << std::endl;
+            }
+        }
+    }
+    
+    //////////////// OUTPUT QUANTITIES TABLE //////////////
+    if (imp == "analytical_Example")   // print Errors only for analytical_Example
+    {
+    std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    std::cout << center("Levels ",tableWidth)       << " | "
+              << center("|q1_ana-q1|",tableWidth)       << " | "
+              << center("|q2_ana-q2|",tableWidth)       << " | "
+              << center("q3",tableWidth)           << " | "
+              << center("|b1_ana-b1|",tableWidth)       << " | "
+              << center("|b2_ana-b2|",tableWidth)       << " | "
+              << center("|b3_ana-b3|",tableWidth)       << " | " << "\n";
+    std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    log       << center("Levels ",tableWidth)       << " | "
+              << center("|q1_ana-q1|",tableWidth)       << " | "
+              << center("|q2_ana-q2|",tableWidth)       << " | "
+              << center("q3",tableWidth)           << " | "
+              << center("|b1_ana-b1|",tableWidth)       << " | "
+              << center("|b2_ana-b2|",tableWidth)       << " | "
+              << center("|b3_ana-b3|",tableWidth)       << " | " << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    }
+    else
+    {
+    std::cout << center("Levels ",tableWidth)       << " | "
+              << center("q1",tableWidth)       << " | "
+              << center("q2",tableWidth)       << " | "
+              << center("q3",tableWidth)           << " | "
+              << center("b1",tableWidth)       << " | "
+              << center("b2",tableWidth)       << " | "
+              << center("b3",tableWidth)       << " | " << "\n";
+    std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";   
+    log       << center("Levels ",tableWidth)       << " | "
+              << center("q1",tableWidth)       << " | "
+              << center("q2",tableWidth)       << " | "
+              << center("q3",tableWidth)           << " | "
+              << center("b1",tableWidth)       << " | "
+              << center("b2",tableWidth)       << " | "
+              << center("b3",tableWidth)       << " | " << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";   
+        
+    }
+    
+    int StorageCount2 = 0;
+    for(auto& v: Storage_Quantities) 
+    {
+        std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth)      << " | ";}, v);
+        std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth)      << " & ";}, v);
+        StorageCount2++;
+        if(StorageCount2 % 7 == 0 )
+        {
+            std::cout << std::endl;
+            log << std::endl;
+        }
+    }
+    std::cout << std::string(tableWidth*7 + 3*7, '-') << "\n";
+    log       << std::string(tableWidth*7 + 3*7, '-') << "\n";  
+
+    log.close();
+}

src/Cell-Problem_muGamma.cc

+#include <config.h>
+#include <array>
+#include <vector>
+#include <fstream>
+
+#include <iostream>
+#include <dune/common/indices.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/math.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/spqr.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/periodicbasis.hh>
+
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/microstructure/prestrain_material_geometry.hh>
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/vtk_filler.hh>    //TEST
+
+
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+
+// #include <dune/fufem/discretizationerror.hh>
+// #include <boost/multiprecision/cpp_dec_float.hpp>
+// #include <boost/any.hpp>
+
+#include <any>
+#include <variant>
+#include <string>
+#include <iomanip>
+
+using namespace Dune;
+using namespace MatrixOperations;
+
+// ------------------------------------------------------------------------
+//
+// This is a stripped-down Version of Cell-Problem.cc 
+// that only computes mu_gamma = q_3 in a faster manner
+// i.e. only the third corrector phi_3 is needed.
+// ------------------------------------------------------------------------ 
+
+
+
+
+
+
+
+
+
+
+
+
+
+//////////////////////////////////////////////////////////////////////
+// Helper functions for Table-Output
+//////////////////////////////////////////////////////////////////////
+/*! Center-aligns string within a field of width w. Pads with blank spaces
+    to enforce alignment. */
+std::string center(const std::string s, const int w) {
+    std::stringstream ss, spaces;
+    int padding = w - s.size();                 // count excess room to pad
+    for(int i=0; i<padding/2; ++i)
+        spaces << " ";
+    ss << spaces.str() << s << spaces.str();    // format with padding
+    if(padding>0 && padding%2!=0)               // if odd #, add 1 space
+        ss << " ";
+    return ss.str();
+}
+
+/* Convert double to string with specified number of places after the decimal
+   and left padding. */
+template<class type>
+std::string prd(const type x, const int decDigits, const int width) {
+    std::stringstream ss;
+//     ss << std::fixed << std::right;
+    ss << std::scientific << std::right;                     // Use scientific Output!
+    ss.fill(' ');        // fill space around displayed #
+    ss.width(width);     // set  width around displayed #
+    ss.precision(decDigits); // set # places after decimal
+    ss << x;
+    return ss.str();
+}
+
+
+
+
+
+template<class Basis>
+auto arbitraryComponentwiseIndices(const Basis& basis,
+                                   const int elementNumber,
+                                   const int leafIdx
+                                   )
+{
+  // (Local Approach -- works for non Lagrangian-Basis too)
+  // Input  : elementNumber & localIdx
+  // Output : determine all Component-Indices that correspond to that basis-function
+  auto localView = basis.localView();
+
+  FieldVector<int,3> row;
+  int elementCounter = 0;
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    if(elementCounter == elementNumber)             // get arbitraryIndex(global) for each Component        ..alternativ:gridView.indexSet
+    {
+      localView.bind(element);
+
+      for (int k = 0; k < 3; k++)
+      {
+        auto rowLocal = localView.tree().child(k).localIndex(leafIdx);    //Input: LeafIndex! TODO bräuchte hier (Inverse )  Local-to-Leaf Idx Map
+        row[k] = localView.index(rowLocal);
+        //         std::cout << "rowLocal:" << rowLocal << std::endl;
+        //         std::cout << "row[k]:" << row[k] << std::endl;
+      }
+    }
+    elementCounter++;
+  }
+  return row;
+}
+
+
+
+
+template<class LocalView, class Matrix, class localFunction1, class localFunction2>
+void computeElementStiffnessMatrix(const LocalView& localView,
+                                   Matrix& elementMatrix,
+                                   const localFunction1& mu,
+                                   const localFunction2& lambda,
+                                   const double gamma
+                                   )
+{
+  // Local StiffnessMatrix of the form:
+  // | phi*phi    m*phi |
+  // | phi *m     m*m   |
+  using Element = typename LocalView::Element;
+  const Element element = localView.element();
+  auto geometry = element.geometry();
+  constexpr int dim = Element::dimension;
+  constexpr int dimWorld = dim;
+  using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+  elementMatrix.setSize(localView.size()+3, localView.size()+3);         //extend by dim ´R_sym^{2x2}
+  elementMatrix = 0;
+
+  // LocalBasis-Offset
+  const int localPhiOffset = localView.size();
+
+  const auto& localFiniteElement = localView.tree().child(0).finiteElement();     
+  const auto nSf = localFiniteElement.localBasis().size();
+//   std::cout << "localView.size(): " << localView.size() << std::endl;
+//   std::cout << "nSf : " << nSf << std::endl;
+
+  ///////////////////////////////////////////////
+  // Basis for R_sym^{2x2}            // wird nicht als Funktion benötigt da konstant...
+  //////////////////////////////////////////////
+  MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+  MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+  MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+  std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+
+//   int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;  // TEST
+  int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+  const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+  for (const auto& quadPoint : quad)
+  {
+    const auto& quadPos = quadPoint.position();
+    const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+    const auto integrationElement = geometry.integrationElement(quadPos);
+
+    std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
+    localFiniteElement.localBasis().evaluateJacobian(quadPos,
+                                                     referenceGradients);
+    // Compute the shape function gradients on the grid element
+    std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+
+    for (size_t i=0; i<gradients.size(); i++)
+      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+    for (size_t k=0; k < dimWorld; k++)
+      for (size_t i=0; i < nSf; i++)
+      {
+        // "phi*phi"-part
+        for (size_t l=0; l< dimWorld; l++)
+          for (size_t j=0; j < nSf; j++ )
+          {
+            // (scaled) Deformation gradient of the ansatz basis function
+            MatrixRT defGradientU(0);
+            defGradientU[k][0] = gradients[i][0];                       // Y
+            defGradientU[k][1] = gradients[i][1];                       //X2
+            defGradientU[k][2] = (1.0/gamma)*gradients[i][2];           //X3
+            // printmatrix(std::cout, defGradientU , "defGradientU", "--");
+            // (scaled)  Deformation gradient of the test basis function
+            MatrixRT defGradientV(0);
+            defGradientV[l][0] = gradients[j][0];                       // Y
+            defGradientV[l][1] = gradients[j][1];                       //X2
+            defGradientV[l][2] = (1.0/gamma)*gradients[j][2];           //X3
+              
+            double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+//             double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);  // also works..
+
+            size_t col = localView.tree().child(k).localIndex(i);                       
+            size_t row = localView.tree().child(l).localIndex(j);
+            elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
+            
+            // "m*phi"  & "phi*m" - part
+            for( size_t m=0; m<3; m++)
+            {
+                double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+
+                auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+                elementMatrix[row][localPhiOffset+m] += value;
+                elementMatrix[localPhiOffset+m][row] += value;
+            }
+          }
+      }
+    // "m*m"-part
+    for(size_t m=0; m<3; m++)
+      for(size_t n=0; n<3; n++)
+      {
+        double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
+        elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
+      }
+  }
+}
+
+
+
+// Get the occupation pattern of the stiffness matrix
+template<class Basis, class ParameterSet>
+void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet& parameterSet)
+{
+  //  OccupationPattern:
+  // | phi*phi    m*phi |
+  // | phi *m     m*m   |
+
+  auto localView = basis.localView();
+  const int phiOffset = basis.dimension();
+
+  nb.resize(basis.size()+3, basis.size()+3);
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    for (size_t i=0; i<localView.size(); i++)
+    {
+      // The global index of the i-th vertex of the element
+      auto row = localView.index(i);
+      for (size_t j=0; j<localView.size(); j++ )
+      {
+        // The global index of the j-th vertex of the element
+        auto col = localView.index(j);
+        nb.add(row[0],col[0]);                   // nun würde auch nb.add(row,col) gehen..
+      }
+    }
+    for (size_t i=0; i<localView.size(); i++)
+    {
+      auto row = localView.index(i);
+
+      for (size_t j=0; j<3; j++ )
+      {
+        nb.add(row,phiOffset+j);               // m*phi part of StiffnessMatrix
+        nb.add(phiOffset+j,row);               // phi*m part of StiffnessMatrix
+      }
+    }
+    for (size_t i=0; i<3; i++ )
+      for (size_t j=0; j<3; j++ )
+      {
+        nb.add(phiOffset+i,phiOffset+j);       // m*m part of StiffnessMatrix
+      }
+  }
+  //////////////////////////////////////////////////////////////////
+  // setOneBaseFunctionToZero
+  //////////////////////////////////////////////////////////////////
+  if(parameterSet.template get<bool>("set_oneBasisFunction_Zero ", true)){
+    FieldVector<int,3> row;
+    unsigned int arbitraryLeafIndex =  parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
+    unsigned int arbitraryElementNumber =  parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
+
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();    // macht keinen Unterschied ob hier k oder 0..
+    const auto nSf = localFiniteElement.localBasis().size();
+    assert(arbitraryLeafIndex < nSf );
+    assert(arbitraryElementNumber  < basis.gridView().size(0));   // "arbitraryElementNumber larger than total Number of Elements"
+
+    //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+    row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
+
+    for (int k = 0; k<3; k++)
+      nb.add(row[k],row[k]);
+  }
+}
+
+
+
+// Compute the source term for a single element
+// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
+template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void computeElementLoadVector( const LocalView& localView,
+                               LocalFunction1& mu,
+                               LocalFunction2& lambda,
+                               const double gamma,
+                               Vector& elementRhs,
+                               const Force& forceTerm
+                               )
+{
+  // | f*phi|
+  // | ---  |
+  // | f*m  |
+  using Element = typename LocalView::Element;
+  const auto element = localView.element();
+  const auto geometry = element.geometry();
+  constexpr int dim = Element::dimension;
+  constexpr int dimWorld = dim;
+
+  using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+  // Set of shape functions for a single element
+  const auto& localFiniteElement= localView.tree().child(0).finiteElement();
+  const auto nSf = localFiniteElement.localBasis().size();
+
+  elementRhs.resize(localView.size() +3);
+  elementRhs = 0;
+
+  // LocalBasis-Offset
+  const int localPhiOffset = localView.size();
+
+  ///////////////////////////////////////////////
+  // Basis for R_sym^{2x2}
+  //////////////////////////////////////////////
+  MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+  MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+  MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+  std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+
+//   int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;  // TEST
+  int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);   // für simplex
+  const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+  for (const auto& quadPoint : quad)
+  {
+    const FieldVector<double,dim>& quadPos = quadPoint.position();
+    const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+    const double integrationElement = geometry.integrationElement(quadPos);
+
+    std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+
+    localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+    std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());  
+    for (size_t i=0; i< gradients.size(); i++)
+      jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+    // "f*phi"-part
+    for (size_t i=0; i < nSf; i++)
+      for (size_t k=0; k < dimWorld; k++)
+      {
+        // Deformation gradient of the ansatz basis function
+        MatrixRT defGradientV(0);
+        defGradientV[k][0] = gradients[i][0];                                 // Y
+        defGradientV[k][1] = gradients[i][1];                                 // X2
+        defGradientV[k][2] = (1.0/gamma)*gradients[i][2];                     // X3
+
+        double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientV, forceTerm(quadPos));
+        size_t row = localView.tree().child(k).localIndex(i);
+        elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
+      }
+
+    // "f*m"-part
+    for (size_t m=0; m<3; m++)
+    {
+      double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], forceTerm(quadPos));
+      elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
+    }
+  }
+}
+
+
+
+template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2, class ParameterSet>
+void assembleCellStiffness(const Basis& basis,
+                           LocalFunction1& muLocal,
+                           LocalFunction2& lambdaLocal,
+                           const double gamma,
+                           Matrix& matrix,
+                           ParameterSet& parameterSet
+                           )
+{
+  std::cout << "assemble Stiffness-Matrix begins." << std::endl;
+
+  MatrixIndexSet occupationPattern;
+  getOccupationPattern(basis, occupationPattern, parameterSet);
+  occupationPattern.exportIdx(matrix);
+  matrix = 0;
+
+  auto localView = basis.localView();
+  const int phiOffset = basis.dimension();
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    muLocal.bind(element);
+    lambdaLocal.bind(element);
+    const int localPhiOffset = localView.size();
+    
+    //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+    Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+    computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
+    //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+    //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
+    //std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
+
+    //////////////////////////////////////////////////////////////////////////////
+    // GLOBAL STIFFNES ASSEMBLY
+    //////////////////////////////////////////////////////////////////////////////
+    for (size_t i=0; i<localPhiOffset; i++)
+    for (size_t j=0; j<localPhiOffset; j++ )
+    {
+        auto row = localView.index(i);
+        auto col = localView.index(j);
+        matrix[row][col] += elementMatrix[i][j];
+    }
+    for (size_t i=0; i<localPhiOffset; i++)
+    for(size_t m=0; m<3; m++)
+    {
+        auto row = localView.index(i);
+        matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m];
+        matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
+    }
+    for (size_t m=0; m<3; m++ )
+    for (size_t n=0; n<3; n++ )
+        matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
+
+    //  printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
+  }
+}
+
+
+template<class Basis, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void assembleCellLoad(const Basis& basis,
+                      LocalFunction1& muLocal,
+                      LocalFunction2& lambdaLocal,
+                      const double gamma,
+                      Vector& b,
+                      const Force& forceTerm
+                      )
+{
+  //     std::cout << "assemble load vector." << std::endl;
+  b.resize(basis.size()+3);
+  b = 0;
+
+  auto localView = basis.localView();
+  const int phiOffset = basis.dimension();
+
+  // Transform G_alpha's to GridViewFunctions/LocalFunctions 
+  auto loadGVF  = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
+  auto loadFunctional = localFunction(loadGVF);      
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    muLocal.bind(element);
+    lambdaLocal.bind(element);
+    loadFunctional.bind(element);
+
+    const int localPhiOffset = localView.size();
+    //         std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+
+    BlockVector<FieldVector<double,1> > elementRhs;
+    computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
+//     printvector(std::cout, elementRhs, "elementRhs", "--");
+//     printvector(std::cout, elementRhs, "elementRhs", "--");
+    //////////////////////////////////////////////////////////////////////////////
+    // GLOBAL LOAD ASSEMBLY
+    //////////////////////////////////////////////////////////////////////////////
+    for (size_t p=0; p<localPhiOffset; p++)
+    {
+      auto row = localView.index(p);
+      b[row] += elementRhs[p];
+    }
+    for (size_t m=0; m<3; m++ )
+      b[phiOffset+m] += elementRhs[localPhiOffset+m];
+      //printvector(std::cout, b, "b", "--");
+  }
+}
+
+
+
+template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
+auto energy(const Basis& basis,
+            LocalFunction1& mu,
+            LocalFunction2& lambda,
+            const MatrixFunction& matrixFieldFuncA,
+            const MatrixFunction& matrixFieldFuncB)
+{
+
+//   TEST HIGHER PRECISION
+//   using float_50 = boost::multiprecision::cpp_dec_float_50;
+//   float_50 higherPrecEnergy = 0.0;
+  double energy = 0.0;
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  constexpr int dimWorld = dim;
+
+  auto localView = basis.localView();
+
+  auto matrixFieldAGVF  = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
+  auto matrixFieldA = localFunction(matrixFieldAGVF);
+  auto matrixFieldBGVF  = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
+  auto matrixFieldB = localFunction(matrixFieldBGVF);
+
+  using GridView = typename Basis::GridView;
+  using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+  using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+//   TEST
+//   FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
+//   printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
+
+  for (const auto& e : elements(basis.gridView()))
+  {
+    localView.bind(e);
+    matrixFieldA.bind(e);
+    matrixFieldB.bind(e);
+    mu.bind(e);
+    lambda.bind(e);
+
+    double elementEnergy = 0.0;
+    //double elementEnergy_HP = 0.0;
+
+    auto geometry = e.geometry();
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+    //int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 );  // TEST
+    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+    const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+
+    for (const auto& quadPoint : quad) 
+    {
+      const auto& quadPos = quadPoint.position();
+      const double integrationElement = geometry.integrationElement(quadPos);
+
+      auto strain1 = matrixFieldA(quadPos);
+      auto strain2 = matrixFieldB(quadPos);
+
+      double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
+
+      elementEnergy += energyDensity * quadPoint.weight() * integrationElement;          
+      //elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
+    }
+    energy += elementEnergy;
+    //higherPrecEnergy += elementEnergy_HP;
+  }
+//   TEST
+//   std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
+  return energy;
+}
+
+
+
+template<class Basis, class Matrix, class Vector, class ParameterSet>   // changed to take only one load... 
+void setOneBaseFunctionToZero(const Basis& basis,
+                              Matrix& stiffnessMatrix,
+                              Vector& load_alpha,
+                              ParameterSet& parameterSet
+                              )
+{
+  std::cout << "Setting one Basis-function to zero" << std::endl;
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+  
+  unsigned int arbitraryLeafIndex =  parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
+  unsigned int arbitraryElementNumber =  parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
+  //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+  FieldVector<int,3> row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
+
+  for (int k = 0; k < dim; k++)
+  {
+    load_alpha[row[k]] = 0.0;
+    auto cIt    = stiffnessMatrix[row[k]].begin();
+    auto cEndIt = stiffnessMatrix[row[k]].end();
+    for (; cIt!=cEndIt; ++cIt)
+      *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
+  }
+}
+
+
+
+template<class Basis>
+auto childToIndexMap(const Basis& basis,
+                     const int k
+                     )
+{
+  // Input  : child/component
+  // Output : determine all Indices that belong to that component
+  auto localView = basis.localView();
+
+  std::vector<int> r = { };
+  //     for (int n: r)
+  //         std::cout << n << ","<< std::endl;
+
+  // Determine global indizes for each component k = 1,2,3.. in order to subtract correct (component of) integral Mean
+  // (global) Indices that correspond to component k = 1,2,3
+  for(const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    const auto& localFiniteElement = localView.tree().child(k).finiteElement();        
+    const auto nSf = localFiniteElement.localBasis().size();                           
+
+    for(size_t j=0; j<nSf; j++)
+    {
+      auto Localidx = localView.tree().child(k).localIndex(j);  // local  indices
+      r.push_back(localView.index(Localidx));                   // global indices
+    }
+  }
+  // Delete duplicate elements
+  // first remove consecutive (adjacent) duplicates
+  auto last = std::unique(r.begin(), r.end());
+  r.erase(last, r.end());
+  // sort followed by unique, to remove all duplicates
+  std::sort(r.begin(), r.end());
+  last = std::unique(r.begin(), r.end());
+  r.erase(last, r.end());
+
+  return r;
+}
+
+
+template<class Basis, class Vector>
+auto integralMean(const Basis& basis,
+                  Vector& coeffVector
+                  )
+{
+  // computes integral mean of given LocalFunction
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+
+  auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
+  auto localfun = localFunction(GVFunction);
+
+  auto localView = basis.localView();
+
+  FieldVector<double,3> r = {0.0, 0.0, 0.0};
+  double area = 0.0;
+
+  // Compute Area integral & integral of FE-function
+  for(const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    localfun.bind(element);
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+    const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+
+    for(const auto& quadPoint : quad)
+    {
+      const auto& quadPos = quadPoint.position();
+      const double integrationElement = element.geometry().integrationElement(quadPos);
+      area += quadPoint.weight() * integrationElement;
+      
+      r += localfun(quadPos) * quadPoint.weight() * integrationElement;
+    }
+  }
+  //     std::cout << "Domain-Area: " << area << std::endl;
+  return (1.0/area) * r ;
+}
+
+
+template<class Basis, class Vector>
+auto subtractIntegralMean(const Basis& basis,
+                          Vector& coeffVector
+                          )
+{
+  // Substract correct Integral mean from each associated component function
+  auto IM = integralMean(basis, coeffVector);
+
+  constexpr int dim = Basis::LocalView::Element::dimension;
+
+  for(size_t k=0; k<dim; k++)
+  {
+    //std::cout << "Integral-Mean: " << IM[k] << std::endl;
+    auto idx = childToIndexMap(basis,k);
+
+    for ( int i : idx)
+      coeffVector[i] -= IM[k];
+  }
+}
+
+
+
+//////////////////////////////////////////////////
+//   Infrastructure for handling periodicity
+//////////////////////////////////////////////////
+
+// Check whether two points are equal on R/Z x R/Z x R
+auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
+                {
+                    return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+                            and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
+                            and (FloatCmp::eq(x[2],y[2]))
+                        );
+                };
+
+
+        
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[])
+{
+  MPIHelper::instance(argc, argv);
+
+  ParameterTree parameterSet;
+  if (argc < 2)
+    ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
+  else
+  {
+    ParameterTreeParser::readINITree(argv[1], parameterSet);
+    ParameterTreeParser::readOptions(argc, argv, parameterSet);
+  }
+
+  // Output setter
+//   std::string outputPath = parameterSet.get("outputPath", "../../outputs/output.txt");
+//   std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt");
+  std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs");
+//   std::string MatlabPath = parameterSet.get("MatlabPath", "/home/klaus/Desktop/DUNE/dune-microstructure/Matlab-Programs");
+//     std::string outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt";
+  std::fstream log;
+  log.open(outputPath + "/output.txt" ,std::ios::out);
+
+  std::cout << "outputPath:" << outputPath << std::endl;
+  
+//   parameterSet.report(log); // short Alternativ
+  
+  constexpr int dim = 3;
+  constexpr int dimWorld = 3;
+
+  ///////////////////////////////////
+  // Get Parameters/Data
+  ///////////////////////////////////
+  double gamma = parameterSet.get<double>("gamma",3.0);   // ratio dimension reduction to homogenization
+  double alpha = parameterSet.get<double>("alpha", 2.0);
+  double theta = parameterSet.get<double>("theta",1.0/4.0); 
+  ///////////////////////////////////
+  // Get Material Parameters
+  ///////////////////////////////////
+  std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
+  double beta = parameterSet.get<double>("beta",2.0); 
+  double mu1 = parameterSet.get<double>("mu1",10.0);;
+  double mu2 = beta*mu1;
+  double lambda1 = parameterSet.get<double>("lambda1",0.0);;
+  double lambda2 = beta*lambda1;
+  
+  // Plate geometry settings
+  double width = parameterSet.get<double>("width", 1.0);   //geometry cell, cross section
+  //     double len  = parameterSet.get<double>("len", 10.0); //length
+  //     double height  = parameterSet.get<double>("h", 0.1); //rod thickness
+  //     double eps  = parameterSet.get<double>("eps", 0.1); //size of perioticity cell
+
+  ///////////////////////////////////
+  // Get Prestrain/Parameters
+  ///////////////////////////////////
+//   unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", "analytical_Example");  //OLD 
+//   std::string prestraintype = parameterSet.get<std::string>("prestrainType", "analytical_Example");
+//   double rho1 = parameterSet.get<double>("rho1", 1.0);
+//   double rho2 = alpha*rho1;
+//   auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
+//   auto B_Term = prestrainImp.getPrestrain(prestraintype);
+  
+  
+  
+  log << "----- Input Parameters -----: " << std::endl;
+//   log << "prestrain imp: " <<  prestraintype << "\nrho1 = " << rho1 << "\nrho2 = " << rho2  << std::endl;
+  log << "alpha: " << alpha << std::endl;
+  log << "gamma: " << gamma << std::endl;
+  log << "theta: " << theta << std::endl;
+  log << "beta: " << beta << std::endl;
+  log << "material parameters: " << std::endl;
+  log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+  log << "lambda1: " << lambda1 <<"\nlambda2: " << lambda2 << std::endl;
+  log << "----------------------------: " << std::endl;
+
+  ///////////////////////////////////
+  // Generate the grid
+  ///////////////////////////////////
+
+  //Corrector Problem Domain:
+  unsigned int cellDomain = parameterSet.get<unsigned int>("cellDomain", 1);
+  // (shifted) Domain (-1/2,1/2)^3
+  FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
+  FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
+  if (cellDomain == 2)
+  {
+    // Domain : [0,1)^2 x (-1/2, 1/2)
+    FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
+    FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
+  }
+
+  std::array<int,2> numLevels = parameterSet.get<std::array<int,2>>("numLevels", {1,3});
+
+  int levelCounter = 0;
+
+  ///////////////////////////////////
+  // Create Data Storage
+  ///////////////////////////////////
+  // Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4  L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
+  std::vector<std::variant<std::string, size_t , double>> Storage_Error;
+  
+  // Storage:: #1 level #2 |q1_a-q1_c| #3 |q2_a-q2_c| #4 |q3_a-q3_c| #5 |b1_a-b1_c| #6 |b2_a-b2_c| #7 |b3_a-b3_c|           
+   std::vector<std::variant<std::string, size_t , double>> Storage_Quantities;         
+
+
+//   for(const size_t &level : numLevels)                             // explixite Angabe der levels.. {2,4}
+  for(size_t level = numLevels[0] ; level <= numLevels[1]; level++)    // levels von bis.. [2,4]
+  {
+    std::cout << " ----------------------------------" << std::endl;
+    std::cout << "Level: " << level << std::endl;
+    std::cout << " ----------------------------------" << std::endl;
+
+    Storage_Error.push_back(level);
+    Storage_Quantities.push_back(level);
+    std::array<int, dim> nElements = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
+
+    std::cout << "Number of Elements in each direction: " << nElements << std::endl;
+    log << "Number of Elements in each direction: " << nElements << std::endl;
+
+    using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+    CellGridType grid_CE(lower,upper,nElements);
+    using GridView = CellGridType::LeafGridView;
+    const GridView gridView_CE = grid_CE.leafGridView();
+
+    using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+    using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+    using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+    using VectorCT = BlockVector<FieldVector<double,1> >;
+    using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
+
+    ///////////////////////////////////
+    //  Create Lambda-Functions for material Parameters depending on microstructure
+    ///////////////////////////////////
+    
+    auto materialImp = IsotropicMaterialImp();
+    auto muTerm = materialImp.getMu(parameterSet);
+    auto lambdaTerm = materialImp.getLambda(parameterSet);
+
+    auto muGridF  = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
+    auto muLocal = localFunction(muGridF);
+    auto lambdaGridF  = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
+    auto lambdaLocal = localFunction(lambdaGridF);
+
+    ///////////////////////////////////
+    // --- Choose Solver ---
+    // 1 : CG-Solver
+    // 2 : GMRES
+    // 3 : QR
+    ///////////////////////////////////
+    unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+    
+    // Print Options 
+    bool print_debug = parameterSet.get<bool>("print_debug", false);
+    
+    
+    bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
+    bool write_prestrainFunctions  = parameterSet.get<bool>("write_prestrainFunctions", false);
+    
+    
+    bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
+    bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
+    bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
+
+    /////////////////////////////////////////////////////////
+    // Choose a finite element space for Cell Problem
+    /////////////////////////////////////////////////////////
+    using namespace Functions::BasisFactory;
+    Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+
+    // Don't do the following in real life: It has quadratic run-time in the number of vertices.
+    for (const auto& v1 : vertices(gridView_CE))
+        for (const auto& v2 : vertices(gridView_CE))
+            if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
+            {
+                periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
+            }
+
+    // First order periodic Lagrange-Basis
+    auto Basis_CE = makeBasis(
+        gridView_CE,
+        power<dim>(       
+        Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
+        flatLexicographic()));     
+
+    log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
+    const int phiOffset = Basis_CE.size();
+
+    /////////////////////////////////////////////////////////
+    // Data structures: Stiffness matrix and right hand side vectors
+    /////////////////////////////////////////////////////////
+    VectorCT load_alpha1, load_alpha2, load_alpha3;
+    MatrixCT stiffnessMatrix_CE;
+
+    bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
+    bool set_oneBasisFunction_Zero = false;
+    bool substract_integralMean = false;
+    if(set_IntegralZero)
+    {
+        set_oneBasisFunction_Zero = true;
+        substract_integralMean = true;
+    }
+
+    /////////////////////////////////////////////////////////
+    // Define "rhs"-functions
+    /////////////////////////////////////////////////////////
+    Func2Tensor x3G_3 = [] (const Domain& x) {                                                                               
+                            return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+                        };
+
+    Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {return -1.0*x3G_3(x);};
+    
+    
+
+    ///////////////////////////////////////////////
+    // Basis for R_sym^{2x2}
+    //////////////////////////////////////////////
+    MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+    MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+    MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+    std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+    //  printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+    //  printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+    //  printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
+    //  log <<  basisContainer[0] << std::endl;
+
+
+    //////////////////////////////////////////////////////////////////////
+    // Determine global indices of components for arbitrary (local) index
+    //////////////////////////////////////////////////////////////////////
+    int arbitraryLeafIndex =  parameterSet.get<unsigned int>("arbitraryLeafIndex", 0);            // localIdx..assert Number < #ShapeFcts 
+    int arbitraryElementNumber =  parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
+
+    FieldVector<int,3> row;
+     row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
+    if(print_debug)
+      printvector(std::cout, row, "row" , "--");
+
+    /////////////////////////////////////////////////////////
+    // Assemble the system
+    /////////////////////////////////////////////////////////
+    assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma,  stiffnessMatrix_CE, parameterSet);
+    assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg);      // only third corrector is needed
+    //   printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
+    //   printvector(std::cout, load_alpha1, "load_alpha1", "--");
+
+
+
+    // set one basis-function to zero
+    if(set_oneBasisFunction_Zero)
+    {
+        setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha3, parameterSet);
+    }
+
+    ////////////////////////////////////////////////////
+    // Compute solution
+    ////////////////////////////////////////////////////
+    VectorCT x_3 = load_alpha3;
+
+    if (Solvertype == 1)  // CG - SOLVER
+    {
+        std::cout << "------------ CG - Solver ------------" << std::endl;
+        MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
+
+
+
+        // Sequential incomplete LU decomposition as the preconditioner
+        SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_CE,1.0);
+        int iter = parameterSet.get<double>("iterations_CG", 999);
+        // Preconditioned conjugate-gradient solver
+        CGSolver<VectorCT> solver(op,
+                                ilu0, //NULL,
+                                1e-8, // desired residual reduction factorlack
+                                iter,   // maximum number of iterations
+                                2);   // verbosity of the solver
+        InverseOperatorResult statistics;
+//         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";
+//         solver.apply(x_2, load_alpha2, statistics);
+        std::cout << "solve linear system for x_3.\n";
+        solver.apply(x_3, load_alpha3, statistics);
+        log << "Solver-type used: " <<" CG-Solver" << std::endl;
+    }
+    ////////////////////////////////////////////////////////////////////////////////////
+
+    else if (Solvertype ==2)  // GMRES - SOLVER
+    {
+
+        std::cout << "------------ GMRES - Solver ------------" << std::endl;
+        // Turn the matrix into a linear operator
+        MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
+
+        // Fancy (but only) way to not have a preconditioner at all
+        Richardson<VectorCT,VectorCT> preconditioner(1.0);
+
+        // Construct the iterative solver
+        RestartedGMResSolver<VectorCT> solver(
+            stiffnessOperator, // Operator to invert
+            preconditioner,    // Preconditioner
+            1e-10,             // Desired residual reduction factor
+            500,               // Number of iterations between restarts,
+                            // here: no restarting
+            500,               // Maximum number of iterations
+            2);                // Verbosity of the solver
+
+        // Object storing some statistics about the solving process
+        InverseOperatorResult statistics;
+
+
+        solver.apply(x_3, load_alpha3, statistics);
+        log << "Solver-type used: " <<" GMRES-Solver" << std::endl;
+    }
+    ////////////////////////////////////////////////////////////////////////////////////
+    else if ( Solvertype == 3)// QR - SOLVER
+    {
+
+        std::cout << "------------ QR - Solver ------------" << std::endl;
+        log << "solveLinearSystems: We use QR solver.\n";
+        //TODO install suitesparse
+        SPQR<MatrixCT> sPQR(stiffnessMatrix_CE);
+        sPQR.setVerbosity(1);
+        InverseOperatorResult statisticsQR;
+
+        sPQR.apply(x_3, load_alpha3, statisticsQR);
+        log << "Solver-type used: " <<" QR-Solver" << std::endl;
+    }
+
+    ////////////////////////////////////////////////////////////////////////////////////
+    // Extract phi_alpha  &  M_alpha coefficients
+    ////////////////////////////////////////////////////////////////////////////////////
+    VectorCT phi_3;
+    phi_3.resize(Basis_CE.size());
+    phi_3 = 0;
+
+    for(size_t i=0; i<Basis_CE.size(); i++)
+    {
+
+        phi_3[i] = x_3[i];
+    }
+
+    FieldVector<double,3> m_1, m_2, m_3;
+
+    for(size_t i=0; i<3; i++)
+    {
+        m_3[i] = x_3[phiOffset+i];
+    }
+    // assemble M_alpha's (actually iota(M_alpha) )
+
+    MatrixRT M_1(0), M_2(0), M_3(0);
+
+    for(size_t i=0; i<3; i++)
+    {
+        M_3 += m_3[i]*basisContainer[i];
+    }
+
+    std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
+//     printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
+    log << "---------- OUTPUT ----------" << std::endl;
+    log << " --------------------" << std::endl;
+    log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
+    log << " --------------------" << std::endl;
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Substract Integral-mean
+    ////////////////////////////////////////////////////////////////////////////                                                                                    
+    using SolutionRange = FieldVector<double,dim>;
+
+    if(substract_integralMean)
+    {
+        std::cout << " --- substracting integralMean --- " << std::endl;
+        subtractIntegralMean(Basis_CE,  phi_3);
+        subtractIntegralMean(Basis_CE,  x_3);
+    }
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Make a discrete function from the FE basis and the coefficient vector
+    ////////////////////////////////////////////////////////////////////////////                                                                                     // 
+//     auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+//         Basis_CE,
+//         phi_3);
+
+
+    /////////////////////////////////////////////////////////
+    // Compute effective quantities: Elastic law & Prestrain
+    /////////////////////////////////////////////////////////
+    VectorCT tmp1;
+    assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3G_3);
+    double GGterm = 0.0;
+    GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3G_3, x3G_3 );
+    auto mu_gamma = (x_3*tmp1) + GGterm;
+    log << "mu_gamma=" << mu_gamma << std::endl;  
+    
+    std::cout << "mu_gamma:" << mu_gamma << std::endl;
+    Storage_Quantities.push_back(mu_gamma);
+  
+    
+  levelCounter++; 
+  } // Level-Loop End
+  
+
+    log.close();
+}

src/CellScript.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import time
+from ClassifyMin import *
+from HelperFunctions import *
+# from scipy.io import loadmat #not Needed anymore?
+import codecs
+import sys
+
+
+# ----------------------------------------------------------------------------------------------------------------------------
+
+
+
+
+
+
+# ----- Setup Paths -----
+InputFile  = "/inputs/cellsolver.parset"
+OutputFile = "/outputs/output.txt"
+# --------- Run  from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+
+print('---- Input parameters: -----')
+mu1 = 10.0
+rho1 = 1.0
+alpha = 2.8
+beta = 2.0
+theta = 1.0/4.0
+gamma = 0.75
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+
+
+
+print('RunCellProblem...')
+RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
+
+print('Read effective quantities...')
+Q, B = ReadEffectiveQuantities()
+print('Q:', Q)
+print('B:', B)
+
+
+
+# Compare symbolicMinimization with Classification 'ClassifyMin' :
+# print('Compare_Classification...')
+# Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
+
+
+
+# ------------- RUN Matlab symbolic minimization program 'symMinimization'
+eng = matlab.engine.start_matlab()
+# s = eng.genpath(path + '/Matlab-Programs')
+s = eng.genpath(path)
+eng.addpath(s, nargout=0)
+# print('current Matlab folder:', eng.pwd(nargout=1))
+eng.cd('Matlab-Programs', nargout=0)  #switch to Matlab-Programs folder
+# print('current Matlab folder:', eng.pwd(nargout=1))
+Inp = False
+Inp_T = True
+print('Run symbolic Minimization...')
+#Arguments: symMinization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath)
+G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4)  #Name of program:symMinimization
+# G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4)  #Optional: add Path
+G = np.asarray(G) #cast Matlab Outout to numpy array
+# --- print Output ---
+print('Minimizer G:')
+print(G)
+print('angle:', angle)
+print('type:', type )
+print('curvature:', kappa)

src/ClassifyMin.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+# print(sys.executable)
+# from subprocess import Popen, PIPE
+
+
+# --------------------------------------------------
+# 'classifyMin' classifies Minimizers by utilizing the result of
+# Lemma1.6
+#
+#
+#
+#
+# 'classifyMin_ana': (Special-Case : Lemma1.4)
+# ..additionally assumes Poisson-ratio=0 => q12==0
+#
+#
+#
+# Output : MinimizingMatrix, Angle, Type, Curvature
+
+
+
+
+
+def harmonicMean(mu_1, beta, theta):
+    return mu_1*(beta/(theta+(1-theta)*beta))
+
+
+def arithmeticMean(mu_1, beta, theta):
+    return mu_1*((1-theta)+theta*beta)
+
+
+def prestrain_b1(rho_1, beta, alpha, theta):
+    return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)))
+
+
+def prestrain_b2(rho_1, beta, alpha, theta):
+    return (3.0*rho_1/(4.0*((1.0-theta) + theta*beta)))*(1-theta*(1+beta*alpha))
+
+
+# Define function to be minimized
+def f(a1, a2, q1, q2, q3, q12, b1, b2):
+    A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
+    B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
+    a = np.array([a1, a2])
+    tmp = np.dot(A, a)
+    tmp2 = np.dot(a, tmp)
+    tmpB = np.dot(B, a)
+    return tmp2 + tmpB + q1*(b1**2) + q2*(b2**2) + q12*b1*b2
+
+# ---- Alternative Version using alpha,beta,theta ,mu_1,rho_1
+def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Output=False):
+
+    # Assumption of Classification-Lemma1.6:
+    #  1. [b3 == 0]
+    #  2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
+    # 3. This additionally assumes that Poisson-Ratio = 0 => q12 == 0
+    q12 = 0.0
+
+    q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
+    q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
+    # print('q1: ', q1)
+    # print('q2: ', q2)
+    b1 = prestrain_b1(rho_1, beta, alpha,theta)
+    b2 = prestrain_b2(rho_1, beta, alpha,theta)
+    return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
+
+
+
+# Matrix Version that just gets matrices Q & B
+def classifyMin_mat(Q,B,print_Cases=False, print_Output=False):
+    q1 = Q[0][0]
+    q2 = Q[1][1]
+    q3 = Q[2][2]
+    q12 = Q[0][1]
+    b1 = B[0]
+    b2 = B[1]
+    b3 = B[2]
+    return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
+
+
+# --------------------------------------------------------------------
+# Classify Type of minimizer  1 = R1 , 2 = R2 , 3 = R3                          # before : destinction between which axis.. (4Types )
+# where
+# R1 : unique local (global) minimizer which is not axial
+# R2 : continuum of local (global) minimizers which are not axial
+# R3 : one or two local (global) minimizers which are axial
+# Partition given by
+# R1 = E1
+# R2 = P1.2
+# R3 = E2 U E3 U P1.1 U P2 U H
+# -------------------------------------------------------------------
+def classifyMin(q1, q2, q3, q12, b1, b2,  print_Cases=False, print_Output=False):   #ClassifyMin_hom?
+    # Assumption of Classification-Lemma1.6:
+    #  1. [b3 == 0]
+    #  2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
+
+    # TODO: check if Q is orthotropic here - assert()
+
+
+    if print_Output: print("Run ClassifyMin...")
+    CaseCount = 0
+    epsilon = sys.float_info.epsilon #Machine epsilon
+
+    B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
+    A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
+    # print('Matrix A:',  A)
+    # print('Matrix B:', B)
+    # print('Matrix rank of A:', np.linalg.matrix_rank(A))
+    # print('shape of A:', A.shape)
+    # print('shape of B:', B.shape)
+    # print('Matrix [A,B]:', np.c_[A, B])
+    # print('Matrix rank of [A,B]:', np.linalg.matrix_rank(np.c_[A, B]))
+    # print('shape of [A,B]:', C.shape)
+    #
+    # x = np.linalg.solve(A, B)     # works only if A is not singular!!
+    # print("Solution LGS:", x)
+    # # print("sym Matrix", sym.Matrix(([A],[B])))
+    #
+    # # Test
+    # C = np.array([[1, 0], [0, 0]])
+    # d = np.array([5, 0])
+    # y = np.linalg.lstsq(C, d)[0]
+    # print("Solution LGS:", y)
+    # T = np.c_[C, d]
+    # print('T:', T)
+    # Trref = sym.Matrix(T).rref()[0]
+    # Out = np.array(Trref, dtype=float)
+    # print('rref:', Out)
+
+    determinant = q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
+    if print_Cases: print("determinant:", determinant)
+
+    # Define values for axial-Solutions (b1*,0) & (0,b2*)
+    b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
+    b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
+
+    # ------------------------------------ Parabolic Case -----------------------------------
+    if abs(determinant) < epsilon:
+        if print_Cases: print('P : parabolic case (determinant equal zero)')
+        # if print_Cases: print('P : parabolic case (determinant equal zero)')
+
+        # check if B is in range of A
+        # check if rank(A) == rank([A,B])
+        # OK this way? (or use Sympy?)
+        if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.c_[A, B]):
+            if print_Cases: print('P1 (B is in the range of A)')
+
+            if (q12+q3)/2.0 <= -1.0*epsilon:
+                print('should not happen(not compatible with det = 0)')
+
+            if (abs(B[0]) < epsilon and abs(B[1]) < epsilon) and (q12+q3)/2.0 >= epsilon:
+                if print_Cases: print('P1.1')
+                a1 = 0.0
+                a2 = 0.0
+                type = 3
+                CaseCount += 1
+
+            if (abs(B[0]) >= epsilon or abs(B[1]) >= epsilon) and (q12+q3)/2.0 >= epsilon:
+                # Continuum of minimizers located along a line of negative slope in Lambda
+                if print_Cases: print('P1.2 (Continuum of minimizers located along a line of negative slope in Lambda) ')
+
+                # Just solve Aa* = b (alternatively using SymPy ?)
+                # we know that A is singular (det A = 0 ) here..
+                # & we know that there are either infinitely many solutions or a unique solution ...
+                # ---- determine one via Least Squares
+                # "If A is square and of full rank, then x (but for round-off error)
+                # is the “exact” solution of the equation. Else, x minimizes the
+                # Euclidean 2-norm || b-Ax ||. If there are multiple minimizing solutions,
+                # the one with the smallest 2-norm  is returned. ""
+                a = np.linalg.lstsq(A, B)[0]                            # TODO check is this Ok ?
+                print("Solution LGS: a =", a)
+                a1 = a[0]
+                a2 = a[1]
+
+                type = 2
+                CaseCount += 1
+        else:
+            if print_Cases: print('P2 (B is not in the range of A)')
+            # local Minimizers occur on the boundary of Lambda...
+            # There are at most two local minima and they are either
+            # (b1_star, 0) or (0, b2_star)
+            # could also outsource this to another function..
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 1
+                CaseCount += 1
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = 0
+                a2 = b2_star
+                type = 3  # 2
+                CaseCount += 1
+            # TODO Problem: angle depends on how you choose... THE angle is not defined for this case
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                # Two Minimizers pick one
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 4
+                CaseCount += 1
+
+    # ------------------------------------ Elliptic Case -----------------------------------
+    if determinant >= epsilon:
+        if print_Cases: print('E : elliptic case (determinant greater zero)')
+        a1_star = -(2*(b1*(q12**2) + 2*b1*q3*q12 - 4*b1*q1*q2 + 4*b2*q2*q3)) / \
+            (4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+        a2_star = -(2*(b2*(q12**2) + 2*b2*q3*q12 + 4*b1*q1*q3 - 4*b2*q1*q2)) / \
+            (4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+        prod = a1_star*a2_star
+
+        if prod >= epsilon:
+            if print_Cases: print('(E1) - inside Lambda ')
+            a1 = a1_star
+            a2 = a2_star
+            type = 1  # non-axial Minimizer
+            CaseCount += 1
+        if abs(prod) < epsilon:   # same as prod = 0 ? or better use <=epsilon ?
+            if print_Cases: print('(E2) - on the boundary of Lambda ')
+            a1 = a1_star
+            a2 = a2_star
+            type = 3  # could check which axis: if a1_star or a2_star close to zero.. ?
+            CaseCount += 1
+
+        if prod <= -1.0*epsilon:
+            if print_Cases: print('(E3) - Outside Lambda ')
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 1
+                CaseCount += 1
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                a1 = 0
+                a2 = b2_star
+                type = 3  # 2
+                CaseCount += 1
+
+            # TODO Problem: angle depends on how you choose... THE angle is not defined for this case
+            if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
+                # Two Minimizers pick one
+                a1 = b1_star
+                a2 = 0.0
+                type = 3  # 4
+                CaseCount += 1
+
+    # ------------------------------------ Hyperbolic Case -----------------------------------
+    if determinant <= -1.0*epsilon:
+        if print_Cases: print('H : hyperbolic case (determinant smaller zero)')
+        # One or two minimizers wich are axial
+        type = 3 # (always type 3)
+        if f(b1_star, 0, q1, q2, q3, q12, b1, b2) < f(0, b2_star, q1, q2, q3, q12, b1, b2):
+            a1 = b1_star
+            a2 = 0.0
+            # type = 3  # 1
+            CaseCount += 1
+        if f(b1_star, 0, q1, q2, q3, q12, b1, b2) > f(0, b2_star, q1, q2, q3, q12, b1, b2):
+            a1 = 0
+            a2 = b2_star
+            # type = 3  # 2
+            CaseCount += 1
+        # TODO can add this case to first or second ..
+        if f(b1_star, 0, q1, q2, q3, q12, b1, b2) == f(0, b2_star, q1, q2, q3, q12, b1, b2):
+            # Two Minimizers pick one
+            a1 = b1_star
+            a2 = 0.0
+            # type = 3  # 4
+            CaseCount += 1
+    # ---------------------------------------------------------------------------------------
+
+    if (CaseCount > 1):
+        print('Error: More than one Case happened!')
+
+    # compute a3
+    # a3 = math.sqrt(2.0*a1*a2)   # never needed?
+
+    # compute the angle <(e,e_1) where Minimizer = kappa* (e (x) e)
+    e = [math.sqrt((a1/(a1+a2))), math.sqrt((a2/(a1+a2)))]
+    angle = math.atan2(e[1], e[0])
+
+    # compute kappa
+    kappa = (a1 + a2)
+
+    # Minimizer G
+    Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]],dtype=object)
+    # Minimizer = np.array([[a1, math.sqrt(a1*a2)], [math.sqrt(a1*a2), a2]])
+
+    if print_Output:
+        print('--- Output ClassifyMin ---')
+        print("Minimizing Matrix G:")
+        print(Minimizer)
+        print("angle = ", angle)
+        print("type: ", type)
+        print("kappa = ", kappa)
+
+    return Minimizer, angle, type, kappa
+# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+
+# ---------------------------------------------- Main ---------------------
+
+# --- Input Parameters ----
+# mu_1 = 1.0
+# rho_1 = 1.0
+# alpha = 9.0
+# beta = 2.0
+# theta = 1.0/8.0
+
+# # define q1, q2 , mu_gamma, q12
+# # 1. read from Cell-Output
+# # 2. define values from analytic formulas (expect for mu_gamma)
+# q1 = (1.0/6.0)*harmonicMean(mu_1, beta, theta)
+# q2 = (1.0/6.0)*arithmeticMean(mu_1, beta, theta)
+# # TEST
+# q12 = 0.0  # (analytical example)
+# # q12 = 12.0  # (analytical example)
+
+# set mu_gamma  to value or read from Cell-Output
+# mu_gamma = q1   # TODO read from Cell-Output
+# b1 = prestrain_b1(rho_1, beta, alpha, theta)
+# b2 = prestrain_b2(rho_1, beta, alpha, theta)
+
+
+# print('---- Input parameters: -----')
+# print('mu_1: ', mu_1)
+# print('rho_1: ', rho_1)
+# print('alpha: ', alpha)
+# print('beta: ', beta)
+# print('theta: ', theta)
+# print("q1: ", q1)
+# print("q2: ", q2)
+# print("mu_gamma: ", mu_gamma)
+# print("q12: ", q12)
+# print("b1: ", b1)
+# print("b2: ", b2)
+# print('----------------------------')
+# # print("machine epsilon", sys.float_info.epsilon)
+#
+#
+# # ------- Options --------
+# print_Cases = True
+# print_Output = True
+
+# G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12,  b1, b2, print_Cases, print_Output)
+#
+# G, angle, type, kappa = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+#
+# Out = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+
+# print('TEST:')
+# Out = classifyMin_ana(alpha, beta, theta)
+
+# print('Out[0]', Out[0])
+# print('Out[1]', Out[1])
+# print('Out[2]', Out[2])
+# print('Out[3]', Out[3])
+
+
+# #supress certain Outout..
+# _,_,T,_ = classifyMin_ana(alpha, beta, theta, mu_gamma, q12, print_Cases, print_Output)
+# print('Output only type..:', T)
+
+# Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+# print("Test", Test)
+
+# -----------------------------------------------------------------------------------------------------------------------------------------------------------------

src/Compute_MuGamma.cc

+#include <config.h>
+
+#include <vector>
+
+#include <dune/geometry/quadraturerules.hh>
+
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/istl/io.hh>
+#include <dune/istl/matrix.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/matrixmarket.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>  // use for LagrangeBasis
+
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/math.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+
+#include <dune/functions/functionspacebases/periodicbasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+
+#include <dune/microstructure/prestrain_material_geometry.hh>
+// #include <dune/microstructure/matrix_operations.hh>
+
+
+// #include <math.h>
+#include <cmath>
+// #include <numbers>
+
+using namespace Dune;
+
+
+
+// ------------------------------------------------------------------------
+//
+// Solving the 2D "Poisson-Type" Equation ( (38) in the draft)
+// in order to compute mu_Gamma in a faster manner 
+//
+// ------------------------------------------------------------------------ 
+
+
+
+template<class LocalView, class Matrix, class LocalFunction>
+void assembleElementStiffnessMatrix(const LocalView& localView,
+                                    Matrix& elementMatrix,
+                                    const LocalFunction& mu,
+                                    const double gamma
+                                   )
+{
+    using Element = typename LocalView::Element;
+    constexpr int dim = Element::dimension;
+    const Element element = localView.element();
+    auto geometry = element.geometry();
+
+
+    const auto& localFiniteElement = localView.tree().finiteElement();
+    const auto nSf = localFiniteElement.localBasis().size();
+    
+    
+//     std::cout << "nSf:" << nSf << std::endl;
+    
+    elementMatrix.setSize(localView.size(),localView.size());
+    elementMatrix = 0;
+
+//     int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1) + 5;   // doesnt change anything
+    int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+//     int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+    const auto& quadRule = QuadratureRules<double, dim>::rule(element.type(),orderQR);
+    
+//     std::cout << "QuadRule-Order:" << orderQR << std::endl;
+    
+//     double muValue = 0.0;
+//     std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+
+    for (const auto& quadPoint : quadRule)
+    {
+        
+    
+        const auto& quadPos = quadPoint.position();
+        
+//         std::cout << " quadPOS: " << quadPos << std::endl;
+                // TEST 
+//         double muValue = mu(quadPos);
+//         std::cout << "muValue:" << muValue << std::endl;
+
+        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+        const double integrationElement = geometry.integrationElement(quadPos);
+
+
+        
+        std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+        localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+        
+        // Compute the shape function gradients on the grid element
+        std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+        
+        for (size_t i=0; i<gradients.size(); i++)
+          jacobian.mv(referenceGradients[i][0], gradients[i]);                                 //TODO ? passt..
+   
+        
+//         // print all gradients  //TEST 
+//         FieldVector<double,dim> tmp = {0.0 , 0.0};
+//         for (size_t i=0; i < nSf; i++)
+//         {
+//             printvector(std::cout, gradients[i], "gradients[i]", "--" ); 
+// 
+//         
+//             tmp[0] += gradients[i][0];
+//             tmp[1] += gradients[i][1];
+//         
+// //         tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
+// //         tmp[1] += coeffVector[globalIdx]*gradients[i][2];
+// //         printvector(std::cout, tmp, "tmp", "--" ); 
+// 
+//         }
+//         printvector(std::cout, tmp, "sum of basisfunction Gradients", "--" ); 
+        
+        for (size_t p=0; p<elementMatrix.N(); p++)
+        {
+            for (size_t q=0; q<elementMatrix.M(); q++)
+            {
+//                 auto localRow = localView.tree().localIndex(p);  // VERTAUSCHT?!?!
+//                 auto localCol = localView.tree().localIndex(q);
+                auto localRow = localView.tree().localIndex(q);
+                auto localCol = localView.tree().localIndex(p);
+                
+//                 auto value = mu(quadPos)*(2.0* gradients[p][0] * gradients[q][0])+ mu(quadPos)*((1.0/(std::pow(gamma,2)))*(gradients[p][1] * gradients[q][1]));
+//                 auto value = (mu(quadPos)*gradients[p][0] * gradients[q][0])+ ((1.0/gamma)*(gradients[p][1] * gradients[q][1]));
+                
+
+                
+//                 std::cout << "gradients[p][0]" << gradients[p][0] << std::endl; 
+//                 std::cout << "gradients[q][0]" << gradients[q][0] << std::endl; 
+//                 std::cout << "(1.0/std::pow(gamma,2))*gradients[p][1]" << (1.0/std::pow(gamma,2))*gradients[p][1]<< std::endl; 
+                
+                
+
+//                 auto value3 = mu(quadPos)*(1.0/gamma)*gradients[p][2];   // 3D-Version
+//                 auto value4 = value3*gradients[q][2];                    // 3D-Version
+                
+                auto value1 = 2.0*mu(quadPos)* gradients[p][0] *gradients[q][0];         //#1 
+//                 auto value1 = 2.0*mu(quadPos)*sqrt(2.0)* gradients[p][0] *gradients[q][0];                          // TEST TODO warum passt es damit besser bei gamma = 1.0 ??? 
+//                 auto value2 = mu(quadPos)*(1.0/std::pow(gamma,2))*gradients[p][1] * gradients[q][1] ;  
+                
+                
+            
+                
+                
+                auto value2 = 2.0*mu(quadPos)*(1.0/std::pow(gamma,2))*gradients[p][1] * gradients[q][1] ;  //#1  TEST FAKTOR 2 hat hier gefehlt?!?! 
+                
+                auto value = value1 + value2;
+                
+                elementMatrix[localRow][localCol] += value * quadPoint.weight() * integrationElement;
+            }
+        }
+    }
+}
+
+
+
+
+// Compute the source term for a single element
+template<class LocalView,class Vector, class LocalFunction, class Force>
+void assembleElementVolumeTerm(const LocalView& localView,
+                               Vector& elementRhs,
+                               LocalFunction& mu,
+                               const Force& forceTerm)
+{
+    using Element = typename LocalView::Element;
+    auto element = localView.element();
+    auto geometry = element.geometry();
+    
+
+    constexpr int dim = Element::dimension;
+    // Set of shape functions for a single element
+    const auto& localFiniteElement = localView.tree().finiteElement();
+    const auto nSf = localFiniteElement.localBasis().size();
+    // Set all entries to zero
+    elementRhs.resize(localFiniteElement.size());
+    elementRhs = 0;
+    
+    // A quadrature rule
+//     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
+    int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+//     int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+    
+//     std::cout << "elementRhs.size():" << elementRhs.size() << std::endl;
+    
+    
+    
+    const auto& quadRule = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+    for (const auto& quadPoint : quadRule)
+    {
+        const FieldVector<double,dim>& quadPos = quadPoint.position();
+        
+        const double integrationElement = element.geometry().integrationElement(quadPos);
+//         double functionValue = forceTerm(element.geometry().global(quadPos));
+
+        // Evaluate all shape function values at this point
+//         std::vector<FieldVector<double,1> > shapeFunctionValues;
+//         localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+        
+        const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+        
+        std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+        localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+        
+        // Compute the shape function gradients on the grid element
+        std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+        for (size_t i=0; i<gradients.size(); i++)
+          jacobian.mv(referenceGradients[i][0], gradients[i]);
+        
+        
+        
+        // Actually compute the vector entries
+        for (size_t p=0; p<nSf; p++)
+        {
+            auto localIndex = localView.tree().localIndex(p);
+//             elementRhs[localIndex] += shapeFunctionValues[p] *  forceTerm(quadPos) * quadPoint.weight() * integrationElement;
+//             auto value = shapeFunctionValues[p]*  (sqrt(2.0)*mu(quadPos) *forceTerm(quadPos));
+//             auto value = -1.0*gradients[p][0] *  (sqrt(2.0)*mu(quadPos) *forceTerm(quadPos));             //NEGATIVE!!! TODO
+//             auto value = -2.0*mu(quadPos)*(sqrt(2.0)*forceTerm(quadPos))*gradients[p][0]  ;   
+//             auto value = -1.0*mu(quadPos)*forceTerm(quadPos)*gradients[p][0]  ;
+            
+            
+            
+            
+//             auto value = -2.0*sqrt(2.0)*mu(quadPos)*forceTerm(quadPos)*gradients[p][0];   //#1
+            
+            
+            auto value = 2.0*sqrt(2.0)*mu(quadPos)*forceTerm(quadPos)*gradients[p][0];  // TEST
+
+            
+            
+            
+            
+//             auto value = 2.0*mu(quadPos)*sqrt(2.0)*forceTerm(quadPos)*gradients[p][0]; 
+//             auto value = -2.0*mu(quadPos)*sqrt(2.0)*quadPos[1]*gradients[p][0];    //TEST
+//             std::cout << "value:" << value << std::endl;
+            
+//             std::cout << "forceTerm(quadPos):" << forceTerm(quadPos) << std::endl;
+//             std::cout << "quadPos:" << quadPos << std::endl;
+//             std::cout << "integrationElement:" << integrationElement << std::endl;
+            
+//             auto value = -1.0*sqrt(2.0)*forceTerm(quadPos)*gradients[p][0]  ;   //TEST
+            
+//             auto value = -1.0*mu(quadPos)*sqrt(2.0)*forceTerm(quadPos)*forceTerm(quadPos)*gradients[p][0]; //TEST
+            
+            elementRhs[localIndex] += value * quadPoint.weight() * integrationElement;
+        }
+    }
+}
+// Get the occupation pattern of the stiffness matrix
+template<class Basis>
+void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
+{
+    nb.resize(basis.size(), basis.size());
+    auto gridView = basis.gridView();
+    // A loop over all elements of the grid
+    auto localView = basis.localView();
+    for (const auto& element : elements(gridView))
+    {
+        localView.bind(element);
+        for (size_t i=0; i<localView.size(); i++)
+        {
+            // The global index of the i-th vertex of the element
+            auto row = localView.index(i);
+            for (size_t j=0; j<localView.size(); j++ )
+            {
+                // The global index of the j-th vertex of the element
+                auto col = localView.index(j);
+                nb.add(row,col);
+            }
+        }
+    }
+}
+/** \brief Assemble the Laplace stiffness matrix on the given grid view */
+template<class Basis, class Matrix, class Vector, class LocalFunction, class Force>
+void assemblePoissonProblem(const Basis& basis,
+                            Matrix& matrix,
+                            Vector& b,
+                            LocalFunction& muLocal,
+                            const Force& forceTerm,
+                            const double gamma)
+{
+    auto gridView = basis.gridView();
+
+
+    MatrixIndexSet occupationPattern;
+    getOccupationPattern(basis, occupationPattern);
+    occupationPattern.exportIdx(matrix);
+    matrix = 0;
+    
+    
+    auto loadGVF  = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
+    auto loadFunctional = localFunction(loadGVF);      
+
+    b.resize(basis.dimension());
+    b = 0;
+
+    auto localView = basis.localView();
+    
+    for (const auto& element : elements(gridView))
+    {
+
+        localView.bind(element);
+        muLocal.bind(element);
+        loadFunctional.bind(element);
+        
+        
+//         Dune::Matrix<double> elementMatrix;
+        Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+//         Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+        assembleElementStiffnessMatrix(localView, elementMatrix, muLocal, gamma);
+
+//         std::cout << "elementMatrix.N() "<<  elementMatrix.N() << std::endl;
+//         std::cout << "elementMatrix.M() "  << elementMatrix.M() << std::endl;
+        
+        
+        for(size_t p=0; p<elementMatrix.N(); p++)
+        {
+            auto row = localView.index(p);
+            for (size_t q=0; q<elementMatrix.M(); q++ )
+            {
+                auto col = localView.index(q);
+                matrix[row][col] += elementMatrix[p][q];
+            }
+        }
+
+        
+//         BlockVector<double> elementRhs;
+        BlockVector<FieldVector<double,1> > elementRhs;
+        assembleElementVolumeTerm(localView, elementRhs, muLocal, loadFunctional);
+        for (size_t p=0; p<elementRhs.size(); p++)
+        {
+            // The global index of the p-th vertex of the element
+            auto row = localView.index(p);
+            b[row] += elementRhs[p];
+        }
+    }
+}
+
+
+
+
+
+template<class Basis, class Vector, class LocalFunc1, class LocalFunc2>
+double computeMuGamma(const Basis& basis,
+                      Vector& coeffVector,
+                      const double gamma, 
+                      LocalFunc1& mu,
+                      LocalFunc2& Func
+                     )
+{
+
+
+//   constexpr int dim = Basis::LocalView::Element::dimension;
+  
+  double output = 0.0;
+  double Term1 = 0.0;
+  double Term2 = 0.0;
+  double Term11 = 0.0;
+  constexpr int dim = 2;
+//   constexpr int dim = 3;
+  auto localView = basis.localView();
+  
+//   auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<double>(basis, coeffVector);
+//   auto localSol = localFunction(solutionFunction);
+  
+//   auto loadGVF  = Dune::Functions::makeGridViewFunction(x3Fun, basis.gridView());
+//   auto x3Functional = localFunction(loadGVF);  
+  
+  
+  
+  
+  double area = 0.0;
+  
+  for(const auto& element : elements(basis.gridView()))
+  {
+      
+//         std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+        double elementEnergy1 = 0.0;
+        double elementEnergy2 = 0.0;
+        
+        localView.bind(element);
+        mu.bind(element);
+    //     x3Functional.bind(element);
+        Func.bind(element);
+        
+        auto geometry = element.geometry();
+        const auto& localFiniteElement = localView.tree().finiteElement();
+        const auto nSf = localFiniteElement.localBasis().size();
+
+//         int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
+        int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+    //     int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+        const auto& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);            // TODO WARUM HIER KEINE COLOR ?? ERROR
+        
+    //     std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+
+        for(const auto& quadPoint : quad)
+        {
+            const auto& quadPos = quadPoint.position();
+            //       const FieldVector<double,dim>& quadPos = quadPoint.position();
+            //       std::cout << " quadPOS: " << quadPos << std::endl;
+
+            
+            const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+            const double integrationElement = geometry.integrationElement(quadPos);
+            
+            area += quadPoint.weight() * integrationElement;
+            
+            std::vector< FieldMatrix<double, 1, dim>> referenceGradients;
+            localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+            
+            // Compute the shape function gradients on the grid element
+            std::vector<FieldVector<double,dim>> gradients(referenceGradients.size());
+            
+            for (size_t i=0; i<gradients.size(); i++)
+                jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+            
+            FieldVector<double,dim> tmp(0.0);
+        //       std::cout << "integrationElement :" << integrationElement << std::endl;
+        //       std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
+            
+            for (size_t i=0; i < nSf; i++)
+            {
+                size_t localIdx = localView.tree().localIndex(i);  // hier i:leafIdx
+                size_t globalIdx = localView.index(localIdx);
+        //         printvector(std::cout, gradients[i], "gradients[i]", "--" ); 
+        //         std::cout << "coeffVector[globalIdx]:" << coeffVector[globalIdx] << std::endl;
+                tmp[0] += coeffVector[globalIdx]*gradients[i][0];
+                tmp[1] += coeffVector[globalIdx]*gradients[i][1];
+        //         tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
+        //         tmp[1] += coeffVector[globalIdx]*gradients[i][2];
+        //         printvector(std::cout, tmp, "tmp", "--" ); 
+            }
+        //       printvector(std::cout, tmp, "gradient_w", "--" );
+            
+        //       auto value1 = std::pow( ((quadPos[1]/sqrt(2.0))+ (tmp[0]/2.0)) ,2); //TEST
+            
+            
+            
+            auto q1 = (Func(quadPos)/sqrt(2.0));   //#1
+            auto q2 = (tmp[0]/2.0);                //#1
+            
+            
+//             auto q2 = (tmp[0]/sqrt(2.0));  // TEST !!!!!!!!!!
+            
+            
+            
+            
+
+            // CHECK : BEITRÄGE CHECKEN!!!!
+            
+//             std::cout << "Beitrag1: " << q2 << std::endl;
+//             std::cout << "Beitrag2: " << (tmp[1]/(2.0*gamma))  << std::endl;
+//             
+//             
+//             
+            
+            auto q3 = q1 - q2;              // TEST
+            
+//             auto q3 = q1 + q2;              // #1      
+            auto value1 = std::pow(q3,2);
+
+            
+            
+//             auto value2 = std::pow( (tmp[1]/(2.0*gamma) )  , 2);
+            auto value2 = std::pow( (tmp[1]/(sqrt(2.0)*gamma) )  , 2);
+
+            
+            
+//             auto value = 2.0*mu(quadPos)*(2.0*value1 + value2);
+            
+
+            
+            elementEnergy1 +=  2.0*mu(quadPos)* 2.0*value1  * quadPoint.weight() * integrationElement;
+            elementEnergy2 +=  2.0*mu(quadPos)* value2  * quadPoint.weight() * integrationElement;
+        //       std::cout << "output:" << output << std::endl;
+        }
+//         std::cout << "elementEnergy:" << elementEnergy << std::endl;
+
+        Term1 += elementEnergy1;
+        Term2 += elementEnergy2;
+        
+//         std::cout << "output: " << output << std::endl;
+  } 
+  std::cout << "Term1: " << Term1 << std::endl;
+  std::cout << "Term2: " << Term2  << std::endl;
+  output = Term1 + Term2;
+  std::cout << "output: " << output << std::endl;
+  std::cout << "Domain-Area: " << area << std::endl;
+  return (1.0/area) * output;
+}
+
+
+// //  -----------------------------------------------------------
+/*
+template<class Basis, class Vector, class LocalFunc1, class LocalFunc2>
+double computeMuGamma(const Basis& basis,
+                      Vector& coeffVector,
+                      const double gamma, 
+                      LocalFunc1& mu,
+                      LocalFunc2& Func
+                     )
+{
+
+
+//   constexpr int dim = Basis::LocalView::Element::dimension;
+  
+  double output = 0.0;
+  constexpr int dim = 2;
+//   constexpr int dim = 3;
+  auto localView = basis.localView();
+  
+//   auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<double>(basis, coeffVector);
+//   auto localSol = localFunction(solutionFunction);
+  
+//   auto loadGVF  = Dune::Functions::makeGridViewFunction(x3Fun, basis.gridView());
+//   auto x3Functional = localFunction(loadGVF);  
+  
+  
+  
+  
+  double area = 0.0;
+  
+  for(const auto& element : elements(basis.gridView()))
+  {
+      
+//         std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+        double elementEnergy = 0.0;
+        
+        localView.bind(element);
+        mu.bind(element);
+    //     x3Functional.bind(element);
+        Func.bind(element);
+        
+        auto geometry = element.geometry();
+        const auto& localFiniteElement = localView.tree().finiteElement();
+        const auto nSf = localFiniteElement.localBasis().size();
+
+    //     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
+        int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+    //     int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+        const auto& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);            // TODO WARUM HIER KEINE COLOR ?? ERROR
+        
+    //     std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+
+        for(const auto& quadPoint : quad)
+        {
+            
+            const auto& quadPos = quadPoint.position();
+            //       const FieldVector<double,dim>& quadPos = quadPoint.position();
+            //       std::cout << " quadPOS: " << quadPos << std::endl;
+
+            
+            const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+            const double integrationElement = geometry.integrationElement(quadPos);
+            
+            area += quadPoint.weight() * integrationElement;
+            
+            std::vector< FieldMatrix<double, 1, dim>> referenceGradients;
+            localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+            
+            // Compute the shape function gradients on the grid element
+            std::vector<FieldVector<double,dim>> gradients(referenceGradients.size());
+            
+            for (size_t i=0; i<gradients.size(); i++)
+                jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+            
+//             FieldVector<double,dim> tmp = {0.0 , 0.0};
+            FieldVector<double,dim> tmp(0.0);
+        //       FieldVector<double,dim> tmp = {0.0 ,0.0,  0.0};  //3D-Version 
+            
+        //       std::cout << "integrationElement :" << integrationElement << std::endl;
+        //       std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
+            
+            for (size_t i=0; i < nSf; i++)
+            {
+                size_t localIdx = localView.tree().localIndex(i);  // hier i:leafIdx
+                size_t globalIdx = localView.index(localIdx);
+                
+        //         printvector(std::cout, gradients[i], "gradients[i]", "--" ); 
+        //         std::cout << "coeffVector[globalIdx]:" << coeffVector[globalIdx] << std::endl;
+                
+                tmp[0] += coeffVector[globalIdx]*gradients[i][0];
+                tmp[1] += coeffVector[globalIdx]*gradients[i][1];
+                
+        //         tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
+        //         tmp[1] += coeffVector[globalIdx]*gradients[i][2];
+        //         printvector(std::cout, tmp, "tmp", "--" ); 
+
+            }
+        //       printvector(std::cout, tmp, "gradient_w", "--" );
+            
+            
+        //       auto value1 = std::pow( ((quadPos[1]/sqrt(2.0))+ (tmp[0]/2.0)) ,2); //TEST
+            
+            
+            
+            auto q1 = (Func(quadPos)/sqrt(2.0));
+            auto q2 = (tmp[0]/2.0);
+            
+            
+//             auto q2 = (tmp[0]/sqrt(2.0));  // TEST !!!!!!!!!!
+            
+            
+            
+            
+
+            // CHECK : BEITRÄGE CHECKEN!!!!
+            
+            std::cout << "Beitrag1: " << q2 << std::endl;
+            std::cout << "Beitrag2: " << (tmp[1]/(2.0*gamma))  << std::endl;
+            
+            
+            
+            
+            auto q3 = q1 + q2;
+            auto value1 = std::pow(q3,2);
+//             auto value1 = std::pow( ((Func(quadPos)/sqrt(2.0)) + (tmp[0]/2.0)) , 2);
+            
+            
+            auto value2 = std::pow( (tmp[1]/(2.0*gamma) )  , 2);
+//             auto value2 = std::pow( (tmp[1]/(sqrt(2.0)*gamma) )  , 2);   //TEST 
+            
+            
+        //       auto value2 =  (1.0/(std::pow(gamma,2)))* std::pow(tmp[1],2)/4.0 ; //TEST
+            
+            auto value = 2.0*mu(quadPos)*(2.0*value1 + value2);
+            
+        //       std::cout << "quadPos[1]:" << quadPos[1]<< std::endl;
+        //       std::cout << "Func(quadPos):" << Func(quadPos) << std::endl;
+        //       std::cout << "sqrt(2.0):" << sqrt(2.0) << std::endl;
+        //       std::cout << "tmp[0]:" << tmp[0] << std::endl;
+        //       std::cout << "tmp[1]:" << tmp[1] << std::endl;
+        //       std::cout << "value1:" << value1 << std::endl;
+        //       std::cout << "value2:" << value2 << std::endl;
+        //       std::cout << "value:" << value << std::endl;
+            
+            
+        //       auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) );
+            
+            
+            
+        //       auto value = 2.0*mu(quadPos)*(2.0* (((x3Functional(quadPos)*x3Functional(quadPos))/2.0) + std::pow( (tmp[0]/2.0) ,2)) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) ); //TEST
+            
+        //       auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) ) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) ; //TEST
+        //       auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2)) + (1.0/gamma)*std::pow( (tmp[1]/2.0) ,2) ;  //TEST
+        //       auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2)  + (1.0/gamma)*std::pow( (tmp[1]/2.0) ,2) ) ;  //TEST
+            
+            elementEnergy +=  value * quadPoint.weight() * integrationElement;
+        //       std::cout << "output:" << output << std::endl;
+        }
+//         std::cout << "elementEnergy:" << elementEnergy << std::endl;
+        output += elementEnergy;
+//         std::cout << "output: " << output << std::endl;
+  }
+  std::cout << "Domain-Area: " << area << std::endl;
+  return (1.0/area) * output;
+}
+// -------------------------------------------------------------------------
+*/
+
+
+  // Check whether two points are equal on R/Z x R/Z
+  auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
+  {
+    return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+           and (FloatCmp::eq(x[1],y[1])) );
+  };
+
+
+
+int main(int argc, char *argv[])
+{
+
+  MPIHelper::instance(argc, argv);
+
+    
+    
+  ParameterTree parameterSet;
+  if (argc < 2)
+    ParameterTreeParser::readINITree("../../inputs/computeMuGamma.parset", parameterSet);
+  else
+  {
+    ParameterTreeParser::readINITree(argv[1], parameterSet);
+    ParameterTreeParser::readOptions(argc, argv, parameterSet);
+  }
+    /////////////////////////////////
+    // SET OUTPUT 
+    /////////////////////////////////
+    std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs");
+    std::fstream log;
+    log.open(outputPath + "/outputMuGamma.txt" ,std::ios::out);
+    std::cout << "outputPath:" << outputPath << std::endl;
+  
+  
+    //////////////////////////////////
+    // Generate the grid
+    //////////////////////////////////
+    constexpr int dim = 2;
+    
+    // QUAD-GRID
+    FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0});
+    FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0});
+//     std::array<int, dim> nElements = {16,16};
+    
+    std::array<int,2> nElements = parameterSet.get<std::array<int,2>>("nElements", {32,32});
+    std::cout << "Number of Elements in each direction: " << nElements << std::endl;
+    log << "Number of Elements in each direction: " << nElements << std::endl;
+    
+    
+    using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+    CellGridType grid_CE(lower,upper,nElements);
+    using GridView = CellGridType::LeafGridView;
+    const GridView gridView = grid_CE.leafGridView();
+    using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+    
+    // ----------- INPUT PARAMETERS -----------------------------
+    std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
+    log << "material_prestrain used: "<< imp  << std::endl;
+    double gamma   = parameterSet.get<double>("gamma", 1.0);   
+    double theta   = parameterSet.get<double>("theta", 1.0/4.0); 
+    double mu1     = parameterSet.get<double>("mu1", 10.0);
+    double beta    = parameterSet.get<double>("beta", 2.0); 
+    double mu2     = beta*mu1;
+    std::cout << "Gamma:" << gamma << std::endl;
+    std::cout << "Theta:" << theta  << std::endl;
+    std::cout << "mu1:" << mu1 << std::endl;
+    std::cout << "mu2:" << mu2 << std::endl;
+    std::cout << "beta:" << beta  << std::endl;
+    log << "----- Input Parameters -----: " << std::endl;
+    log << "gamma: " << gamma << std::endl;
+    log << "theta: " << theta << std::endl;
+    log << "beta: " << beta << std::endl;
+    log << "material parameters: " << std::endl;
+    log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+
+//     auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+// 
+// //             if (abs(z[0]) > (theta/2.0)) 
+//             if (abs(z[0]) >= (theta/2.0))   
+//                 return mu1;
+//             else
+//                 return mu2;
+//             };
+
+    auto materialImp = IsotropicMaterialImp<dim>();
+    auto muTerm = materialImp.getMu(parameterSet);
+    auto muGridF  = Dune::Functions::makeGridViewFunction(muTerm, gridView);
+    auto muLocal  = localFunction(muGridF);
+    
+
+    /////////////////////////////////////////////////////////
+    // Stiffness matrix and right hand side vector
+    /////////////////////////////////////////////////////////
+    using Vector = BlockVector<FieldVector<double,1> >; 
+    using Matrix = BCRSMatrix<FieldMatrix<double,1,1> >;
+    Matrix stiffnessMatrix;
+    Vector b;
+
+    /////////////////////////////////////////////////////////
+    // Assemble the system
+    /////////////////////////////////////////////////////////
+    using namespace Functions::BasisFactory;
+    Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+
+    // Don't do the following in real life: It has quadratic run-time in the number of vertices.
+    for (const auto& v1 : vertices(gridView))
+      for (const auto& v2 : vertices(gridView))
+        if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))                                      
+          periodicIndices.unifyIndexPair({gridView.indexSet().index(v1)}, {gridView.indexSet().index(v2)});
+
+    auto basis = makeBasis(gridView, Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices)); // flatLexicographic()?
+
+
+    auto forceTerm = [](const FieldVector<double,dim>& x){return x[1];}; //2D-Version
+    auto ForceGridF  = Dune::Functions::makeGridViewFunction(forceTerm, gridView);
+    auto ForceLocal  = localFunction(ForceGridF);
+    
+    assemblePoissonProblem(basis, stiffnessMatrix, b, muLocal, forceTerm, gamma);
+//     printmatrix(std::cout, stiffnessMatrix, "StiffnessMatrix", "--");
+//     printvector(std::cout, b, "b", "--");
+
+
+
+
+
+
+    ///////////////////////////
+    // Compute solution
+    ///////////////////////////
+
+    Vector x(basis.size());
+    x = b;
+    
+     std::cout << "------------ CG - Solver ------------" << std::endl;
+    MatrixAdapter<Matrix, Vector, Vector> op(stiffnessMatrix);
+
+
+
+    // Sequential incomplete LU decomposition as the preconditioner
+    SeqILU<Matrix, Vector, Vector> ilu0(stiffnessMatrix,1.0);
+    int iter = parameterSet.get<double>("iterations_CG", 999);
+    // Preconditioned conjugate-gradient solver
+    CGSolver<Vector> solver(op,
+                            ilu0, //NULL,
+                            1e-8, // desired residual reduction factorlack
+                            iter,   // maximum number of iterations
+                            2);   // verbosity of the solver
+    InverseOperatorResult statistics;
+    // Solve!
+    solver.apply(x, b, statistics);
+    
+//     std::cout << "------------ GMRES - Solver ------------" << std::endl;
+//     // Turn the matrix into a linear operator
+//     MatrixAdapter<Matrix, Vector, Vector> op(stiffnessMatrix);
+// 
+//     // Fancy (but only) way to not have a preconditioner at all
+//     Richardson<Vector,Vector> preconditioner(1.0);
+// 
+//     // Construct the iterative solver
+//     RestartedGMResSolver<Vector> solver(
+//         op, // Operator to invert
+//         preconditioner,    // Preconditioner
+//         1e-10,             // Desired residual reduction factor
+//         500,               // Number of iterations between restarts,
+//                         // here: no restarting
+//         500,               // Maximum number of iterations
+//         2);                // Verbosity of the solver
+// 
+//     // Object storing some statistics about the solving process
+//     InverseOperatorResult statistics;
+// 
+//     // solve for different Correctors (alpha = 1,2,3)
+//     solver.apply(x, b, statistics);
+//     
+    
+// -----------------------------------------------------------------------------------------------------
+    
+    
+    using SolutionRange = double;
+    auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(basis, x);
+    
+    
+    //  -------- PRINT OUTPUT --------
+//     printvector(std::cout, x, "coeffVector", "--" );
+    std::cout << "Gamma:" << gamma << std::endl;
+    double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, ForceLocal);
+    std::cout << "mu_gamma:" << mu_gamma << std::endl;
+    log << "----- OUTPUT: -----: " << std::endl;
+    log << "mu_gamma=" << mu_gamma << std::endl;   
+    
+//     std::cout.precision(10);
+//     std::cout << "mu_gamma:" << std::fixed << mu_gamma  << std::endl;
+    
+    
+
+//     Vector zeroVec(basis.size());
+//     zeroVec = 0;
+//     std::cout << "TEST computeMuGamma:" << computeMuGamma(basis, zeroVec, gamma, muLocal, ForceLocal)<< std::endl;
+
+    std::cout << " --- print analytic solutions --- " << std::endl;
+    
+    double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+    double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
+    double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);  
+    double am = mu1*((1-theta)+theta*beta)     *(1.0/6.0);
+    
+    std::cout << "q1 : "     << q1 << std::endl;
+    std::cout << "q2 : "     << q2 << std::endl;
+    std::cout << "q3 should be between q1 and q2"  << std::endl;
+    std::cout << "hm : "     << hm << std::endl;
+    std::cout << "am : "     << am << std::endl;
+    
+    
+    if(mu_gamma > q2)
+    {
+        std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
+    }
+
+    
+    std::cout << "beta : "     << beta << std::endl;
+    std::cout << "theta : "     << theta << std::endl;
+    std::cout << "Gamma : "     << gamma << std::endl;
+    std::cout << "mu_gamma:" << mu_gamma << std::endl;
+    
+    
+
+    // Output result
+
+    std::string VTKOutputName =  outputPath + "/Compute_MuGamma-Result";
+    
+    VTKWriter<GridView> vtkWriter(gridView);
+//     vtkWriter.addVertexData(x, "solution");                      //--- Anpassen für P2
+    vtkWriter.addVertexData(
+        solutionFunction,
+        VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
+//         VTK::FieldInfo("solution", VTK::FieldInfo::Type::vector, dim));
+
+    vtkWriter.write( VTKOutputName );
+    std::cout << "wrote data to file: " + VTKOutputName  << std::endl;     
+    
+    
+    
+    using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+    VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0},{64,64});
+    using GridViewVTK = VTKGridType::LeafGridView;
+    const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+    
+    auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+    
+    std::vector<double> mu_CoeffP0;
+    Functions::interpolate(scalarP0FeBasis, mu_CoeffP0, muTerm);
+    auto mu_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, mu_CoeffP0);
+        
+     VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+
+       MaterialVtkWriter.addCellData(
+            mu_DGBF_P0,
+            VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));  
+     
+     MaterialVtkWriter.write(outputPath + "/MaterialFunctions" );
+        std::cout << "wrote data to file:" + outputPath  +  "/MaterialFunctions" << std::endl;  
+     
+}

src/HelperFunctions.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import time
+from ClassifyMin import *
+# from scipy.io import loadmat #not Needed anymore?
+import codecs
+import sys
+
+
+def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
+    # Read Output Matrices (effective quantities)
+    # From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
+    # -- Read Matrix Qhom
+    X = []
+    # with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
+    with codecs.open(QFilePath, encoding='utf-8-sig') as f:
+        for line in f:
+            s = line.split()
+            X.append([float(s[i]) for i in range(len(s))])
+    Q = np.array([[X[0][2], X[1][2], X[2][2]],
+                  [X[3][2], X[4][2], X[5][2]],
+                  [X[6][2], X[7][2], X[8][2]] ])
+
+    # -- Read Beff (as Vector)
+    X = []
+    # with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
+    with codecs.open(BFilePath, encoding='utf-8-sig') as f:
+        for line in f:
+            s = line.split()
+            X.append([float(s[i]) for i in range(len(s))])
+    B = np.array([X[0][2], X[1][2], X[2][2]])
+    return Q, B
+
+
+
+def SetParametersCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+    with open(InputFilePath, 'r') as file:
+        filedata = file.read()
+    filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+    filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+    filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+    filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+    filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+    filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+    f = open(InputFilePath,'w')
+    f.write(filedata)
+    f.close()
+
+#TODO combine these...
+
+def SetParametersComputeMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+    with open(InputFilePath, 'r') as file:
+        filedata = file.read()
+    filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+    filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+    filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+    filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+    filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+    f = open(InputFilePath,'w')
+    f.write(filedata)
+    f.close()
+
+
+
+
+def RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+        # with open(InputFilePath, 'r') as file:
+        #     filedata = file.read()
+        # filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+        # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+        # filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+        # filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+        # filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+        # filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+        # f = open(InputFilePath,'w')
+        # f.write(filedata)
+        # f.close()
+        SetParametersCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath)
+        # --- Run Cell-Problem
+        # Optional: Check Time
+        # t = time.time()
+        subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+                                             capture_output=True, text=True)
+        # print('elapsed time:', time.time() - t)
+        # --------------------------------------------------------------------------------------
+
+
+
+
+
+
+
+
+
+def GetCellOutput(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+                OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+        RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath)
+        print('Read effective quantities...')
+        Q, B = ReadEffectiveQuantities()
+        # print('Q:', Q)
+        # print('B:', B)
+        return Q, B
+
+
+# unabhängig von alpha...
+def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+                OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+    # ------------------------------------ get mu_gamma ------------------------------
+    # ---Scenario 1.1: extreme regimes
+    if gamma == '0':
+        # print('extreme regime: gamma = 0')
+        mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
+        # print("mu_gamma:", mu_gamma)
+    elif gamma == 'infinity':
+        # print('extreme regime: gamma = infinity')
+        mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta)   # = q1
+        # print("mu_gamma:", mu_gamma)
+    else:
+        # --- Scenario 1.2:  compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
+        # print("Run computeMuGamma for Gamma = ", gamma)
+
+        SetParametersComputeMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath)
+        # with open(InputFilePath, 'r') as file:
+        #     filedata = file.read()
+        # filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+        # # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+        # filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+        # filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+        # filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+        # filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+        # f = open(InputFilePath,'w')
+        # f.write(filedata)
+        # f.close()
+        # --- Run Cell-Problem
+
+        # Check Time
+        # t = time.time()
+        # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+        #                                      capture_output=True, text=True)
+        # --- Run Cell-Problem_muGama   -> faster
+        # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+        #                                              capture_output=True, text=True)
+        # --- Run Compute_muGamma (2D Problem much much faster)
+
+        subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+                                                             capture_output=True, text=True)
+        # print('elapsed time:', time.time() - t)
+
+        #Extract mu_gamma from Output-File                                           TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+        with open(OutputFilePath, 'r') as file:
+            output = file.read()
+        tmp = re.search(r'(?m)^mu_gamma=.*',output).group()                           # Not necessary for Intention of Program t output Minimizer etc.....
+        s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+        mu_gamma = float(s[0])
+        # print("mu_gamma:", mu_gammaValue)
+    # --------------------------------------------------------------------------------------
+    return mu_gamma
+
+
+
+
+
+
+def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+        # ---------------------------------------------------------------
+        # Comparison of the analytical Classification 'ClassifyMin'
+        # and the symbolic Minimizatio + Classification 'symMinimization'
+        # ----------------------------------------------------------------
+        comparison_successful = True
+        eps = 1e-8
+
+        # 1. Substitute Input-Parameters for the Cell-Problem
+        with open(InputFilePath, 'r') as file:
+            filedata = file.read()
+        filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+        filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+        filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+        filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+        filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+        filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+        f = open(InputFilePath,'w')
+        f.write(filedata)
+        f.close()
+        # 2. --- Run Cell-Problem
+        print('Run Cell-Problem ...')
+        # Optional: Check Time
+        # t = time.time()
+        subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+                                             capture_output=True, text=True)
+
+
+        # 3. --- Run Matlab symbolic minimization program: 'symMinimization'
+        eng = matlab.engine.start_matlab()
+        # s = eng.genpath(path + '/Matlab-Programs')
+        s = eng.genpath(os.path.dirname(os.getcwd()))
+        eng.addpath(s, nargout=0)
+        # print('current Matlab folder:', eng.pwd(nargout=1))
+        eng.cd('Matlab-Programs', nargout=0)  #switch to Matlab-Programs folder
+        # print('current Matlab folder:', eng.pwd(nargout=1))
+        Inp = False
+        print('Run symbolic Minimization...')
+        G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4)  #Name of program:symMinimization
+        # G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4)  #Optional: add Path
+        G = np.asarray(G) #cast Matlab Outout to numpy array
+        # --- print Output ---
+        print('Minimizer G:')
+        print(G)
+        print('angle:', angle)
+        print('type:', type )
+        print('curvature:', kappa)
+
+        # 4. --- Read the effective quantities (output values of the Cell-Problem)
+        # Read Output Matrices (effective quantities)
+        print('Read effective quantities...')
+        Q, B = ReadEffectiveQuantities()
+        print('Q:', Q)
+        print('B:', B)
+        q1 = Q[0][0]
+        q2 = Q[1][1]
+        q3 = Q[2][2]
+        q12 = Q[0][1]
+        b1 = B[0]
+        b2 = B[1]
+        b3 = B[2]
+        print("q1:", q1)
+        print("q2:", q2)
+        print("q3:", q3)
+        print("q12:", q12)
+        print("b1:", b1)
+        print("b2:", b2)
+        print("b3:", b3)
+
+        # --- Check Assumptions:
+        # Assumption of Classification-Lemma1.6:  [b3 == 0] & [Q orthotropic]
+        # Check if b3 is close to zero
+        assert (b3 < eps), "ClassifyMin only defined for b3 == 0 !"
+
+        # Check if Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
+        assert(Q[2][0] < eps and Q[0][2] < eps and Q[1][2] < eps and Q[2][1] < eps), "Q is not orthotropic !"
+
+        # 5. --- Get output from the analytical Classification 'ClassifyMin'
+        G_ana, angle_ana, type_ana, kappa_ana = classifyMin(q1, q2, q3, q12, b1, b2)
+        print('Minimizer G_ana:')
+        print(G_ana)
+        print('angle_ana:', angle_ana)
+        print('type_ana:', type_ana )
+        print('curvature_ana:', kappa_ana)
+
+        # 6. Compare
+        # print('DifferenceMatrix:', G_ana - G )
+        # print('MinimizerError (FrobeniusNorm):', np.linalg.norm(G_ana - G , 'fro') )
+
+        if np.linalg.norm(G_ana - G , 'fro') > eps :
+            comparison_successful = False
+            print('Difference between Minimizers is too large !')
+        if type != type_ana :
+            comparison_successful = False
+            print('Difference in type !')
+        if abs(angle-angle_ana) > eps :
+            comparison_successful = False
+            print('Difference in angle is too large!')
+        if abs(kappa-kappa_ana) > eps :
+            comparison_successful = False
+            print('Difference in curvature is too large!')
+
+
+        if comparison_successful:
+            print('Comparison successful')
+        else:
+            print('Comparison unsuccessful')
+
+        return comparison_successful
+
+
+
+
+# ----------------------------------------------------------------------------------------------------------------------------

src/PhaseDiagram.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+
+import time
+# print(sys.executable)
+
+# --------------------------------------------------------------------
+# START :
+# INPUT (Parameters):   alpha, beta, theta, gamma, mu1, rho1
+#
+# -Option 1 : (Case lambda = 0 => q12 = 0)
+#   compute q1,q2,b1,b2 from Formula
+#       Option 1.1 :
+#           set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
+#       Option 1.2 :
+#           compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
+# -Option 2 :
+#   compute Q_hom & B_eff by running 'Cell-Problem'
+#
+# -> CLASSIFY ...
+#
+# OUTPUT: Minimizer G, angle , type, curvature
+# -----------------------------------------------------------------------
+#
+#
+# def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset",
+#                 OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ):
+#     # ------------------------------------ get mu_gamma ------------------------------
+#     # ---Scenario 1.1: extreme regimes
+#     if gamma == '0':
+#         print('extreme regime: gamma = 0')
+#         mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
+#         print("mu_gamma:", mu_gamma)
+#     elif gamma == 'infinity':
+#         print('extreme regime: gamma = infinity')
+#         mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta)   # = q1
+#         print("mu_gamma:", mu_gamma)
+#     else:
+#         # --- Scenario 1.2:  compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
+#         # print("Run computeMuGamma for Gamma = ", gamma)
+#         with open(InputFilePath, 'r') as file:
+#             filedata = file.read()
+#         filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+#         # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+#         filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+#         filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+#         filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+#         filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+#         f = open(InputFilePath,'w')
+#         f.write(filedata)
+#         f.close()
+#         # --- Run Cell-Problem
+#
+#         # Check Time
+#         # t = time.time()
+#         # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+#         #                                      capture_output=True, text=True)
+#         # --- Run Cell-Problem_muGama   -> faster
+#         # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+#         #                                              capture_output=True, text=True)
+#         # --- Run Compute_muGamma (2D Problem much much faster)
+#
+#         subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+#                                                              capture_output=True, text=True)
+#         # print('elapsed time:', time.time() - t)
+#
+#         #Extract mu_gamma from Output-File                                           TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+#         with open(OutputFilePath, 'r') as file:
+#             output = file.read()
+#         tmp = re.search(r'(?m)^mu_gamma=.*',output).group()                           # Not necessary for Intention of Program t output Minimizer etc.....
+#         s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+#         mu_gamma = float(s[0])
+#         # print("mu_gamma:", mu_gammaValue)
+#     # --------------------------------------------------------------------------------------
+#     return mu_gamma
+#
+
+
+
+# ----------- SETUP PATHS
+# InputFile  = "/inputs/cellsolver.parset"
+# OutputFile = "/outputs/output.txt"
+InputFile  = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run  from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+
+# -------------------------- Input Parameters --------------------
+mu1 = 10.0               # TODO : here must be the same values as in the Parset for computeMuGamma
+rho1 = 1.0
+alpha = 2.0
+beta = 2.0
+# beta = 10.0
+theta = 1.0/4.0
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+gamma = 0.25
+# gamma = 1.0
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+
+# gamma_min = 0.05
+# gamma_max = 2.0
+
+gamma_min = 1
+gamma_max = 1
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=1)
+
+# Gamma_Values = np.linspace(gamma_min, gamma_max, num=13)    # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+
+for gamma in Gamma_Values:
+
+
+    # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
+    # # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
+    # print('Test MuGamma:', muGamma)
+
+    # ------- Options --------
+    # print_Cases = True
+    # print_Output = True
+
+                            #TODO
+    generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
+
+    # make_3D_plot = True
+    # make_3D_PhaseDiagram = True
+    make_2D_plot = False
+    # make_2D_PhaseDiagram = False
+    make_3D_plot = False
+    make_3D_PhaseDiagram = False
+    # make_2D_plot = True
+    make_2D_PhaseDiagram = True
+
+
+    # --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
+    # q1 = harmonicMean(mu1, beta, theta)
+    # q2 = arithmeticMean(mu1, beta, theta)
+    # --- Set q12
+    # q12 = 0.0  # (analytical example)              # TEST / TODO read from Cell-Output
+
+
+
+
+
+    # b1 = prestrain_b1(rho1, beta, alpha, theta)
+    # b2 = prestrain_b2(rho1, beta, alpha, theta)
+    #
+    # print('---- Input parameters: -----')
+    # print('mu1: ', mu1)
+    # print('rho1: ', rho1)
+    # print('alpha: ', alpha)
+    # print('beta: ', beta)
+    # print('theta: ', theta)
+    # print("q1: ", q1)
+    # print("q2: ", q2)
+    # print("mu_gamma: ", mu_gamma)
+    # print("q12: ", q12)
+    # print("b1: ", b1)
+    # print("b2: ", b2)
+    # print('----------------------------')
+    # print("machine epsilon", sys.float_info.epsilon)
+
+    # G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12,  b1, b2, print_Cases, print_Output)
+    # Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+    # print("Test", Test)
+
+    # ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------
+
+    # SamplePoints_3D = 10 # Number of sample points in each direction
+    # SamplePoints_2D = 10 # Number of sample points in each direction
+    SamplePoints_3D = 300 # Number of sample points in each direction
+    SamplePoints_2D = 5 # Number of sample points in each direction
+
+
+
+    if make_3D_PhaseDiagram:
+        alphas_ = np.linspace(-20, 20, SamplePoints_3D)
+        betas_  = np.linspace(0.01,40.01,SamplePoints_3D)
+        thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
+        # print('type of alphas', type(alphas_))
+        # print('Test:', type(np.array([mu_gamma])) )
+        alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+        classifyMin_anaVec = np.vectorize(classifyMin_ana)
+
+        # Get MuGamma values ...
+        GetMuGammaVec = np.vectorize(GetMuGamma)
+        muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
+        # Classify Minimizers....
+        G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas,  mu1, rho1)   # Sets q12 to zero!!!
+        # print('size of G:', G.shape)
+        # print('G:', G)
+        # Out = classifyMin_anaVec(alphas,betas,thetas)
+        # T = Out[2]
+        # --- Write to VTK
+        GammaString = str(gamma)
+
+        VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString
+        gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+        print('Written to VTK-File:', VTKOutputName )
+
+    if make_2D_PhaseDiagram:
+        # alphas_ = np.linspace(-20, 20, SamplePoints_2D)
+        # thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
+        alphas_ = np.linspace(9, 10, SamplePoints_2D)
+        thetas_ = np.linspace(0.075,0.14,SamplePoints_2D)
+
+            # alphas_ = np.linspace(8, 12, SamplePoints_2D)
+            # thetas_ = np.linspace(0.05,0.2,SamplePoints_2D)
+        # betas_  = np.linspace(0.01,40.01,1)
+        #fix to one value:
+        betas_ = 2.0;
+        alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
+
+
+        if generalCase:              #TODO
+            classifyMin_matVec = np.vectorize(classifyMin_mat)
+            GetCellOutputVec = np.vectorize(GetCellOutput)
+            Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+
+
+            print('type of Q:', type(Q))
+            print('Q:', Q)
+            G, angles, Types, curvature = classifyMin_matVec(Q,B)
+
+        else:
+            classifyMin_anaVec = np.vectorize(classifyMin_ana)
+            GetMuGammaVec = np.vectorize(GetMuGamma)
+            muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
+            G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas,  mu1, rho1)    # Sets q12 to zero!!!
+            # print('size of G:', G.shape)
+            # print('G:', G)
+            # print('Types:', Types)
+            # Out = classifyMin_anaVec(alphas,betas,thetas)
+            # T = Out[2]
+            # --- Write to VTK
+            # VTKOutputName = + path + "./PhaseDiagram2DNEW"
+            GammaString = str(gamma)
+
+        VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
+        gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
+        print('Written to VTK-File:', VTKOutputName )
+
+
+    # --- Make 3D Scatter plot
+    if(make_3D_plot or make_2D_plot):
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+        colors = cm.plasma(Types)
+        # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
+        # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
+        if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
+        if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
+        # cbar=plt.colorbar(pnt3d)
+        # cbar.set_label("Values (units)")
+        ax.set_xlabel('alpha')
+        ax.set_ylabel('beta')
+        if make_3D_plot: ax.set_zlabel('theta')
+        plt.show()
+        # plt.savefig('common_labels.png', dpi=300)
+        # print('T:', T)
+        # print('Type 1 occured here:', np.where(T == 1))
+        # print('Type 2 occured here:', np.where(T == 2))
+
+
+# print(alphas_)
+# print(betas_)
+
+
+
+
+
+# ALTERNATIVE
+# colors = ("red", "green", "blue")
+# groups = ("Type 1", "Type2", "Type3")
+#
+# # Create plot
+# fig = plt.figure()
+# ax = fig.add_subplot(1, 1, 1)
+#
+# for data, color, group in zip(Types, colors, groups):
+#     # x, y = data
+#     ax.scatter(alphas, thetas, alpha=0.8, c=color, edgecolors='none', label=group)
+#
+# plt.title('Matplot scatter plot')
+# plt.legend(loc=2)
+# plt.show()

src/Plotq3-Angle.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+import time
+
+
+def find_nearest(array, value):
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return array[idx]
+
+
+def find_nearestIdx(array, value):
+    array = np.asarray(array)
+    idx = (np.abs(array - value)).argmin()
+    return idx
+
+
+InputFile  = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- Run  from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+print('---- Input parameters: -----')
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10
+# beta = 40.0
+# theta = 0.125
+#
+
+
+mu1 = 10.0
+rho1 = 1.0
+# alpha = 10.02021333
+alpha = 10.0
+beta = 2.0
+theta = 0.125
+# theta = 0.124242
+# gamma = 0.75
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+# gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+# gamma = 0.25
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+
+# ----------------------------------------------------------------
+
+
+gamma_min = 0.01
+gamma_max = 1.0
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=100)    # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+
+# # --- Options
+# # make_Plot = False
+make_Plot = True
+
+# Get values for mu_Gamma
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1, InputFilePath ,OutputFilePath )
+print('muGammas:', muGammas)
+
+q12 = 0.0
+q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+print('q1: ', q1)
+print('q2: ', q2)
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
+q3_star = math.sqrt(q1*q2)
+print('q3_star:', q3_star)
+
+# TODO these have to be compatible with input parameters!!!
+# compute certain ParameterValues that this makes sense
+# b1 = q3_star
+# b2 = q1
+print('b1: ', b1)
+print('b2: ', b2)
+
+# return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
+
+
+
+# classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+classifyMin_anaVec = np.vectorize(classifyMin_ana)
+G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+# _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+print('angles:', angles)
+
+idx = find_nearestIdx(muGammas, q3_star)
+print('GammaValue Idx closest to q_3^*', idx)
+gammaClose = Gamma_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose)
+
+
+
+
+
+
+# Make Plot
+if make_Plot:
+    plt.figure()
+    # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+    plt.title(r'angle$-\gamma$-Plot')
+    plt.plot(Gamma_Values, angles)
+    plt.scatter(Gamma_Values, angles)
+    # plt.plot(muGammas, angles)
+    # plt.scatter(muGammas, angles)
+    # plt.axis([0, 6, 0, 20])
+    # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+    # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+    # plt.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+    # plt.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+    # plt.axvline(x = q3_star, color = 'r', linestyle = 'dashed', label='$\gamma^*$')
+
+
+    # Plot Gamma Value that is closest to q3_star
+    plt.axvline(x = gammaClose, color = 'b', linestyle = 'dashed', label='$\gamma^*$')
+
+    plt.axvspan(gamma_min, gammaClose, color='red', alpha=0.5)
+    plt.axvspan(gammaClose, gamma_max, color='green', alpha=0.5)
+
+    plt.xlabel("$\gamma$")
+    plt.ylabel("angle")
+    plt.legend()
+    plt.show()

src/makePlot.py

+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+import math
+import os
+import subprocess
+import fileinput
+import re
+import matlab.engine
+from HelperFunctions import *
+# from subprocess import Popen, PIPE
+#import sys
+
+###################### makePlot.py #########################
+#  Generalized Plot-Script giving the option to define
+#  quantity of interest and the parameter it depends on
+#  to create a plot
+#
+#  Input: Define y & x for "x-y plot" as Strings
+#  - Run the 'Cell-Problem' for the different Parameter-Points
+#  (alternatively run 'Compute_MuGamma' if quantity of interest
+#   is q3=muGamma for a significant Speedup)
+
+###########################################################
+
+# TODO
+# - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3
+# - Also Add option to plot Minimization Output
+
+
+# ----- Setup Paths -----
+InputFile  = "/inputs/cellsolver.parset"
+OutputFile = "/outputs/output.txt"
+# path = os.getcwd()
+# InputFilePath = os.getcwd()+InputFile
+# OutputFilePath = os.getcwd()+OutputFile
+# --------- Run  from src folder:
+path_parent = os.path.dirname(os.getcwd())
+os.chdir(path_parent)
+path = os.getcwd()
+print(path)
+InputFilePath = os.getcwd()+InputFile
+OutputFilePath = os.getcwd()+OutputFile
+print("InputFilepath: ", InputFilePath)
+print("OutputFilepath: ", OutputFilePath)
+print("Path: ", path)
+
+#---------------------------------------------------------------
+
+print('---- Input parameters: -----')
+mu1 = 10.0
+lambda1 = 10.0
+rho1 = 1.0
+alpha = 2.8
+beta = 2.0
+theta = 1.0/4.0
+
+
+gamma = 1.0/4.0
+
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+print('gamma:', gamma)
+print('----------------------------')
+
+
+# TODO? : Ask User for Input ...
+# function = input("Enter value you want to plot (y-value):\n")
+# print(f'You entered {function}')
+# parameter = input("Enter Parameter this value depends on (x-value) :\n")
+# print(f'You entered {parameter}')
+
+# Add Option to change NumberOfElements used for computation of Cell-Problem
+
+
+# --- Define Quantity of interest:
+# Options: 'q1', 'q2', 'q3', 'q12' ,'q21', 'q31', 'q13' , 'q23', 'q32' , 'b1', 'b2' ,'b3'
+# TODO: EXTRA (MInimization Output) 'Minimizer (norm?)' 'angle', 'type', 'curvature'
+yName = 'q12'
+# yName = 'b1'
+yName = 'q3'
+yName = 'angle'
+
+# --- Define Parameter this function/quantity depends on:
+# Options: mu1 ,lambda1, rho1 , alpha, beta, theta, gamma
+# xName = 'theta'
+# xName = 'gamma'
+# xName = 'lambda1'
+xName = 'alpha'
+
+
+# --- define Interval of x-values:
+xmin = 0.0
+xmax = 10.0
+# xmin = 0.01
+# xmax = 3.0
+numPoints = 5
+X_Values = np.linspace(xmin, xmax, num=numPoints)
+print(X_Values)
+
+
+Y_Values = []
+
+
+
+# --- Options
+RUN = True
+# RUN = False
+# make_Plot = False
+make_Plot = True
+
+if RUN:
+    for x in X_Values:
+
+        if yName =='q3' or yName == 'mu_gamma':   # if only q3 is needed .. compute 2D problem (much faster)
+
+            # fist replace old parameters in parset
+            SetParametersComputeMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath)
+
+            # replace the sought after x value in the parset
+            with open(path+"/inputs/computeMuGamma.parset", 'r') as file:
+                filedata = file.read()
+            filedata = re.sub('(?m)^'+xName+'=.*',xName+'='+str(x),filedata)
+            f = open(path+"/inputs/computeMuGamma.parset",'w')
+            f.write(filedata)
+            f.close()
+            #Run Compute_MuGamma since its much faster
+            print("Run Compute_MuGamma for " + xName + " =" , x)
+            subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+                                                                 capture_output=True, text=True)
+            #Extract mu_gamma from Output-File
+            with open(path + '/outputs/outputMuGamma.txt', 'r') as file:
+                output = file.read()
+            tmp = re.search(r'(?m)^mu_gamma=.*',output).group()                           # Not necessary for Intention of Program t output Minimizer etc.....
+            s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+            Y_Values.append(float(s[0]))
+        else :     # Run full Cell-Problem...
+            # fist replace old parameters in parset
+            SetParametersCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
+
+            with open(InputFilePath, 'r') as file:
+                filedata = file.read()
+            filedata = re.sub('(?m)^'+xName+'=.*',xName+'='+str(x),filedata)
+            f = open(InputFilePath,'w')
+            f.write(filedata)
+            f.close()
+            # Run Cell-Problem
+            print("Run Cell-Problem for " + xName + " =" , x)
+            subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+                                                 capture_output=True, text=True)
+            #Extract quantity from Output-File
+            Q, B = ReadEffectiveQuantities()
+            if yName == 'q1':                   # TODO: Better use dictionary?...
+                print('q1 used')
+                Y_Values.append(Q[0][0])
+            elif yName =='q2':
+                print('q2 used')
+                Y_Values.append(Q[1][1])
+            # elif yName =='q3':
+            #     print('q3 used')
+            #     Y_Values.append(Q[0][0])
+            elif yName =='q12':
+                print('q12 used')
+                Y_Values.append(Q[0][1])
+            elif yName =='q21':
+                print('q21 used')
+                Y_Values.append(Q[1][0])
+            elif yName =='q13':
+                print('q13 used')
+                Y_Values.append(Q[0][2])
+            elif yName =='q31':
+                print('q31 used')
+                Y_Values.append(Q[2][0])
+            elif yName =='q23':
+                print('q23 used')
+                Y_Values.append(Q[1][2])
+            elif yName =='q32':
+                print('q32 used')
+                Y_Values.append(Q[2][1])
+            elif yName =='b1':
+                print('b1 used')
+                Y_Values.append(B[0])
+            elif yName =='b2':
+                print('b2 used')
+                Y_Values.append(B[1])
+            elif yName =='b3':
+                print('b3 used')
+                Y_Values.append(B[2])
+            elif yName == 'angle' or yName =='type' or yName =='curvature':
+                # ------------- Run Matlab symbolic minimization program 'symMinimization'
+                eng = matlab.engine.start_matlab()
+                # s = eng.genpath(path + '/Matlab-Programs')
+                s = eng.genpath(path)
+                eng.addpath(s, nargout=0)
+                # print('current Matlab folder:', eng.pwd(nargout=1))
+                eng.cd('Matlab-Programs', nargout=0)  #switch to Matlab-Programs folder
+                # print('current Matlab folder:', eng.pwd(nargout=1))
+                Inp = False
+                Inp_T = True
+                print('Run symbolic Minimization...')
+                #Arguments: symMinization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath)
+                G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4)  #Name of program:symMinimization
+                # G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4)  #Optional: add Path
+                G = np.asarray(G) #cast Matlab Outout to numpy array
+                # --- print Output ---
+                print('Minimizer G:')
+                print(G)
+                print('angle:', angle)
+                print('type:', type )
+                print('curvature:', kappa)
+                if yName =='angle':
+                    print('angle used')
+                    Y_Values.append(angle)
+                if yName =='type':
+                    print('angle used')
+                    Y_Values.append(angle)
+                if yName =='kappa':
+                    print('angle used')
+                    Y_Values.append(angle)
+    # ------------end of for-loop -----------------
+    print("(Output) Values of " + yName + ": ", Y_Values)
+# ----------------end of if-statement -------------
+
+
+
+
+
+
+
+
+# ---------------- Create Plot -------------------
+plt.figure()
+# plt.title(r''+ yName + '-Plot')
+plt.plot(X_Values, Y_Values)
+plt.scatter(X_Values, Y_Values)
+# plt.axis([0, 6, 0, 20])
+plt.xlabel(xName)
+plt.ylabel(yName)
+
+
+# plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+# plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+# plt.legend()
+plt.show()
+# #---------------------------------------------------------------