EffectiveQuantitiesComputer.hh

#ifndef DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH

#include <filesystem>


#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/CorrectorComputer.hh>

#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;
using std::shared_ptr;
using std::make_shared;
using std::string;
using std::cout;
using std::endl;

// template <class Basis>
// class EffectiveQuantitiesComputer : public CorrectorComputer<Basis,Material> {

template <class Basis, class Material>
class EffectiveQuantitiesComputer {

public:
	static const int dimworld = 3;
	// static const int nCells = 4;
    
	static const int dim = Basis::GridView::dimension;

	using Domain = typename CorrectorComputer<Basis,Material>::Domain; 

	using VectorRT = typename CorrectorComputer<Basis,Material>::VectorRT;
	using MatrixRT = typename CorrectorComputer<Basis,Material>::MatrixRT;

	using Func2Tensor = typename CorrectorComputer<Basis,Material>::Func2Tensor;
	using FuncVector = typename CorrectorComputer<Basis,Material>::FuncVector;

	using VectorCT = typename CorrectorComputer<Basis,Material>::VectorCT;

	using HierarchicVectorView = typename CorrectorComputer<Basis,Material>::HierarchicVectorView;

protected:

	CorrectorComputer<Basis,Material>& correctorComputer_; 
	Func2Tensor prestrain_;
    const Material& material_;

public:
	VectorCT B_load_TorusCV_;				//<B, Chi>_L2 
	// FieldMatrix<double, dim, dim> Q_;	    //effective moduli <LF_i, F_j>_L2
	// FieldVector<double, dim> Bhat_;			//effective loads induced by prestrain <LF_i, B>_L2
	// FieldVector<double, dim> Beff_;		//effective strains Mb = ak
	MatrixRT Q_;	    //effective moduli <LF_i, F_j>_L2
	VectorRT Bhat_;			//effective loads induced by prestrain <LF_i, B>_L2
	VectorRT Beff_;		//effective strains Mb = ak


	// corrector parts
	VectorCT phi_E_TorusCV_;		//phi_i * (a,K)_i
	VectorCT phi_perp_TorusCV_;
	VectorCT phi_TorusCV_;
	VectorCT phi_1_;		//phi_i * (a,K)_i
	VectorCT phi_2_;
	VectorCT phi_3_;
	
	// is this really interesting???
	// double phi_E_L2norm_;
	// double phi_E_H1seminorm_;

	// double phi_perp_L2norm_;
	// double phi_perp_H1seminorm_;

	// double phi_L2norm_;
	// double phi_H1seminorm_;

	// double Chi_E_L2norm_;		
	// double Chi_perp_L2norm_;
	// double Chi_L2norm_;


	double B_energy_;			 // < B, B >_L 		B = F + Chi_perp + B_perp 
	double F_energy_;			 // < F, F >_L
	double Chi_perp_energy_;	 // < Chi_perp, Chi_perp >_L
	double B_perp_energy_; 		 // < B_perp, B_perp >_L

	//Chi(phi) is only implicit computed, can we store this?

 
	///////////////////////////////
	// constructor
	///////////////////////////////
	// EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer, Func2Tensor prestrain)
    //     : correctorComputer_(correctorComputer), prestrain_(prestrain)
	EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer, 
                                Func2Tensor prestrain,
                                const Material& material)
        : correctorComputer_(correctorComputer), 
          prestrain_(prestrain),
          material_(material)
    { 
    	
    	// computePrestressLoadCV();
	  	// computeEffectiveStrains();
        // Q_ = 0;
        // Q_ = {{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
    	// compute_phi_E_TorusCV();
    	// compute_phi_perp_TorusCV();
    	// compute_phi_TorusCV();

    	// computeCorrectorNorms();
    	// computeChiNorms();
    	// computeEnergiesPrestainParts();	

    	// writeInLogfile();
    } 


    ///////////////////////////////
    // getter
    ///////////////////////////////
	CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}

	const shared_ptr<Basis> getBasis()  
	{
		return correctorComputer_.getBasis();
	}

    auto getQeff(){return Q_;}
    auto getBeff(){return Beff_;}


  // -----------------------------------------------------------------
  // --- Compute Effective Quantities
    void computeEffectiveQuantities()
    {

        // Get everything.. better TODO: with Inheritance?
        // auto test = correctorComputer_.getLoad_alpha1();
        // auto phiContainer = correctorComputer_.getPhicontainer();
        auto MContainer = correctorComputer_.getMcontainer();
        auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
        auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
        auto mu_ = *correctorComputer_.getMu();
        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(),
                                            correctorComputer_.getCorr_phi3()
										    };

        auto prestrainGVF  = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
        auto prestrainFunctional = localFunction(prestrainGVF);   

        Q_ = 0 ;
        Bhat_ = 0;
    
        for(size_t a=0; a < 3; a++)
        for(size_t b=0; b < 3; b++)
        {
            double energy = 0.0;
            double prestrain = 0.0;
            auto localView = basis.localView();
            // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
            auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
            //   auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,phiContainer[b]));
            auto localfun_a = localFunction(GVFunc_a);
            //   auto localfun_b = localFunction(GVFunc_b);

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

            auto matrixFieldG1GVF  = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[a], basis.gridView());
            auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
            auto matrixFieldG2GVF  = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[b], basis.gridView());
            auto matrixFieldG2 = localFunction(matrixFieldG2GVF);

            auto muGridF  = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
            auto mu = localFunction(muGridF);
            auto lambdaGridF  = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
            auto lambda= localFunction(lambdaGridF);

            // using GridView = typename Basis::GridView;

            for (const auto& e : elements(basis.gridView()))
            {
                localView.bind(e);
                matrixFieldG1.bind(e);
                matrixFieldG2.bind(e);
                localfun_a.bind(e);
                // DerPhi2.bind(e);
                mu.bind(e);
                lambda.bind(e);
                prestrainFunctional.bind(e);

                double elementEnergy = 0.0;
                double elementPrestrain = 0.0;

                auto geometry = e.geometry();
                const auto& localFiniteElement = localView.tree().child(0).finiteElement();

            //     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 );  // TEST
                int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
            //     int orderQR = 0;
            //     int orderQR = 1;
            //     int orderQR = 2;
            //     int orderQR = 3;
                const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);

                for (const auto& quadPoint : quad) 
                {
                    const auto& quadPos = quadPoint.position();
                    const double integrationElement = geometry.integrationElement(quadPos);
                    
                    auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a];
                    
                    
                    auto G1 = matrixFieldG1(quadPos);
                    auto G2 = matrixFieldG2(quadPos);
                //       auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
                //       auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
                    
                    auto X1 = G1 + Chi1;
                    //   auto X2 = G2 + Chi2;
                    
                    
                    double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
                    elementEnergy += energyDensity * quadPoint.weight() * integrationElement;      // quad[quadPoint].weight() ???
                    if (b==0)
                    {
                        elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos)) * quadPoint.weight() * integrationElement;
                    }
                }
                energy += elementEnergy;
                prestrain += elementPrestrain;
            
            }
            Q_[a][b] = energy;    
            if (b==0)
                Bhat_[a] = prestrain;
        }
        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)
        //////////////////////////////
         
        // std::cout << "------- Information about Q matrix -----" << std::endl;        // TODO
        // MatrixInfo<MatrixRT> matrixInfo(Q_,true,2,1);
        // std::cout << "----------------------------------------" << std::endl;
        Q_.solve(Beff_,Bhat_);
        if(parameterSet.get<bool>("print_debug", false))
            printvector(std::cout, Beff_, "Beff_", "--");
        
    
        //LOG-Output
        auto& log = *(correctorComputer_.getLog());
        log << "--- Prestrain Output --- " << std::endl;
        log << "Bhat_: " << Bhat_ << std::endl;
        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 ;
    }


  // -----------------------------------------------------------------
  // --- write Data to Matlab / Optimization-Code
    void writeToMatlab(std::string outputPath)
    {
        std::cout << "write effective quantities to Matlab folder..." << std::endl;
        //writeMatrixToMatlab(Q, "../../Matlab-Programs/QMatrix.txt");
        writeMatrixToMatlab(Q_, outputPath + "/QMatrix.txt");
        // write effective Prestrain in Matrix for Output
        FieldMatrix<double,1,3> BeffMat;
        BeffMat[0] = Beff_;
        writeMatrixToMatlab(BeffMat, outputPath + "/BMatrix.txt");
        return;
    }




    template<class MatrixFunction>
    double energySP(const MatrixFunction& matrixFieldFuncA,
                    const MatrixFunction& matrixFieldFuncB)
    {
        double energy = 0.0;
        auto mu_ = *correctorComputer_.getMu();
        auto lambda_ = *correctorComputer_.getLambda();
        auto gamma = correctorComputer_.getGamma();
        auto basis = *correctorComputer_.getBasis();
        auto localView = basis.localView();

        auto matrixFieldAGVF  = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
        auto matrixFieldA = localFunction(matrixFieldAGVF);
        auto matrixFieldBGVF  = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
        auto matrixFieldB = localFunction(matrixFieldBGVF);
        auto muGridF  = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
        auto mu = localFunction(muGridF);
        auto lambdaGridF  = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
        auto lambda= localFunction(lambdaGridF);
        for (const auto& e : elements(basis.gridView()))
        {
            localView.bind(e);
            matrixFieldA.bind(e);
            matrixFieldB.bind(e);
            mu.bind(e);
            lambda.bind(e);

            double elementEnergy = 0.0;

            auto geometry = e.geometry();
            const auto& localFiniteElement = localView.tree().child(0).finiteElement();

            int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
            const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
            for (const auto& quadPoint : quad) 
            {
                const auto& quadPos = quadPoint.position();
                const double integrationElement = geometry.integrationElement(quadPos);
                double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), matrixFieldA(quadPos), matrixFieldB(quadPos));
                elementEnergy += energyDensity * quadPoint.weight() * integrationElement;          
            }
            energy += elementEnergy;
        }
        return energy;
    }



    // --- Alternative that does not use orthogonality relation (75) in the paper
    // void computeFullQ()
    // {
    //     auto MContainer = correctorComputer_.getMcontainer();
    //     auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
    //     auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
    //     auto mu_ = *correctorComputer_.getMu();
    //     auto lambda_ = *correctorComputer_.getLambda();
    //     auto gamma = correctorComputer_.getGamma();
    //     auto basis = *correctorComputer_.getBasis();

	// 	shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(), 
    //                                         correctorComputer_.getCorr_phi2(),
    //                                         correctorComputer_.getCorr_phi3()
	// 									    };

    //     auto prestrainGVF  = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
    //     auto prestrainFunctional = localFunction(prestrainGVF);   

    //     Q_ = 0 ;
    //     Bhat_ = 0;
    
    //     for(size_t a=0; a < 3; a++)
    //     for(size_t b=0; b < 3; b++)
    //     {
    //         double energy = 0.0;
    //         double prestrain = 0.0;
    //         auto localView = basis.localView();
    //         // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
    //         auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
    //         auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[b]));
    //         auto localfun_a = localFunction(GVFunc_a);
    //         auto localfun_b = localFunction(GVFunc_b);

    //         ///////////////////////////////////////////////////////////////////////////////
    //         auto matrixFieldG1GVF  = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[a], basis.gridView());
    //         auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
    //         auto matrixFieldG2GVF  = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[b], basis.gridView());
    //         auto matrixFieldG2 = localFunction(matrixFieldG2GVF);

    //         auto muGridF  = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
    //         auto mu = localFunction(muGridF);
    //         auto lambdaGridF  = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
    //         auto lambda= localFunction(lambdaGridF);

    //         // using GridView = typename Basis::GridView;

    //         for (const auto& e : elements(basis.gridView()))
    //         {
    //             localView.bind(e);
    //             matrixFieldG1.bind(e);
    //             matrixFieldG2.bind(e);
    //             localfun_a.bind(e);
    //             localfun_b.bind(e);
    //             mu.bind(e);
    //             lambda.bind(e);
    //             prestrainFunctional.bind(e);

    //             double elementEnergy = 0.0;
    //             double elementPrestrain = 0.0;

    //             auto geometry = e.geometry();
    //             const auto& localFiniteElement = localView.tree().child(0).finiteElement();

    //         //     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 );  // TEST
    //             int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
    //         //     int orderQR = 0;
    //         //     int orderQR = 1;
    //         //     int orderQR = 2;
    //         //     int orderQR = 3;
    //             const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);

    //             for (const auto& quadPoint : quad) 
    //             {
    //                 const auto& quadPos = quadPoint.position();
    //                 const double integrationElement = geometry.integrationElement(quadPos);
                    
    //                 auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a] + matrixFieldG1(quadPos);
    //                 auto Chi2 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_b(quadPos))) + *MContainer[b] + matrixFieldG2(quadPos);
                    
    //                 // auto G1 = matrixFieldG1(quadPos);
    //                 // auto G2 = matrixFieldG2(quadPos);
    //             //       auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
    //             //       auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
                    
    //                 // auto X1 = G1 + Chi1;
    //                 //   auto X2 = G2 + Chi2;
                    
                    
    //                 double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), Chi1, Chi2);
    //                 elementEnergy += energyDensity * quadPoint.weight() * integrationElement;      // quad[quadPoint].weight() ???
    //             }
    //             energy += elementEnergy;
    //             prestrain += elementPrestrain;
            
    //         }
    //         Q_[a][b] = energy;    
    //         if (b==0)
    //             Bhat_[a] = prestrain;
    //     }
    //     printmatrix(std::cout, Q_, "Matrix Q_", "--");
    //     printvector(std::cout, Bhat_, "Bhat_", "--");
    //     ///////////////////////////////
    //     // Compute effective Prestrain B_eff (by solving linear system)
    //     //////////////////////////////
    //     // std::cout << "------- Information about Q matrix -----" << std::endl;        // TODO
    //     // MatrixInfo<MatrixRT> matrixInfo(Q_,true,2,1);
    //     // std::cout << "----------------------------------------" << std::endl;
    //     Q_.solve(Beff_,Bhat_);
    //     printvector(std::cout, Beff_, "Beff_", "--");
        
    //     //LOG-Output
    //     auto& log = *(correctorComputer_.getLog());
    //     log << "--- Prestrain Output --- " << std::endl;
    //     log << "Bhat_: " << Bhat_ << std::endl;
    //     log << "Beff_: " << Beff_ <<  " (Effective Prestrain)" << std::endl;
    //     log << "------------------------ " << std::endl;
    //     return ;
    // }


}; // end class



#endif