Cell-Problem-New.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/utility/structuredgridfactory.hh> //TEST
#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/functions/functionspacebases/hierarchicvectorwrapper.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/microstructure/CorrectorComputer.hh>    
#include <dune/microstructure/EffectiveQuantitiesComputer.hh>  
#include <dune/microstructure/prestrainedMaterial.hh>  

#include <dune/solvers/solvers/umfpacksolver.hh>  //TEST 
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number 

// #include <dune/fufem/discretizationerror.hh>
#include <dune/fufem/dunepython.hh>
#include <python2.7/Python.h>

// #include <boost/multiprecision/cpp_dec_float.hpp>
#include <any>
#include <variant>
#include <string>
#include <iomanip>   // needed when working with relative paths e.g. from python-scripts

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();
}

//////////////////////////////////////////////////
//   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]))
                        );
                };



// // a function:
// int half(int x, int y) {return x/2+y/2;}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char *argv[])
{
  MPIHelper::instance(argc, argv);
  
  Dune::Timer globalTimer;

  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");

  //--- setup Log-File
  std::fstream log;
  log.open(outputPath + "/output.txt" ,std::ios::out);

  std::cout << "outputPath:" << outputPath << std::endl;
  
//   parameterSet.report(log); // short Alternativ
  
    //--- Get Path for Material/Geometry functions
    std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath");
    //--- Start Python interpreter
    Python::start();
    Python::Reference main = Python::import("__main__");
    Python::run("import math");
    Python::runStream()
        << std::endl << "import sys"
        << std::endl << "sys.path.append('" << geometryFunctionPath << "')"
        << std::endl;


  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
  ///////////////////////////////////
  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",1.0);;
  double mu2 = beta*mu1;
  double lambda1 = parameterSet.get<double>("lambda1",0.0);;
  double lambda2 = beta*lambda1;
  

  if(imp == "material_neukamm")
  {
      std::cout <<"mu: " << parameterSet.get<std::array<double,3>>("mu", {1.0,3.0,2.0}) << std::endl;
      std::cout <<"lambda: " << parameterSet.get<std::array<double,3>>("lambda", {1.0,3.0,2.0}) << std::endl;
  }
  else
  {
      std::cout <<"mu: "     <<  parameterSet.get<double>("mu1",1.0) << std::endl;
      std::cout <<"lambda: " <<  parameterSet.get<double>("lambda1",0.0) << std::endl;
  }
  
  ///////////////////////////////////
  // Get Prestrain/Parameters
  ///////////////////////////////////
  auto prestrainImp = PrestrainImp<dim>(); //NEW 
  auto B_Term = prestrainImp.getPrestrain(parameterSet);

  log << "----- Input Parameters -----: " << 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 (-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});

  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:: | 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}
  for(size_t level = numLevels[0] ; level <= numLevels[1]; level++)     // levels von bis.. [2,4]
  {
    std::cout << " ----------------------------------" << std::endl;
    std::cout << "GridLevel: " << 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 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();
    if(print_debug)
       std::cout << "Host grid has " << gridView_CE.size(dim) << " vertices." << std::endl;
    
    // //not needed
    using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
    using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
    using Func2Tensor = std::function< MatrixRT(const Domain&) >;
    // 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 muGridF  = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
    auto muLocal = localFunction(muGridF);
    auto lambdaGridF  = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
    auto lambdaLocal = localFunction(lambdaGridF);


    //--- Choose a finite element space for Cell Problem
    using namespace Functions::BasisFactory;
    Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;

    //--- get PeriodicIndices for periodicBasis (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)});
            }

    //--- setup first order periodic Lagrange-Basis
    auto Basis_CE = makeBasis(
        gridView_CE,
        power<dim>(                                                                             // eig dimworld?!?! 
        Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
        flatLexicographic()
        //blockedInterleaved()   // Not Implemented
        ));     
    if(print_debug)
       std::cout << "power<periodic> basis has " << Basis_CE.dimension() << " degrees of freedom" << std::endl;
    




    //TEST 
    //Read from Parset...
    // int Phases = parameterSet.get<int>("Phases", 3);


    std::string materialFunctionName = parameterSet.get<std::string>("materialFunction", "material");
    Python::Module module = Python::import(materialFunctionName);
    // auto indicatorFunction = Python::make_function<double>(module.get("f"));
    // Func2Tensor indicatorFunction = Python::make_function<double>(module.get("f"));
    // auto materialFunction_ = Python::make_function<double>(module.get("f"));
    // auto materialFunction_ = Python::make_function<double>(module.get("f"));
    auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));

    int Phases;
    module.get("Phases").toC<int>(Phases);
    std::cout << "Number of Phases used:" << Phases << std::endl;


    // std::cout << typeid(mu_).name() << '\n';


    //---- Get mu/lambda values (for isotropic material) from Material-file
    FieldVector<double,3> mu_(0);
    module.get("mu_").toC<FieldVector<double,3>>(mu_);
    printvector(std::cout, mu_ , "mu_", "--");
    FieldVector<double,3> lambda_(0);
    module.get("lambda_").toC<FieldVector<double,3>>(lambda_);
    printvector(std::cout, lambda_ , "lambda_", "--");


    //////////////////////////////////////////////////////////////////////////////////////////////////////////
    // TESTING PRESTRAINEDMATERIAL.HH:
    using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
    
    auto material_ = prestrainedMaterial(gridView_CE,parameterSet);
    // Func2Tensor elasticityTensor = material_.getElasticityTensor();
    // auto elasticityTensor = material_.getElasticityTensor();
    // Func2TensorParam elasticityTensor_ = *material_.getElasticityTensor();
    auto elasticityTensor_ = material_.getElasticityTensor();

    Func2TensorParam TestTensor = Python::make_function<MatrixRT>(module.get("H"));

    // std::cout << "decltype(elasticityTensor_) " << decltype(elasticityTensor_) << std::endl;


    std::cout <<"typeid(elasticityTensor).name() :" <<  typeid(elasticityTensor_).name() << '\n';
    std::cout << "typeid(TestTensor).name() :" << typeid(TestTensor).name() << '\n';

    using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
    // std::cout << "Import NOW:" << std::endl;
    // MatrixFunc symTest = Python::make_function<MatrixRT>(module.get("sym"));

    // using MatrixDomainFunc = std::function< MatrixRT(const MatrixRT&,const Domain&)>;
   // // MatrixFunc elasticityTensor = Python::make_function<MatrixRT>(module.get("L"));

    // auto elasticityTensorGVF  = Dune::Functions::makeGridViewFunction(elasticityTensor , Basis_CE.gridView());
    // auto localElasticityTensor = localFunction(elasticityTensorGVF);

    // Func2Tensor forceTerm = [] (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?
    //                     };

    // auto loadGVF  = Dune::Functions::makeGridViewFunction(forceTerm, Basis_CE.gridView());
    // auto loadFunctional = localFunction(loadGVF);  

    MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
    // auto xu = symTest(G1_);
    // std::cout << "TEST NOW:" << std::endl;
    // printmatrix(std::cout, symTest(G1_), "symTest(G1_)", "--");

    // auto TestTensorGVF  = Dune::Functions::makeGridViewFunction(TestTensor , Basis_CE.gridView());
    // auto localTestTensor = localFunction(TestTensorGVF );


    // printmatrix(std::cout, elasticityTensor(G1_), "elasticityTensor(G1_)", "--");
    // auto temp = elasticityTensor(G1_);

    

    for (const auto& element : elements(Basis_CE.gridView()))
    {
  
        int orderQR = 2;
        const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
        for (const auto& quadPoint : quad)
        {
            const auto& quadPos = quadPoint.position();
            // std::cout << "quadPos : " << quadPos  << std::endl;
            auto temp = TestTensor(G1_, element.geometry().global(quadPos));
            auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
            // std::cout << "material_.applyElasticityTensor:" << std::endl;
            auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos));
            // printmatrix(std::cout, tmp3, "tmp3", "--");
        }
    }



    // for (auto&& vertex : vertices(gridView_CE))
    // {
    //     std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
    //     auto tmp = vertex.geometry().corner(0);
    //     auto temp = elasticityTensor(tmp);
    //     // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
    //     // printmatrix(std::cout, localElasticityTensor(G1_,tmp), "localElasticityTensor(vertex.geometry().corner(0))", "--");
    // }

    

    // std::function<int(int,int)> fn1 = half; 
    // std::cout << "fn1(60,20): " << fn1(60,20) << '\n';


    // std::cout << typeid(elasticityTensorGVF).name() << '\n';
    // std::cout << typeid(localElasticityTensor).name() << '\n';


    // ParameterTree parameterSet_2;
    // ParameterTreeParser::readINITree(geometryFunctionPath + "/"+ materialFunctionName + ".py", parameterSet_2);

    // auto lu = parameterSet_2.get<FieldVector<double,3>>("lu", {1.0,3.0,2.0});
    // std::cout <<"lu[1]: " << lu[1]<< std::endl;
    // std::cout <<"lu: " << parameterSet_2.get<std::array<double,3>>("lu", {1.0,3.0,2.0}) << std::endl;
    
    // auto mU_ = module.evaluate(parameterSet_2.get<std::string>("lu", "[1,2,3]"));
    // std::cout << "typeid(mU_).name()" << typeid(mU_.operator()()).name() << '\n';


    // for (auto&& vertex : vertices(gridView_CE))
    // {
    //     std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
    //     // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
    //     printvector(std::cout, materialFunction_(vertex.geometry().corner(0)), "materialFunction_(vertex.geometry().corner(0))", "--");
    // }
    // std::cout << "materialFunction_({0.0,0.0,0.0})", materialFunction_({0.0,0.0,0.0}) << std::endl;
    


    // --------------------------------------------------------------

    //TODO// Phasen anhand von Mu bestimmen?
    //TODO: DUNE_THROW(Exception, "Inconsistent choice of interpolation method");   if number of Phases != mu/lambda parameters

    //FÜR L GARNICHT NÖTIG DENN RÜCKGABETYPE IS IMMER MATRIXRT!?!:
    // BEi materialfunction (isotopic) reicht auch FieldVector<double,2> für lambda/mu

    // switch (Phases)
    // {
    //     case 1: //homogeneous material
    //     {
    //       std::cout << "Phase - 1" << std::endl;
    //       auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
    //       break;
    //     }
    //     case 2: 
    //     {
    //       std::cout << "Phase - 1" << std::endl;
    //       auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
    //       break;
    //     }
    //     case 3:
    //     {
    //       std::cout << "Phase - 3" << std::endl;
    //       auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
    //       break;
    //     }
    // }
    

    // switch (Phases)
    // {
    //     case 1: //homogeneous material
    //     {
    //       std::cout << "Phases - 1" << std::endl;
    //       std::array<double,1> mu_ = parameterSet.get<std::array<double,1>>("mu", {1.0});
    //       Python::Module module = Python::import(materialFunction);
    //       auto indicatorFunction = Python::make_function<double>(module.get("f"));  // get indicator function
    //       auto muTerm = [mu_] (const Domain& x) {return mu_;};
    //       break;
    //     }
    //     case 2: 
    //     {
    //       std::cout << "Phases - 2" << std::endl;
    //       std::array<double,2> mu_ = parameterSet.get<std::array<double,2>>("mu", {1.0,3.0});
    //       Python::Module module = Python::import(materialFunction);
    //       auto indicatorFunction = Python::make_function<double>(module.get("f"));  // get indicator function
    //       auto muTerm = [mu_,indicatorFunction] (const Domain& x) 
	// 	      {
    //         if (indicatorFunction(x) == 1) 
    //           return mu_[0];
    //         else 
    //           return mu_[1];
    //         };
    //       break;
    //     }
    //     case 3:
    //     {
    //       std::cout << "Phases - 3" << std::endl;
    //       std::array<double,3> mu_ = parameterSet.get<std::array<double,3>>("mu", {1.0,3.0, 5.0});
    //       Python::Module module = Python::import(materialFunction);
    //       auto indicatorFunction = Python::make_function<double>(module.get("f"));  // get indicator function
    //       auto muTerm = [mu_,indicatorFunction] (const Domain& x) 
    //             {
    //               if (indicatorFunction(x) == 1) 
    //                 return mu_[0];
    //               else if (indicatorFunction(x) == 2) 
    //                 return mu_[1];
    //               else 
    //                 return mu_[2];
    //             };
    //       break;
    //     }
    // }

    

    
    //TEST 
//     std::cout << "Test crossSectionDirectionScaling:" << std::endl;
/*    
    MatrixRT T {{1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}};
    printmatrix(std::cout, T, "Matrix T", "--");
    
    auto ST = crossSectionDirectionScaling((1.0/5.0),T);
    printmatrix(std::cout, ST, "scaled Matrix T", "--");*/
    
    //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);
//                         };


    //------------------------------------------------------------------------------------------------
    //--- compute Correctors
    // auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
    auto correctorComputer = CorrectorComputer(Basis_CE, material_, muTerm, lambdaTerm, gamma, log, parameterSet);
    correctorComputer.solve();



//////////////////////////////////////////////////
   

    //--- 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:
    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,material_);
    effectiveQuantitiesComputer.computeEffectiveQuantities();


      }
 /*

    //--- 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 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<< std::fixed << std::setprecision(6) << "q_onetwo=" << Qeff[0][1] << std::endl;
    // -------------------------------------------


    //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?
    //                     };

    // 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 << "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;
    log << "b3=" << Beff[2] << std::endl;
    log << "mu_gamma=" << Qeff[2][2] << std::endl;           // added for Python-Script




   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});
//         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},{40,40,40});
        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},nElements);
        
        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});
// //     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},{40,40,40});
//     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},nElements);
//     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");;
    
    
//     // TEST 
//     auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
//     auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());

//     std::vector<double> B_CoeffP0;
//     Functions::interpolate(scalarP0FeBasis, B_CoeffP0, B_Term);
//     auto B_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, B_CoeffP0);
        

        
//     VTKWriter<GridView> PrestrainVtkWriter(gridView_VTK);
         
//     PrestrainVtkWriter.addCellData(
//             B_DGBF_P0,
//             VTK::FieldInfo("B_P0", VTK::FieldInfo::Type::scalar, 1));     
        
//     PrestrainVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
//     std::cout << "wrote data to file:" + outputPath +"/PrestrainFunctions-level" + std::to_string(level) << std::endl; 

//    }
  
  


  levelCounter++; 
  } // Level-Loop End
  



    //////////////////////////////////////////
    //--- Print Storage
    int tableWidth = 12;
    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) 
    {
        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 % 9 == 0 )
        {
            std::cout << std::endl;
            log << std::endl;
        }
    }
    std::cout << std::string(tableWidth*9 + 3*9, '-') << "\n";
    log       << std::string(tableWidth*9 + 3*9, '-') << "\n";  

    log.close();

    std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;

    */


}