Forked from
Klaus Böhnlein / dune-microstructure
474 commits behind, 269 commits ahead of the upstream repository.
-
Klaus Böhnlein authoredKlaus Böhnlein authored
EffectiveQuantitiesComputer.hh 20.30 KiB
#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 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 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 Func2int = std::function< int(const Domain&) >;
using MatrixPhase = std::function< MatrixRT(const int&) >;
protected:
CorrectorComputer<Basis,Material>& correctorComputer_;
// Func2Tensor prestrain_;
// MatrixPhase prestrain_;
const Material& material_;
// Func2int indicatorFunction_;
public:
MatrixRT Qeff_; // Effective moduli
VectorRT Bhat_;
VectorRT Beff_; // Effective prestrain
///////////////////////////////
// constructor
///////////////////////////////
EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
const Material& material)
: correctorComputer_(correctorComputer),
material_(material)
{
// prestrain_ = material_.getPrestrainFunction();
}
///////////////////////////////
// Getter
///////////////////////////////
CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}
const shared_ptr<Basis> getBasis() {return correctorComputer_.getBasis();}
auto getQeff(){return Qeff_;} auto getBeff(){return Beff_;}
// -----------------------------------------------------------------
// --- Compute Effective Quantities
void computeEffectiveQuantities()
{
std::cout << "starting effective quantities computation..." << std::endl;
Dune::Timer effectiveQuantitiesTimer;
// 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();
Dune::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);
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
Qeff_ = 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(Dune::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& element : elements(basis.gridView()))
{
localView.bind(element);
matrixFieldG1.bind(element);
matrixFieldG2.bind(element);
localfun_a.bind(element);
// DerPhi2.bind(e);
// mu.bind(e);
// lambda.bind(e); // prestrainFunctional.bind(e);
localIndicatorFunction.bind(element);
double elementEnergy = 0.0;
double elementPrestrain = 0.0;
auto geometry = element.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 auto& quad = Dune::QuadratureRules<double, dim>::rule(element.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 = matrixToVoigt(G1 + Chi1);
// auto X2 = G2 + Chi2;
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)), matrixToVoigt(sym(G2)));
// double energyDensity= scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
if (b==0)
{
// std::cout << "WORKS TILL HERE" << std::endl;
// auto prestrainPhaseValue_old = prestrain_(localIndicatorFunction(quadPos));
auto prestrainPhaseValue = material_.prestrain(localIndicatorFunction(quadPos),element.geometry().global(quadPos));
// if (localIndicatorFunction(quadPos) != 3)
// {
// std::cout << "element.geometry().global(quadPos):"<< element.geometry().global(quadPos) << std::endl;
// printmatrix(std::cout, prestrainPhaseValue_old, "prestrainPhaseValue_old", "--");
// printmatrix(std::cout, prestrainPhaseValue, "prestrainPhaseValue", "--");
// }
auto value = voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),matrixToVoigt(sym(prestrainPhaseValue)));
elementPrestrain += value * quadPoint.weight() * integrationElement;
// elementPrestrain += scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue)) * quadPoint.weight() * integrationElement;
// elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos)) * quadPoint.weight() * integrationElement;
}
}
energy += elementEnergy;
prestrain += elementPrestrain;
}
Qeff_[a][b] = energy;
if (b==0)
Bhat_[a] = prestrain;
}
if(parameterSet.get<bool>("print_debug", false))
{
printmatrix(std::cout, Qeff_, "Qeff_", "--");
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(Qeff_,true,2,1);
// std::cout << "----------------------------------------" << std::endl;
Qeff_.solve(Beff_,Bhat_);
if(parameterSet.get<bool>("print_debug", false))
printvector(std::cout, Beff_, "Beff_", "--");
//LOG-Output
auto& log = *(correctorComputer_.getLog());
log << "--- Effective moduli --- " << std::endl;
log << "Qeff_: \n" << Qeff_ << std::endl;
log << "------------------------ " << std::endl;
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;
std::cout << "Time required to solve for effective quantities: " << effectiveQuantitiesTimer.elapsed() << std::endl;
return ;
}
// -----------------------------------------------------------------
// --- write Data to Matlab / Optimization-Code
void writeToMatlab(std::string outputPath)
{
std::cout << "write effective quantities as .txt files to output folder..." << std::endl;
//writeMatrixToMatlab(Q, "../../Matlab-Programs/QMatrix.txt");
writeMatrixToMatlab(Qeff_, outputPath + "/QMatrix.txt");
// write effective Prestrain in Matrix for Output
Dune::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 localIndicatorFunction = material_.getLocalIndicatorFunction();
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);
localIndicatorFunction.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 auto& quad = Dune::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= scalarProduct(material_.applyElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
// double energyDensity= scalarProduct(material_.ElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(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);
// Qeff_ = 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;
// }
// Qeff_[a][b] = energy;
// if (b==0)
// Bhat_[a] = prestrain;
// }
// printmatrix(std::cout, Qeff_, "Matrix Qeff_", "--");
// 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(Qeff_,true,2,1);
// // std::cout << "----------------------------------------" << std::endl;
// Qeff_.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