CorrectorComputer.hh

#ifndef DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH

#include <map>

#include <dune/common/parametertree.hh>
// #include <dune/common/float_cmp.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh> 
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh> 
#include <dune/microstructure/matrix_operations.hh>

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

#include <dune/microstructure/voigthelper.hh>

using namespace MatrixOperations;

using std::shared_ptr;
using std::make_shared;
using std::fstream;


template <class Basis, class Material> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
class CorrectorComputer {

public:
	static const int dimworld = 3; //GridView::dimensionworld;
	static const int dim = Basis::GridView::dimension; //const int dim = Domain::dimension;
	
	using GridView = typename Basis::GridView;
	using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
	using VectorRT = Dune::FieldVector< double, dimworld>;
	using MatrixRT = Dune::FieldMatrix< double, dimworld, dimworld>;
	using FuncScalar = std::function< double(const Domain&) >;
	using FuncVector = std::function< VectorRT(const Domain&) >;
	using Func2Tensor = std::function< MatrixRT(const Domain&) >;
  using Func2int = std::function< int(const Domain&) >;
	using VectorCT = Dune::BlockVector<Dune::FieldVector<double,1> >;
	using MatrixCT = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
	using ElementMatrixCT = Dune::Matrix<double>;

protected:
//private:
	const Basis& basis_; 


  const Material& material_;

  double gamma_;

	fstream& log_;      // Output-log
	const Dune::ParameterTree& parameterSet_;

	// const FuncScalar mu_; 
  // const FuncScalar lambda_; 

  MatrixCT stiffnessMatrix_; 
	VectorCT load_alpha1_,load_alpha2_,load_alpha3_; //right-hand side(load) vectors

  VectorCT x_1_, x_2_, x_3_;            // (all) Corrector coefficient vectors 
  VectorCT phi_1_, phi_2_, phi_3_;      // Corrector phi_i coefficient vectors 
  Dune::FieldVector<double,3> m_1_, m_2_, m_3_;  // Corrector m_i coefficient vectors

  // (assembled) corrector matrices M_i
  std::array<MatrixRT, 3 > mContainer;
  const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};

  // ---- Basis for R_sym^{2x2}
  // Could really be constexpr static, but then we would need to write the basis here directly in Voigt notation.
  // However, I suspect that our Voigt representation may still change in the future.
  const std::array<VoigtVector<double,3>,3 > matrixBasis_;

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

  const std::array<Func2Tensor, 3> x3MatrixBasisContainer_ = {x3G_1_, x3G_2_, x3G_3_};

    // --- Offset between basis indices 
  const int phiOffset_;

public: 
    ///////////////////////////////
    // constructor
    ///////////////////////////////
    CorrectorComputer( const Basis& basis, 
            const Material& material,
            double gamma,
            std::fstream& log, 
            const Dune::ParameterTree& parameterSet)
          : basis_(basis), 
            material_(material),
            gamma_(gamma),
            log_(log),
            parameterSet_(parameterSet),
            matrixBasis_(std::array<VoigtVector<double,3>,3>{matrixToVoigt(Dune::FieldMatrix<double,3,3>({{1, 0, 0}, {0, 0, 0}, {0, 0, 0}})),
                                                             matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
                                                             matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))}),
            phiOffset_(basis.size())
    {

      // assemble();

      //write VTK call here...
    } 


  // -----------------------------------------------------------------
  // --- Assemble Corrector problems
  void assemble()
  {
      Dune::Timer StiffnessTimer;
      assembleCellStiffness(stiffnessMatrix_);
      std::cout << "Stiffness assembly Timer: " << StiffnessTimer.elapsed() << std::endl;

      assembleCellLoad(load_alpha1_ ,x3G_1_);
      assembleCellLoad(load_alpha2_ ,x3G_2_);
      assembleCellLoad(load_alpha3_ ,x3G_3_);
  };



    ///////////////////////////////
    // Getter
    ///////////////////////////////
    const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}

    Dune::ParameterTree getParameterSet() const {return parameterSet_;}

    fstream* getLog(){return &log_;}

    double getGamma(){return gamma_;}

    shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
    shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
    shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
    shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}

    // shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
    // shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}

    shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}

    // --- Get Correctors
    // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
    // auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
    auto getMcontainer(){return mContainer;}
    shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}

    
    // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
    auto getMatrixBasiscontainer(){return make_shared<std::array<VoigtVector<double,3>,3 >>(matrixBasis_);}
    // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
    auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}

    
    // shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
    shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
    shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
    shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}




  // Get the occupation pattern of the stiffness matrix
  void getOccupationPattern(Dune::MatrixIndexSet& nb)
  {
    //  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_.get<bool>("set_oneBasisFunction_Zero ", true)){
      Dune::FieldVector<int,3> row;
      unsigned int arbitraryLeafIndex =  parameterSet_.get<unsigned int>("arbitraryLeafIndex", 0);
      unsigned int arbitraryElementNumber =  parameterSet_.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  < (std::size_t)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(arbitraryElementNumber,arbitraryLeafIndex);

      for (int k = 0; k<3; k++)
        nb.add(row[k],row[k]);
    }
    std::cout << "finished occupation pattern\n";
  }

  // template<class localFunction1, class localFunction2>
  // void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
  //                                   ElementMatrixCT& elementMatrix,
  //                                   const localFunction1& mu,
  //                                   const localFunction2& lambda
  //                                   )

  /** \brief Compute the stiffness matrix for one element
   *
   * \param quadRule The quadrature rule to be used
   * \param phaseAtQuadPoint The material phase (i.e., the type of material) at each quadrature point
   */
  void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
                                     const Dune::QuadratureRule<double,dim>& quadRule,
                                     const std::vector<int>& phaseAtQuadPoint,
                                    ElementMatrixCT& elementMatrix
                                    )
  {
    // Local StiffnessMatrix of the form:
    // | phi*phi    m*phi |
    // | phi *m     m*m   |
    auto element = localView.element();
    auto geometry = element.geometry();
  //   using MatrixRT = FieldMatrix< double, dimworld, dimworld>;

    elementMatrix.setSize(localView.size()+3, localView.size()+3);         //extend by dim ´R_sym^{2x2}
    elementMatrix = 0;


    // auto elasticityTensor = material_.getElasticityTensor();

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

    int QPcounter= 0;
    for (const auto& quadPoint : quadRule)
    {
      QPcounter++;
      const auto& quadPos = quadPoint.position();
      const auto geometryJacobianInverse = geometry.jacobianInverse(quadPos);
      const auto integrationElement = geometry.integrationElement(quadPos);

      std::vector<Dune::FieldMatrix<double,1,dim> > jacobians;
      localFiniteElement.localBasis().evaluateJacobian(quadPos, jacobians);

      for (size_t i=0; i< jacobians.size(); i++)
        jacobians[i] = jacobians[i] * geometryJacobianInverse;

      // Compute the scaled deformation gradient for each shape function
      std::vector<std::array<VoigtVector<double,dim>,dimworld> > deformationGradient(nSf);

      for (size_t i=0; i < nSf; i++)
      {
        for (size_t k=0; k< dimworld; k++)
        {
          MatrixRT defGradientV(0);
          defGradientV[k] = jacobians[i][0];

          deformationGradient[i][k] = symVoigt(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
        }
      }

      // Compute the material phase that the current quadrature point is in
      const auto phase = phaseAtQuadPoint[QPcounter-1];

      for (size_t l=0; l< dimworld; l++)
      for (size_t j=0; j < nSf; j++ )
      {
          size_t row = localView.tree().child(l).localIndex(j);
          
          // "phi*phi"-part
          for (size_t k=0; k < dimworld; k++)
          for (size_t i=0; i < nSf; i++)
          {
              // auto etmp = material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)); 
              // auto etmp = elasticityTensor(defGradientU,element.geometry().global(quadPos)); 
              // auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos); 
              // printmatrix(std::cout, etmp, "etmp", "--");
              // double energyDensity= scalarProduct(etmp,defGradientV);
              // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),defGradientV);

              // auto tmp_1 = scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
              // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
              // if (abs(tmp_1-tmp_2)>1e-2)
              // {
              //   std::cout << "difference :" <<  abs(tmp_1-tmp_2) << std::endl;
              //   std::cout << "scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))<< std::endl;
              //   std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)<< std::endl;
              // }
  

              // Python::Module module = Python::import("material_neukamm");
              // auto PythonindicatorFunction = Python::make_function<int>(module.get("f"));
              // if (PythonindicatorFunction(element.geometry().global(quadPos)) != material_.applyIndicatorFunction(element.geometry().global(quadPos)))
              // { 
              //     std::cout <<" element.geometry().global(quadPos): " << element.geometry().global(quadPos) << std::endl;
              //     std::cout << "PythonindicatorFunction(element.geometry().global(quadPos)): " << PythonindicatorFunction(element.geometry().global(quadPos)) << std::endl;
              //     // std::cout << "mu(quadPos):" << mu(quadPos) << std::endl;
              //     std::cout << "indicatorFunction(element.geometry().global(quadPos)): " << indicatorFunction(element.geometry().global(quadPos)) << std::endl;
              //     std::cout << "material_.applyIndicatorFunction(element.geometry().global(quadPos)): " << material_.applyIndicatorFunction(element.geometry().global(quadPos)) << std::endl;
              //     std::cout << "indicatorFunction_material_1: " << indicatorFunction_material_1(element.geometry().global(quadPos)) << std::endl;
              // }
              // double energyDensity= scalarProduct(material_.ElasticityTensorGlobal(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));


            //  std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
            //  std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
              // }


              // std::cout << "--------------- \n";
              // printmatrix(std::cout, defGradientU, "defGradientU", "--");
              // printmatrix(std::cout, material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos))", "--");
              // printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");


              double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
              // double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
              // double energyDensity= scalarProduct(elasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
              // double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
              // double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV));
              // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
              // double energyDensity= scalarProduct(material_.localElasticityTensor(defGradientU,quadPos,element),sym(defGradientV));
              // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,quadPos),sym(defGradientV));
              // double energyDensity= scalarProduct(material_.applyElasticityTensorLocal(defGradientU,quadPos),defGradientV);

            
              // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
  //             double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
  //             double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);  // also works..

              size_t col = localView.tree().child(k).localIndex(i);                       
              
              elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
          }
              
          // "m*phi"  & "phi*m" - part
          for( size_t m=0; m<3; m++)
          {


              // auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) ;
              // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
              // if (abs(tmp_1-tmp_2)>1e-2)
              // {
              //   std::cout << "difference :" <<  abs(tmp_1-tmp_2) << std::endl;
              //   std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
              //   std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
              // }
                

              // double energyDensityGphi = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));





              double energyDensityGphi = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),deformationGradient[j][l]);
              // double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
              // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
              // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) <<std::endl;
              // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
              // double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
              // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
              // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(defGradientV));
              // double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
              // double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(defGradientV));
              // double energyDensityGphi= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(defGradientV));
              // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
              // double energyDensityGphi= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),defGradientV);
              // double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
  //             double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST 
              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++)                //TODO ist symmetric.. reicht die hälfte zu berechnen!!!
        for(size_t n=0; n<3; n++)
        {
            
  //        std::cout << "QPcounter: " << QPcounter << std::endl; 
  //         std::cout << "m= " << m  << "   n = " << n << std::endl;
  //         printmatrix(std::cout, basisContainer[m] , "basisContainer[m]", "--");
  //         printmatrix(std::cout, basisContainer[n] , "basisContainer[n]", "--");
  //         std::cout  << "integrationElement:" << integrationElement << std::endl;
  //         std::cout << "quadPoint.weight(): "<< quadPoint.weight() << std::endl;
  //         std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
  //         std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;

          // auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
          // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
          // if (abs(tmp_1-tmp_2)>1e-2)
          // {
          //     std::cout << "difference :" <<  abs(tmp_1-tmp_2) << std::endl;  
          //     std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])):" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])) << std::endl;
          //     std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])<< std::endl;

          // }
          // double energyDensityGG = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));




          double energyDensityGG = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),matrixBasis_[n]);
          // double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
          // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
          // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
          // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
          // double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
          // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(basisContainer[n]));
          // double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
          // double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(basisContainer[n]));
          // double energyDensityGG= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(basisContainer[n]));
          // double energyDensityGG= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),basisContainer[n]);

          // double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
  //         double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST 
          elementMatrix[localPhiOffset+m][localPhiOffset+n]  += energyDensityGG * quadPoint.weight() * integrationElement;                           // += !!!!! (Fixed-Bug)
          
  //         std::cout  << "energyDensityGG:" << energyDensityGG<< std::endl;
  //         std::cout << "product " << energyDensityGG * quadPoint.weight() * integrationElement << std::endl;
  //         printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
        }
    }
  //   std::cout << "Number of QuadPoints:" << QPcounter << std::endl;
  //   printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
  }



  // Compute the source term for a single element
  // < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
  //   template<class LocalFunction1, class LocalFunction2, class Vector, class LocalForce>
  // void computeElementLoadVector( const typename Basis::LocalView& localView,
  //                               LocalFunction1& mu,
  //                               LocalFunction2& lambda,
  //                               Vector& elementRhs,
  //                               const LocalForce& forceTerm
  //                               )
  template<class Vector, class LocalForce>
  void computeElementLoadVector( const typename Basis::LocalView& localView,
                                Vector& elementRhs,
                                const LocalForce& 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>;

    // auto elasticityTensor = material_.getElasticityTensor();
    // auto indicatorFunction = material_.getIndicatorFunction();
    auto localIndicatorFunction = material_.getLocalIndicatorFunction();
    localIndicatorFunction.bind(element);
    // // auto indicatorFunction = material_.getLocalIndicatorFunction();
    // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
    // auto indicatorFunction = localFunction(indicatorFunctionGVF);
    // indicatorFunction.bind(element);
    // Func2int indicatorFunction = *material_.getIndicatorFunction();


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

  //   int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;  // TEST
  //   int orderQR = 0;
  //   int orderQR = 1;
  //   int orderQR = 2;
  //   int orderQR = 3;
    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);  
    const auto& quad = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
  //   std::cout << "Quad-Rule order used: " << orderQR << std::endl;

    for (const auto& quadPoint : quad)
    {
      const Dune::FieldVector<double,dim>& quadPos = quadPoint.position();
      const auto geometryJacobianInverse = geometry.jacobianInverse(quadPos);
      const double integrationElement = geometry.integrationElement(quadPos);

      std::vector<Dune::FieldMatrix<double,1,dim> > jacobians;
      localFiniteElement.localBasis().evaluateJacobian(quadPos, jacobians);

      for (size_t i=0; i< jacobians.size(); i++)
        jacobians[i] = jacobians[i] * geometryJacobianInverse;
      
      //TEST
  //     std::cout << "forceTerm(element.geometry().global(quadPos)):"  << std::endl;
  //     std::cout << forceTerm(element.geometry().global(quadPos)) << std::endl;
  //     std::cout << "forceTerm(quadPos)" << std::endl;
  //     std::cout << forceTerm(quadPos) << std::endl;
  //     
  //     //TEST QUadrature 
  //     std::cout << "quadPos:" << quadPos << std::endl;
  //     std::cout << "element.geometry().global(quadPos):" << element.geometry().global(quadPos) << std::endl;
  // // //     
  // // //     
  //     std::cout << "quadPoint.weight() :" << quadPoint.weight()  << std::endl;
  //     std::cout << "integrationElement(quadPos):" << integrationElement << std::endl;
  //     std::cout << "mu(quadPos) :" << mu(quadPos) << std::endl;
  //     std::cout << "lambda(quadPos) :" << lambda(quadPos) << std::endl;
      

      // "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 tmpDefGradientV(0);
          tmpDefGradientV[k][0] = jacobians[i][0][0];                                 // Y
          tmpDefGradientV[k][1] = jacobians[i][0][1];                                 // X2
  //         tmpDefGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2];                     //
          tmpDefGradientV[k][2] = jacobians[i][0][2];                     // X3
          
          VoigtVector<double,3> defGradientV = symVoigt(crossSectionDirectionScaling((1.0/gamma_),tmpDefGradientV));

          // double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));



          double energyDensity= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),defGradientV);
          // double energyDensity= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
          // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
          // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
          // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV));
          // double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
          // double energyDensity= scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),sym(defGradientV));
          // double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(defGradientV));
          // double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),defGradientV);
          // double energyDensity= scalarProduct(material_.applyElasticityTensorLocal((-1.0)*forceTerm(quadPos),quadPos),defGradientV);

          // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); 
  //         double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST 
  //         double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
  //         double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST 
          
          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= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));


        double energyDensityfG= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixBasis_[m]);
        // double energyDensityfG= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
        // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
        // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
        // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m]));
        // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
        // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(basisContainer[m]));
        // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),basisContainer[m]);
        // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),basisContainer[m]);
        // double energyDensityfG = scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),basisContainer[m]);
        // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
  //       double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST 
  //       double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
  //       double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST 
        elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
  //       std::cout << "energyDensityfG * quadPoint.weight() * integrationElement: " << energyDensityfG * quadPoint.weight() * integrationElement << std::endl;   
      }
    }
  }


  void assembleCellStiffness(MatrixCT& matrix)
  {
    std::cout << "assemble Stiffness-Matrix begins." << std::endl;

    Dune::MatrixIndexSet occupationPattern;
    getOccupationPattern(occupationPattern);
    occupationPattern.exportIdx(matrix);
    matrix = 0;

    auto localView = basis_.localView();
    const int phiOffset = basis_.dimension();

    // A cache for the element stiffness matrices, if desired
    bool cacheElementMatrices = true;
    std::map<std::vector<int>, ElementMatrixCT> elementMatrixCache;

    // auto muGridF  = makeGridViewFunction(mu_, basis_.gridView());
    // auto muLocal = localFunction(muGridF);
    // auto lambdaGridF  = makeGridViewFunction(lambda_, basis_.gridView());
    // auto lambdaLocal = localFunction(lambdaGridF);

    for (const auto& element : elements(basis_.gridView()))
    {
      localView.bind(element);
      // muLocal.bind(element);
      // lambdaLocal.bind(element);
      const auto localPhiOffset = localView.size();
  //     Dune::Timer Time;
      //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;

      // Determine a suitable quadrature rule
      const auto& localFiniteElement = localView.tree().child(0).finiteElement();
      int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
      const auto& quadRule = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);

      // Determine the material phase at each quadrature point
      auto localIndicatorFunction = material_.getLocalIndicatorFunction();
      localIndicatorFunction.bind(element);

      std::vector<int> phaseAtQuadPoint(quadRule.size());

      for (std::size_t i=0; i<quadRule.size(); ++i)
        phaseAtQuadPoint[i] = localIndicatorFunction(quadRule[i].position());

      ElementMatrixCT elementMatrix;

      // We know that all elements of the grid have the same geometry.  Therefore, the element stiffness matrices
      // depend only on the material phases at the quadrature points.  In practice, a few particular phase distributions
      // dominate (identical phase at all quadrature points).  This makes caching the element matrices worthwhile.
      //
      // In principle the cache implemented here can get arbitrarily large.  If that ever becomes a problem
      // we need to implement a smarter cache.
      if (cacheElementMatrices)
      {
        auto cachedMatrix = elementMatrixCache.find(phaseAtQuadPoint);
        if (cachedMatrix != elementMatrixCache.end())
        {
          // This copies the matrix.  A smarter implementation would avoid this.
          elementMatrix = (*cachedMatrix).second;
        }
        else
        {
          computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
          // This copies the matrix again.
          elementMatrixCache.insert({phaseAtQuadPoint, elementMatrix});
        }
      }
      else // No caching
      {
        computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
      }
      
      
  //     std::cout << "local assembly time:" << Time.elapsed() << std::endl;
      
      //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
      //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
      //std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
      
      //TEST
      //Check Element-Stiffness-Symmetry:
      if(parameterSet_.get<bool>("print_debug", false))
      {
          for (size_t i=0; i<localPhiOffset; i++)
          for (size_t j=0; j<localPhiOffset; j++ )
          {
              if(abs(elementMatrix[i][j] - elementMatrix[j][i]) > 1e-12 )
                  std::cout << "ELEMENT-STIFFNESS MATRIX NOT SYMMETRIC!!!" << 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", "--");
    }
  }


  void assembleCellLoad(VectorCT& b,
                        const Func2Tensor&  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);      

    // auto muGridF  = makeGridViewFunction(mu_, basis_.gridView());
    // auto muLocal = localFunction(muGridF);
    // auto lambdaGridF  = makeGridViewFunction(lambda_, basis_.gridView());
    // auto lambdaLocal = localFunction(lambdaGridF);


  //   int counter = 1;
    for (const auto& element : elements(basis_.gridView()))
    {
      localView.bind(element);
      // muLocal.bind(element);
      // lambdaLocal.bind(element);
      loadFunctional.bind(element);

      const auto localPhiOffset = localView.size();
      //         std::cout << "localPhiOffset : " << localPhiOffset << std::endl;

      VectorCT elementRhs;
  //     std::cout << "----------------------------------Element : " <<  counter << std::endl; //TEST
  //     counter++;
      // computeElementLoadVector(localView, muLocal, lambdaLocal, elementRhs, loadFunctional);
      computeElementLoadVector(localView, elementRhs, loadFunctional);
  //     computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
  //     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", "--");
    }
  }

  // -----------------------------------------------------------------
  // --- Functions for global integral mean equals zero constraint
  auto arbitraryComponentwiseIndices(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();

    Dune::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;
  }

  void setOneBaseFunctionToZero()
  {
    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)
    Dune::FieldVector<std::size_t,3> row = arbitraryComponentwiseIndices(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;
    }
  }


  auto childToIndexMap(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;
  }


  auto integralMean(VectorCT& coeffVector)
  {
    auto GVFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,dim>>(basis_,coeffVector);
    auto localfun = localFunction(GVFunction);

    auto localView = basis_.localView();

    Dune::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 auto& quad = Dune::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 ;
  }


  auto subtractIntegralMean(VectorCT& coeffVector)
  {
    // Subtract correct integral mean from each associated component function
    auto IM = integralMean(coeffVector);

    for(size_t k=0; k<dim; k++)
    {
      //std::cout << "Integral-Mean: " << IM[k] << std::endl;
      auto idx = childToIndexMap(k);
      for ( int i : idx)
        coeffVector[i] -= IM[k];
    }
  }


  // -----------------------------------------------------------------
  // --- Solving the Corrector-problem:
  auto solve()
  {
      bool set_oneBasisFunction_Zero = parameterSet_.get<bool>("set_oneBasisFunction_Zero", false);
      bool substract_integralMean = false;
      if(parameterSet_.get<bool>("set_IntegralZero", false))
      {
          set_oneBasisFunction_Zero = true;
          substract_integralMean = true;
      }
      // set one basis-function to zero
      if(set_oneBasisFunction_Zero)
          setOneBaseFunctionToZero();

      //TEST: Compute Condition Number  (Can be very expensive !) 
      const bool verbose = true;
      const unsigned int arppp_a_verbosity_level = 2;
      const unsigned int pia_verbosity_level = 1;
      if(parameterSet_.get<bool>("print_conditionNumber", false))
      {
        MatrixInfo<MatrixCT> matrixInfo(stiffnessMatrix_,verbose,arppp_a_verbosity_level,pia_verbosity_level);
        std::cout << "Get condition number of Stiffness_CE: " << matrixInfo.getCond2(true) << std::endl;
      }
      ///////////////////////////////////
      // --- Choose Solver ---
      // 1 : CG-Solver
      // 2 : GMRES
      // 3 : QR (default)
      // 4 : UMFPACK
      ///////////////////////////////////
      unsigned int Solvertype = parameterSet_.get<unsigned int>("Solvertype", 3);
      unsigned int Solver_verbosity = parameterSet_.get<unsigned int>("Solver_verbosity", 2);

      // --- set initial values for solver
      x_1_ = load_alpha1_;
      x_2_ = load_alpha2_;
      x_3_ = load_alpha3_;

      std::cout << "start corrector solver..." << std::endl;
      Dune::Timer SolverTimer;
      if (Solvertype==1)  // CG - SOLVER
      {
          std::cout << "------------ CG - Solver ------------" << std::endl;
          Dune::MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);

          // Sequential incomplete LU decomposition as the preconditioner
          Dune::SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_,1.0);
          int iter = parameterSet_.get<double>("iterations_CG", 999);
          // Preconditioned conjugate-gradient solver
          Dune::CGSolver<VectorCT> solver(op,
                                  ilu0, //NULL,
                                  1e-8, // desired residual reduction factorlack
                                  iter,   // maximum number of iterations
                                  Solver_verbosity,
                                  true    // Try to estimate condition_number
                                  );   // verbosity of the solver
          Dune::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";