diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 1ff3364e6a797c7c7ba6e3726455ba7f29d9edb2..9fdf01aa23bd0ef4e34fedd46ddd7bc095bb496d 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -126,40 +126,6 @@ public:
   {
       assembleCellStiffness(stiffnessMatrix_);
 
-
-      // // lateral stretch strain		E_U1 = e1 ot e1  				= MatrixRT{{1, 0, 0}, {0, 0, 0}, {0, 0, 0}};
-      // 	Func2Tensor neg_E_a  = [] (const Domain& z) { return biotStrainApprox({-1, 0, 0}, {0, 0, 0}, {0, z[dim-2], z[dim-1]}); }; 
-
-      // // twist in x2-x3 plane 		E_K1 = (e1 x (0,x2,x3)^T ot e1) = MatrixRT{{0,-z[2], z[1]}, {0, 0 , 0 }, {0,0,0}}; 
-      // Func2Tensor neg_E_K1 = [] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {-1, 0, 0}, {0, z[dim-2], z[dim-1]}); };
-
-      // //  bend strain in x1-x2 plane	E_K2 = (e2 x (0,x2,x3)^T ot e1) = MatrixRT{{z[2], 0, 0}, {0, 0, 0}, {0, 0, 0}}; 
-      // Func2Tensor neg_E_K2 = [] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {0, -1, 0}, {0, z[dim-2], z[dim-1]}); };
-      
-      // // bend strain in x1-x3 plane	E_K3 = (e3 x (0,x2,x3)^T ot e1) = MatrixRT{{-z[1], 0, 0}, {0, 0, 0}, {0, 0, 0}};
-      // Func2Tensor neg_E_K3 = [] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {0, 0, -1}, {0, z[dim-2], z[dim-1]});  };
-
-      // assembleLoadVector(neg_E_a, load_a_); 
-      // //log_ << "\n\n\n\n\n\n";
-      // assembleLoadVector(neg_E_K1, load_K1_);
-      // assembleLoadVector(neg_E_K2, load_K2_);
-      // assembleLoadVector(neg_E_K3, load_K3_);
-  /////////////////////////////////////////////////////////
-  // Define "rhs"-functions
-  /////////////////////////////////////////////////////////
-
-                      
-
-  // 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);};
-  // Func2Tensor x3G_1neg = [] (const Domain& x) {return -1.0*x3G_1_(x);};
-  // Func2Tensor x3G_2neg = [] (const Domain& x) {return -1.0*x3G_2_(x);};
-  // Func2Tensor x3G_3neg = [] (const Domain& x) {return -1.0*x3G_3_(x);};
-      // assembleCellLoad(load_alpha1_ ,x3G_1neg);
-      // assembleCellLoad(load_alpha2_ ,x3G_2neg);
-      // assembleCellLoad(load_alpha3_ ,x3G_3neg);
-
       assembleCellLoad(load_alpha1_ ,x3G_1_);
       assembleCellLoad(load_alpha2_ ,x3G_2_);
       assembleCellLoad(load_alpha3_ ,x3G_3_);
@@ -172,13 +138,10 @@ public:
     ///////////////////////////////
     const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}
 
-    // std::map<int,int> getIndexMapPeriodicBoundary(){return perIndexMap_;}
-
-    ParameterTree getParameterSet() const {return parameterSet_;} 
+    ParameterTree getParameterSet() const {return parameterSet_;}
 
     fstream* getLog(){return &log_;}
 
-
     double getGamma(){return gamma_;}
 
     shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
@@ -202,6 +165,7 @@ public:
     // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
     auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
 
+
   
 
     
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index b6fde15d5b3b61d6647fb4d09b98395d0c22dcf5..d4f4b2548c1592ee6c93a4588394eea518b1d0aa 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -10,6 +10,7 @@
 #include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number 
 #include <dune/istl/io.hh>
 #include <dune/istl/matrix.hh>
+#include <dune/common/parametertree.hh>
 
 using namespace Dune;
 using namespace MatrixOperations;
@@ -143,6 +144,7 @@ public:
         auto lambda_ = *correctorComputer_.getLambda();
         auto gamma = correctorComputer_.getGamma();
         auto basis = *correctorComputer_.getBasis();
+        ParameterTree parameterSet = correctorComputer_.getParameterSet();
 
 		shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(), 
                                             correctorComputer_.getCorr_phi2(),
@@ -238,9 +240,11 @@ public:
             if (b==0)
                 Bhat_[a] = prestrain;
         }
-        printmatrix(std::cout, Q_, "Matrix Q_", "--");
-        printvector(std::cout, Bhat_, "Bhat_", "--");
-
+        if(parameterSet.get<bool>("print_debug", false))
+        {
+            printmatrix(std::cout, Q_, "Matrix Q_", "--");
+            printvector(std::cout, Bhat_, "Bhat_", "--");
+        }
 
         ///////////////////////////////
         // Compute effective Prestrain B_eff (by solving linear system)
@@ -250,7 +254,8 @@ public:
         // MatrixInfo<MatrixRT> matrixInfo(Q_,true,2,1);
         // std::cout << "----------------------------------------" << std::endl;
         Q_.solve(Beff_,Bhat_);
-        printvector(std::cout, Beff_, "Beff_", "--");
+        if(parameterSet.get<bool>("print_debug", false))
+            printvector(std::cout, Beff_, "Beff_", "--");
         
     
         //LOG-Output
@@ -260,7 +265,6 @@ public:
         log << "Beff_: " << Beff_ <<  " (Effective Prestrain)" << std::endl;
         log << "------------------------ " << std::endl;
 
-
         //   TEST
         //   std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
         return ;
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 5f4db4625b202c12314d0fb10a988784b4c24b47..9f0192d8d07c83221b7cba98b10ee07a1da8d569 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -13,16 +13,10 @@
 outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
 
 
-### DEBUG Option:
-print_debug = true
+# --- DEBUG Option:
+print_debug = true  #(default=false)
 
 
-#############################################
-#  Cell Domain
-#############################################
-# Domain 1: (-1/2, 1/2)^3  , Domain 2 : [0,1)^2 x (-1/2, 1/2)
-cellDomain=1
-
 #############################################
 #  Grid parameters
 #############################################
@@ -106,45 +100,47 @@ material_prestrain_imp= "material_neukamm"   #(Python-indicator-function with sa
 #rF = 0.1 # Fiber-Radius   (default = 0.5*(width/(2.0*nF))  //half of the max-fiber-radius mrF = (width/(2.0*nF)) )
 
 
-# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files (default=false):
-write_materialFunctions = true
-write_prestrainFunctions = true  # VTK norm of B ,
-write_VTK = true
-
-
-# --- (Optional output) L2Error, compute integral mean:
-#write_L2Error = true
-#write_IntegralMean = true
-
 
 #############################################
 #  Assembly options
 #############################################
-#set_IntegralZero = true
-set_IntegralZero = false
-#set_oneBasisFunction_Zero = true
-
-#arbitraryLocalIndex = 7
-#arbitraryElementNumber = 3
-#arbitraryLocalIndex = 0
-#arbitraryElementNumber = 0
+#set_IntegralZero = true            #(default = false)
+#set_oneBasisFunction_Zero = true   #(default = false)
+
+#arbitraryLocalIndex = 7            #(default = 0)
+#arbitraryElementNumber = 3         #(default = 0)
 #############################################
 
 
 #############################################
-#  Solver Type: #1: CG - SOLVER (default), #2: GMRES - SOLVER, #3: QR - SOLVER  #4: UMFPACK - SOLVER
+#  Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
 #############################################
 Solvertype = 3
 Solver_verbosity = 0  #(default = 2)  degree of information for solver output
 
 
+#############################################
+#  Write/Output options                 #(default=false)
+#############################################
 
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+write_materialFunctions = true
+write_prestrainFunctions = true  # VTK norm of B ,
 
+# --- Write Correctos to VTK-File:  
+write_VTK = true
+
+# --- (Optional output) L2Error, integral mean: 
+#write_L2Error = true                   
+#write_IntegralMean = true              
 
-# --- Write corrector-coefficients to log-File (default=false):
-#write_corrector_phi1 = false
-#write_corrector_phi2 = false
-#write_corrector_phi3 = false
+# --- check orthogonality (75) from paper: 
+write_checkOrthogonality = true         
+
+# --- Write corrector-coefficients to log-File:
 #write_corrector_phi1 = true
 #write_corrector_phi2 = true
 #write_corrector_phi3 = true
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+write_toMATLAB = false  
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 2a743c255ae133a5555a5e8eba22ed1eb4cb6035..b7ae0127d8f60344f93fbd3223c1d3c7a5078487 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -12,13 +12,11 @@
 #include <dune/common/math.hh>
 
 
-
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/grid/uggrid.hh>
 #include <dune/grid/yaspgrid.hh>
-#include <dune/grid/utility/structuredgridfactory.hh> //TEST
-
+// #include <dune/grid/utility/structuredgridfactory.hh> //TEST
 #include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
 
 #include <dune/istl/matrix.hh>
@@ -32,18 +30,16 @@
 #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/functions/functionspacebases/hierarchicvectorwrapper.hh>
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
@@ -55,18 +51,13 @@
 #include <dune/microstructure/EffectiveQuantitiesComputer.hh>  
 
 #include <dune/solvers/solvers/umfpacksolver.hh>  //TEST 
-
 #include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number 
 
-#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
-
+// #include <dune/fufem/discretizationerror.hh>
 #include <dune/fufem/dunepython.hh>
 #include <python2.7/Python.h>
 
-// #include <dune/fufem/discretizationerror.hh>
 // #include <boost/multiprecision/cpp_dec_float.hpp>
-// #include <boost/any.hpp>
-
 #include <any>
 #include <variant>
 #include <string>
@@ -75,7 +66,6 @@
 using namespace Dune;
 using namespace MatrixOperations;
 
-
 //////////////////////////////////////////////////////////////////////
 // Helper functions for Table-Output
 //////////////////////////////////////////////////////////////////////
@@ -118,8 +108,6 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
                         );
                 };
 
-
-
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 int main(int argc, char *argv[])
@@ -162,6 +150,17 @@ int main(int argc, char *argv[])
 
   constexpr int dim = 3;
   constexpr int dimWorld = 3;
+
+  // Debug/Print Options 
+  bool print_debug = parameterSet.get<bool>("print_debug", false);
+
+  // VTK-write options
+  bool write_materialFunctions   = parameterSet.get<bool>("write_materialFunctions", false);
+  bool write_prestrainFunctions  = parameterSet.get<bool>("write_prestrainFunctions", false);
+
+
+
+
   ///////////////////////////////////
   // Get Parameters/Data
   ///////////////////////////////////
@@ -223,7 +222,7 @@ int main(int argc, char *argv[])
   ///////////////////////////////////
   //--- 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|           
+  //--- Storage:: | level | q1 | q2 | q3 | q12 | q23 | b1 | b2 | b3 |           
   std::vector<std::variant<std::string, size_t , double>> Storage_Quantities;         
 
   //   for(const size_t &level : numLevels)                             // explixite Angabe der levels.. {2,4}
@@ -236,21 +235,22 @@ int main(int argc, char *argv[])
     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;
+    std::cout << "Number of Grid-Elements in each direction: " << nElements << std::endl;
+    log << "Number of Grid-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();
-    std::cout << "Host grid has " << gridView_CE.size(dim) << " vertices." << std::endl;
+    if(print_debug)
+       std::cout << "Host grid has " << gridView_CE.size(dim) << " vertices." << std::endl;
     
-    //TODO needed?
-    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> >;
+    // //not needed
+    // 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
@@ -265,14 +265,6 @@ int main(int argc, char *argv[])
     auto lambdaLocal = localFunction(lambdaGridF);
 
 
-    // Print Options 
-    bool print_debug = parameterSet.get<bool>("print_debug", false);
-
-    //VTK-Write 
-    bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
-    bool write_prestrainFunctions  = parameterSet.get<bool>("write_prestrainFunctions", false);
-
-
     //--- Choose a finite element space for Cell Problem
     using namespace Functions::BasisFactory;
     Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
@@ -293,25 +285,25 @@ int main(int argc, char *argv[])
         flatLexicographic()
         //blockedInterleaved()   // Not Implemented
         ));     
-    std::cout << "power<periodic> basis has " << Basis_CE.dimension() << " degrees of freedom" << std::endl;
+    if(print_debug)
+       std::cout << "power<periodic> basis has " << Basis_CE.dimension() << " degrees of freedom" << std::endl;
     
 
 
 
 
 
-
-
-
-
     //Read from Parset...
     int Phases = parameterSet.get<int>("Phases", 3);
     std::string materialFunction = parameterSet.get<std::string>("materialFunction", "material");
 
-      Python::Module module = Python::import("material");
-      auto indicatorFunction = Python::make_function<double>(module.get("f"));
+    Python::Module module = Python::import("material");
+    auto indicatorFunction = Python::make_function<double>(module.get("f"));
 
-    std::cout << "Phases:" << Phases << std::endl;
+    std::cout << "Phases:" << Phases << std::endl;  
+    
+    //TODO// Phasen anhand von Mu bestimmen?
+    //TODO: DUNE_THROW(Exception, "Inconsistent choice of interpolation method");   if number of Phases != mu/lambda parameters
 
     
     // switch (Phases)
@@ -377,43 +369,50 @@ int main(int argc, char *argv[])
 
 
 
-
+    //------------------------------------------------------------------------------------------------
     //--- compute Correctors
     auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
-    
     correctorComputer.solve();
-    correctorComputer.computeNorms();
-    correctorComputer.writeCorrectorsVTK(level);
+
+    //--- check Correctors (options):
+    if(parameterSet.get<bool>("write_L2Error", false))
+         correctorComputer.computeNorms();
+    if(parameterSet.get<bool>("write_VTK", false))
+         correctorComputer.writeCorrectorsVTK(level);
     //--- additional Test: check orthogonality (75) from paper:
-    correctorComputer.check_Orthogonality();
-    correctorComputer.checkSymmetry();
+    if(parameterSet.get<bool>("write_checkOrthogonality", false))
+        correctorComputer.check_Orthogonality();
+    //--- Check symmetry of stiffness matrix
+    if(print_debug)
+        correctorComputer.checkSymmetry();
 
     //--- compute effective quantities
     auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term);
     effectiveQuantitiesComputer.computeEffectiveQuantities();
 
-    //TEST
-    // std::cout << "----------computeFullQ-----------"<< std::endl;
+    //--- Test:: Compute Qeff without using the orthogonality (75)... 
+    // only really makes a difference whenever the orthogonality is not satisfied!
+    // std::cout << "----------computeFullQ-----------"<< std::endl;  //TEST
     // effectiveQuantitiesComputer.computeFullQ();
 
     //--- get effective quantities
     auto Qeff = effectiveQuantitiesComputer.getQeff();
     auto Beff = effectiveQuantitiesComputer.getBeff();
-  //  printmatrix(std::cout, Qeff, "Matrix Q_T", "--");
-  //  printvector(std::cout, Beff, "Beff", "--");
+    printmatrix(std::cout, Qeff, "Matrix Qeff", "--");
+    printvector(std::cout, Beff, "Beff", "--");
+
+    //--- write effective quantities to matlab folder (for symbolic minimization)
+    if(parameterSet.get<bool>("write_toMATLAB", false))
+        effectiveQuantitiesComputer.writeToMatlab(outputPath);
 
     std::cout.precision(10);
     std::cout<< "q1 : " << Qeff[0][0] << std::endl;
     std::cout<< "q2 : " << Qeff[1][1] << std::endl;
-    std::cout << "q3 : " << std::fixed << Qeff[2][2] << std::endl;
+    std::cout<< "q3 : " << std::fixed << Qeff[2][2] << std::endl;
     std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Qeff[0][1] << std::endl;
     // -------------------------------------------
 
 
-    //--- write effective quantities to matlab folder (for symbolic minimization)
-    effectiveQuantitiesComputer.writeToMatlab(outputPath);
-
-
     //TEST 
     // 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?
@@ -422,19 +421,21 @@ int main(int argc, char *argv[])
     // double energy = effectiveQuantitiesComputer.energySP(x3G_1,x3G_1);
     // std::cout << "energy:" << energy << std::endl;
 
-
     Storage_Quantities.push_back(Qeff[0][0] );
     Storage_Quantities.push_back(Qeff[1][1] );
     Storage_Quantities.push_back(Qeff[2][2] );
+    Storage_Quantities.push_back(Qeff[0][1] );
+    Storage_Quantities.push_back(Qeff[1][2] );
     Storage_Quantities.push_back(Beff[0]);
     Storage_Quantities.push_back(Beff[1]);
     Storage_Quantities.push_back(Beff[2]);
 
     log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
-    log << "q1=" << Qeff[0][0] << std::endl;
-    log << "q2=" << Qeff[1][1] << std::endl;
-    log << "q3=" << Qeff[2][2] << std::endl;
+    log << "q1="  << Qeff[0][0] << std::endl;
+    log << "q2="  << Qeff[1][1] << std::endl;
+    log << "q3="  << Qeff[2][2] << std::endl;
     log << "q12=" << Qeff[0][1] << std::endl;
+    log << "q23=" << Qeff[1][2] << std::endl;
     log << std::fixed << std::setprecision(6) << "q_onetwo=" << Qeff[0][1] << std::endl;
     log << "b1=" << Beff[0] << std::endl;
     log << "b2=" << Beff[1] << std::endl;
@@ -544,23 +545,27 @@ int main(int argc, char *argv[])
     //////////////////////////////////////////
     //--- Print Storage
     int tableWidth = 12;
-    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";   
+    std::cout << center("Levels ",tableWidth)   << " | "
+              << center("q1",tableWidth)        << " | "
+              << center("q2",tableWidth)        << " | "
+              << center("q3",tableWidth)        << " | "
+              << center("q12",tableWidth)       << " | "
+              << center("q23",tableWidth)       << " | "
+              << center("b1",tableWidth)        << " | "
+              << center("b2",tableWidth)        << " | "
+              << center("b3",tableWidth)        << " | " << "\n";
+    std::cout << std::string(tableWidth*9 + 3*9, '-')    << "\n";
+    log       << std::string(tableWidth*9 + 3*9, '-')    << "\n";   
+    log       << center("Levels ",tableWidth)   << " | "
+              << center("q1",tableWidth)        << " | "
+              << center("q2",tableWidth)        << " | "
+              << center("q3",tableWidth)        << " | "
+              << center("q12",tableWidth)       << " | "
+              << center("q23",tableWidth)       << " | "
+              << center("b1",tableWidth)        << " | "
+              << center("b2",tableWidth)        << " | "
+              << center("b3",tableWidth)        << " | " << "\n";
+    log       << std::string(tableWidth*9 + 3*9, '-')    << "\n";   
   
     int StorageCount2 = 0;
     for(auto& v: Storage_Quantities) 
@@ -568,14 +573,14 @@ int main(int argc, char *argv[])
         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 )
+        if(StorageCount2 % 9 == 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";  
+    std::cout << std::string(tableWidth*9 + 3*9, '-') << "\n";
+    log       << std::string(tableWidth*9 + 3*9, '-') << "\n";  
 
     log.close();
 
diff --git a/src/deprecated_code/cellsolver.parset b/src/deprecated_code/cellsolver.parset
new file mode 100644
index 0000000000000000000000000000000000000000..5f4db4625b202c12314d0fb10a988784b4c24b47
--- /dev/null
+++ b/src/deprecated_code/cellsolver.parset
@@ -0,0 +1,150 @@
+# --- Parameter File as Input for 'Cell-Problem'
+#
+# NOTE: define variables without whitespaces in between! i.e. : gamma=1.0 instead of gamma = 1.0
+# since otherwise these cant be read from other Files!
+# --------------------------------------------------------
+
+
+#############################################
+#  Choose Output-path for logfile
+#############################################
+### Remove/Comment this when running via Python-Script:
+#outputPath=/home/stefan/DUNE/dune-microstructure/outputs
+outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
+
+
+### DEBUG Option:
+print_debug = true
+
+
+#############################################
+#  Cell Domain
+#############################################
+# Domain 1: (-1/2, 1/2)^3  , Domain 2 : [0,1)^2 x (-1/2, 1/2)
+cellDomain=1
+
+#############################################
+#  Grid parameters
+#############################################
+#----------------------------------------------------
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+
+numLevels= 2 2
+#numLevels =  0 0   # computes all levels from first to second entry
+#numLevels =  2 2   # computes all levels from first to second entry
+#numLevels =  1 3   # computes all levels from first to second entry
+#numLevels = 4 4    # computes all levels from first to second entry
+#numLevels = 5 5    # computes all levels from first to second entry
+#numLevels = 6 6    # computes all levels from first to second entry
+#numLevels = 1 6
+
+
+
+########################################################################################
+
+# --- Choose scale ratio gamma:
+gamma=1.0
+
+
+#############################################
+#  Material / Prestrain parameters and ratios
+#############################################
+
+#-- stiffness-ratio (ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case)
+beta=3.0
+
+$--- strength of prestrain
+rho1=1.0
+
+#--- prestrain-contrast (ratio between prestrain parameters rho1 & rho2)
+alpha=8.0
+
+
+#--- Lame-Parameters
+mu1=80.0
+lambda1=80.0
+
+
+# better: pass material parameters as a vector
+# Poisson ratio nu = 1/4 => lambda = 0.4*mu
+#mu=1.2 1.2 1
+#lambda=0.48 0.48 0.4
+#rho=1.0 0
+mu=80 80 60
+lambda=80 80 25
+
+#mu=1 
+#lambda=4 
+
+
+
+# ---volume fraction  (default value = 1.0/4.0)
+theta=0.25
+#theta = 0.3
+#theta = 0.75
+#theta = 0.125
+#theta = 0.5
+
+
+
+#--- choose composite-Implementation:
+#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
+geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
+material_prestrain_imp= "material_neukamm"   #(Python-indicator-function with same name determines material phases)
+#material_prestrain_imp= "two_phase_material_2"  #(Python-indicator-function with same name determines material phases)
+#material_prestrain_imp= "two_phase_material_3"  #(Python-indicator-function with same name determines material phases)
+#material_prestrain_imp= "parametrized_Laminate"
+#material_prestrain_imp= "homogeneous"
+#material_prestrain_imp= "analytical_Example"
+#material_prestrain_imp="isotropic_bilayer"
+#material_prestrain_imp= "circle_fiber"    #TEST
+#material_prestrain_imp= "matrix_material_circles"
+#material_prestrain_imp= "matrix_material_squares"
+#nF = 10 # Number of Fibers (default = 3)
+#rF = 0.1 # Fiber-Radius   (default = 0.5*(width/(2.0*nF))  //half of the max-fiber-radius mrF = (width/(2.0*nF)) )
+
+
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files (default=false):
+write_materialFunctions = true
+write_prestrainFunctions = true  # VTK norm of B ,
+write_VTK = true
+
+
+# --- (Optional output) L2Error, compute integral mean:
+#write_L2Error = true
+#write_IntegralMean = true
+
+
+#############################################
+#  Assembly options
+#############################################
+#set_IntegralZero = true
+set_IntegralZero = false
+#set_oneBasisFunction_Zero = true
+
+#arbitraryLocalIndex = 7
+#arbitraryElementNumber = 3
+#arbitraryLocalIndex = 0
+#arbitraryElementNumber = 0
+#############################################
+
+
+#############################################
+#  Solver Type: #1: CG - SOLVER (default), #2: GMRES - SOLVER, #3: QR - SOLVER  #4: UMFPACK - SOLVER
+#############################################
+Solvertype = 3
+Solver_verbosity = 0  #(default = 2)  degree of information for solver output
+
+
+
+
+
+# --- Write corrector-coefficients to log-File (default=false):
+#write_corrector_phi1 = false
+#write_corrector_phi2 = false
+#write_corrector_phi3 = false
+#write_corrector_phi1 = true
+#write_corrector_phi2 = true
+#write_corrector_phi3 = true
diff --git a/inputs/cellsolver_backup.parset b/src/deprecated_code/cellsolver_backup.parset
similarity index 100%
rename from inputs/cellsolver_backup.parset
rename to src/deprecated_code/cellsolver_backup.parset
diff --git a/inputs/computeMuGamma_backup.parset b/src/deprecated_code/computeMuGamma_backup.parset
similarity index 100%
rename from inputs/computeMuGamma_backup.parset
rename to src/deprecated_code/computeMuGamma_backup.parset