#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