Newer
Older
#ifndef DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
#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>
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&) >;
CorrectorComputer<Basis,Material>& correctorComputer_;
const Material& material_;
///////////////////////////////
// constructor
///////////////////////////////
EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
const Material& material)
: correctorComputer_(correctorComputer),
material_(material)
CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}
const shared_ptr<Basis> getBasis() {return correctorComputer_.getBasis();}
// -----------------------------------------------------------------
// --- Compute Effective Quantities
void computeEffectiveQuantities()
{
std::cout << "starting effective quantities computation..." << std::endl;
// 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();
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);
for (const auto& element : elements(basis.gridView()))
localView.bind(element);
matrixFieldG1.bind(element);
matrixFieldG2.bind(element);
localfun_a.bind(element);
double elementEnergy = 0.0;
double elementPrestrain = 0.0;
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);
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;
}
///////////////////////////////
// Compute effective Prestrain B_eff (by solving linear system)
//////////////////////////////
// std::cout << "------- Information about Q matrix -----" << std::endl; // TODO
// std::cout << "----------------------------------------" << std::endl;
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;
// -----------------------------------------------------------------
// --- 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");
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));
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
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);
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
// 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;
// }
// 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;
// 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 ;
// }