#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/common/parametertree.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/microstructure/matrix_operations.hh>
// #include <dune/microstructure/materialDefinitions.hh> //deprecated
using namespace MatrixOperations;
using namespace Dune::Functions;
using namespace Dune::Functions::BasisFactory;
using std::pow;
using std::abs;
using std::sqrt;
using std::sin;
using std::cos;
using std::shared_ptr;
using std::make_shared;
// // ----------------------------------------------------------------------------------
// template<class Domain>
// int indicatorFunction_material(const Domain& x)
// {
// double theta=0.25;
// if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
// return 1; //#Phase1
// else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
// return 2; //#Phase2
// else
// return 3; //#Phase3
// }
// MatrixRT getL(const int& phase)
// {
// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
// if (phase == 1)
// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); //#Phase1 //Isotrop(mu,lambda)
// else if (phase == 2)
// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
// else
// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id(); //#Phase3
// }
// MatrixRT prestrain_material(const int& phase)
// {
// if (phase == 1)
// return {{1.0, 0.0, 0.0}, //#Phase1
// {0.0, 1.0, 0.0},
// {0.0, 0.0, 1.0}
// };
// else if (phase == 2)
// return {{1.0, 0.0, 0.0}, //#Phase2
// {0.0, 1.0, 0.0},
// {0.0, 0.0, 1.0}
// };
// else
// return {{0.0, 0.0, 0.0}, //#Phase3
// {0.0, 0.0, 0.0},
// {0.0, 0.0, 0.0}
// };
// }
//--------------------------------------------------------
template <class GridView>
class prestrainedMaterial
{
public:
static const int dimworld = 3; //GridView::dimensionworld;
static const int dim = 3; //Domain::dimension;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using LocalDomain = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
using Element = typename GridViewEntitySet<GridView, 0>::Element;
using ScalarRT = Dune::FieldVector< double, 1>;
using VectorRT = Dune::FieldVector< double, dimworld>;
using MatrixRT = Dune::FieldMatrix< double, dimworld, dimworld>;
using FuncScalar = std::function< double(const Domain&) >;
using Func2int = std::function< int(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
using Func2TensorPhase = std::function< MatrixRT(const MatrixRT& ,const int&) >;
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
using MatrixPhase = std::function< MatrixRT(const int&) >;
// using VoigtMatrix = FieldMatrix< double, 6, 6>;
protected:
const GridView& gridView_;
const Dune::ParameterTree& parameterSet_;
std::string materialFunctionName_;
// MatrixFunc L1_;
// MatrixFunc L2_;
// MatrixFunc L3_;
// Elasticity tensors (in voigt notation)
Dune::FieldMatrix<double,6,6> L1_;
Dune::FieldMatrix<double,6,6> L2_;
Dune::FieldMatrix<double,6,6> L3_;
Dune::FieldMatrix<double,6,6> L4_; //not used yet
Func2Tensor prestrain1_;
Func2Tensor prestrain2_;
Func2Tensor prestrain3_;
Func2Tensor prestrain4_; //not used yet
Python::Module module_;
// const FuncScalar mu_;
// const FuncScalar lambda_;
// const int phases_;
// Func2Tensor elasticityTensor_;
// Func2TensorParam elasticityTensor_;
Func2TensorPhase elasticityTensor_;
MatrixPhase prestrain_;
Func2int indicatorFunction_;
GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
LocalFunction<int(const Domain&), typename GridViewEntitySet<GridView, 0>::Element > localIndicatorFunction_;
//Transformation matrix for Voigt notation
Dune::FieldMatrix<double,6,6> D_ = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
{0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
{0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
{0.0, 0.0, 0.0, sqrt(2.0) , 0.0, 0.0},
{0.0, 0.0, 0.0, 0.0 , sqrt(2.0), 0.0},
{0.0, 0.0, 0.0, 0.0 , 0.0, sqrt(2.0)}
};
Dune::FieldMatrix<double,6,6> D_inv = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
{0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
{0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
{0.0, 0.0, 0.0, 1.0/sqrt(2.0) , 0.0, 0.0},
{0.0, 0.0, 0.0, 0.0 , 1.0/sqrt(2.0), 0.0},
{0.0, 0.0, 0.0, 0.0 , 0.0, 1.0/sqrt(2.0)}
};
public:
///////////////////////////////
// constructor
///////////////////////////////
prestrainedMaterial(const GridView gridView,
const Dune::ParameterTree& parameterSet)
: gridView_(gridView),
parameterSet_(parameterSet)
{
std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
// setupMaterial(materialFunctionName_); //old-Version
module_ = Python::import(materialFunctionName_); //TODO use this instead?
// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
//elastic moduli in the order (E1,E2,E3,G12,G23,G31,nu12,nu13,nu23)
//-- Walnut see [Dimwoodie; Timber its nature and behavior p.109]
// const FieldVector<double,9> walnut_parameters={11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37};
//-- Norway spruce see [Dimwoodie; Timber its nature and behavior p.109]
// const FieldVector<double,9> Nspruce_parameters={10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31};
//frame vectors
VectorRT e1 = {1.0,0.0,0.0};
VectorRT e2 = {0.0,1.0,0.0};
VectorRT e3 = {0.0,0.0,1.0};
// L1_ = orthotropic(e1,e2,e3,walnut_parameters);
// L2_ = orthotropic(e1,e2,e3,Nspruce_parameters);
// L3_ = orthotropic(e1,e2,e3,Nspruce_parameters);
setup(materialFunctionName_);
// setup("material");
indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
// L1_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
// L2_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63, 0.37});
// L3_ = transversely_isotropic(e1,e2,e3, {10.7e3,710 ,500,0.51 ,0.31});
// L1_ = isotropic(mu[0],lambda[0]);
// L2_ = isotropic(mu[1],lambda[1]);
// L3_ = isotropic(mu[2],lambda[2]);
// L1_ = {{lambda[0]+2.0*mu[0], lambda[0], lambda[0], 0.0 , 0.0, 0.0},
// {lambda[0], lambda[0]+2*mu[0], lambda[0], 0.0 ,0.0 ,0.0},
// {lambda[0], lambda[0], lambda[0]+2.0*mu[0], 0.0, 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 2.0*mu[0], 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[0], 0.0},
// { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[0]}
// };
// L2_ = {{lambda[1]+2.0*mu[1], lambda[1], lambda[1], 0.0 , 0.0, 0.0},
// {lambda[1], lambda[1]+2*mu[1], lambda[1], 0.0 ,0.0 ,0.0},
// {lambda[1], lambda[1], lambda[1]+2.0*mu[1], 0.0, 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 2.0*mu[1], 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[1], 0.0},
// { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[1]}
// };
// L3_ = {{lambda[2]+2.0*mu[2], lambda[2], lambda[2], 0.0 , 0.0, 0.0},
// {lambda[2], lambda[2]+2.0*mu[2], lambda[2], 0.0 ,0.0 ,0.0},
// {lambda[2], lambda[2], lambda[2]+2.0*mu[2], 0.0, 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 2.0*mu[2], 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[2], 0.0},
// { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[2]}
// };
// printmatrix(std::cout, D_, "D_", "--");
// printmatrix(std::cout, D_inv, "D_inv", "--");
// printmatrix(std::cout, L1_, "L1_", "--");
// L1_.leftmultiply(D_inv);
// printmatrix(std::cout, L1_, "L1_.leftmultiply(D_inv)", "--");
// L1_.rightmultiply(D_);
// printmatrix(std::cout, L1_, "L1_.rightmultiply(D_)", "--");
// Python::Module module = Python::import(materialFunctionName_);
// elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
// elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
// auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
// module.get("Phases").toC<int>(Phases_);
// L1_ = Python::make_function<MatrixRT>(module.get("L1"));
// L2_ = Python::make_function<MatrixRT>(module.get("L2"));
// L3_ = Python::make_function<MatrixRT>(module.get("L3"));
// bool isotropic_ = true; // read from module File
// auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
}
Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
{
std::string materialParameters_string;
std::cout << "Setup phase "<< phase << " as " << phaseType << " material." << std::endl;
switch (phase)
{
case 1:
{
materialParameters_string = "materialParameters_phase1";
break;
}
case 2:
{
materialParameters_string = "materialParameters_phase2";
break;
}
case 3:
{
materialParameters_string= "materialParameters_phase3";
break;
}
default:
DUNE_THROW(Dune::Exception, "Phase number not implemented.");
break;
}
if(phaseType == "isotropic") //TODO SetupPHASE (phase_type)
{
Dune::FieldVector<double,2> materialParameters(0);
module.get(materialParameters_string).toC<Dune::FieldVector<double,2>>(materialParameters);
return isotropic(materialParameters[0],materialParameters[1]);
}
if(phaseType == "orthotropic") //TODO SetupPHASE (phase_type)
{
// TODO: Generalize here:
//frame vectors
VectorRT e1 = {1.0,0.0,0.0};
VectorRT e2 = {0.0,1.0,0.0};
VectorRT e3 = {0.0,0.0,1.0};
Dune::FieldVector<double,9> materialParameters(0);
module.get(materialParameters_string).toC<Dune::FieldVector<double,9>>(materialParameters);
return orthotropic(e1,e2,e3,materialParameters);
}
if(phaseType == "transversely_isotropic") //TODO SetupPHASE (phase_type)
{
// TODO: Generalize here:
//frame vectors
VectorRT e1 = {1.0,0.0,0.0};
VectorRT e2 = {0.0,1.0,0.0};
VectorRT e3 = {0.0,0.0,1.0};
Dune::FieldVector<double,5> materialParameters(0);
module.get(materialParameters_string).toC<Dune::FieldVector<double,5>>(materialParameters);
return transversely_isotropic(e1,e2,e3,materialParameters);
}
if(phaseType == "general_anisotropic") //TODO SetupPHASE (phase_type)
{
// TODO: Generalize here:
//frame vectors
VectorRT e1 = {1.0,0.0,0.0};
VectorRT e2 = {0.0,1.0,0.0};
VectorRT e3 = {0.0,0.0,1.0};
Dune::FieldMatrix<double,6,6> materialParameters(0); //compliance matric
module.get(materialParameters_string).toC<Dune::FieldMatrix<double,6,6>>(materialParameters);
printmatrix(std::cout, materialParameters, "Compliance matrix used:", "--");
return general_anisotropic(e1,e2,e3,materialParameters);
}
else
DUNE_THROW(Dune::Exception, " Unknown material class for phase ");
}
//---function that determines elasticity Tensor
// void setup(const std::string name, const int Phases) // cant use materialFunctionName_ here!?
void setup(const std::string name) // cant use materialFunctionName_ here!?
{
Python::Module module = Python::import(name);
indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
std::cout << "---- setup material... " << std::endl;
int Phases;
module.get("Phases").toC<int>(Phases);
std::cout << "Number of Phases: " << Phases << std::endl;
switch (Phases)
{
case 1:
{
std::string phase1_type;
module.get("phase1_type").toC<std::string>(phase1_type);
L1_ = setupPhase(phase1_type, module,1);
prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
break;
}
case 2:
{
std::string phase1_type;
module.get("phase1_type").toC<std::string>(phase1_type);
std::string phase2_type;
module.get("phase2_type").toC<std::string>(phase2_type);
L1_ = setupPhase(phase1_type, module,1);
L2_ = setupPhase(phase2_type, module,2);
prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
break;
}
case 3:
{
std::string phase1_type;
module.get("phase1_type").toC<std::string>(phase1_type);
std::string phase2_type;
module.get("phase2_type").toC<std::string>(phase2_type);
std::string phase3_type;
module.get("phase3_type").toC<std::string>(phase3_type);
L1_ = setupPhase(phase1_type, module,1);
L2_ = setupPhase(phase2_type, module,2);
L3_ = setupPhase(phase3_type, module,3);
prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
prestrain3_ = Python::make_function<MatrixRT>(module.get("prestrain_phase3"));
// if (phase1_type == "isotropic" ) //TODO SetupPHASE (phase_type)
// {
// FieldVector<double,2> materialParameters_phase1(0);
// module.get("materialParameters_phase1").toC<FieldVector<double,2>>(materialParameters_phase1);
// L1_ = isotropic(materialParameters_phase1[0],materialParameters_phase1[1]);
// }
// if (phase2_type == "isotropic" ) //TODO SetupPHASE (phase_type)
// {
// FieldVector<double,2> materialParameters_phase2(0);
// module.get("materialParameters_phase2").toC<FieldVector<double,2>>(materialParameters_phase2);
// L2_ = isotropic(materialParameters_phase2[0],materialParameters_phase2[1]);
// }
// if (phase3_type == "isotropic" ) //TODO SetupPHASE (phase_type)
// {
// FieldVector<double,2> materialParameters_phase3(0);
// module.get("materialParameters_phase3").toC<FieldVector<double,2>>(materialParameters_phase3);
// L3_ = isotropic(materialParameters_phase3[0],materialParameters_phase3[1]);
// }
break;
}
default:
break;
}
std::cout << "Material setup done." << std::endl;
}
///////////////////////////////////////////////////////////////////////
// //---function that determines elasticity Tensor
// void setupMaterial(const std::string name) // cant use materialFunctionName_ here!?
// {
// // std::cout << "Using material definition:" << name << std::endl;
// if(name == "material_neukamm")
// {
// indicatorFunction_ = indicatorFunction_material_neukamm<Domain>;
// elasticityTensor_ = material_neukamm;
// prestrain_ = prestrain_material_neukamm;
// }
// else if(name == "two_phase_material_1")
// {
// indicatorFunction_ = indicatorFunction_two_phase_material_1<Domain>;
// elasticityTensor_ = two_phase_material_1;
// prestrain_ = prestrain_two_phase_material_1;
// }
// else if(name == "two_phase_material_2")
// {
// indicatorFunction_ = indicatorFunction_two_phase_material_2<Domain>;
// elasticityTensor_ = two_phase_material_2;
// prestrain_ = prestrain_two_phase_material_2;
// }
// else if(name == "homogeneous")
// {
// indicatorFunction_ = indicatorFunction_material_homogeneous<Domain>;
// elasticityTensor_ = material_homogeneous;
// prestrain_ = prestrain_homogeneous;
// }
// else
// DUNE_THROW(Exception, "There exists no material with this name ");
// }
////////////////////////////////////////////////////////
// MATERIAL CLASSES DEFINITIONS:
////////////////////////////////////////////////////////
//--- Definition of isotropic elasticity tensor (in voigt notation)
Dune::FieldMatrix<double,6,6> isotropic(const double& mu, const double& lambda)
{
return {{lambda+2.0*mu, lambda , lambda , 0.0 , 0.0 , 0.0 },
{lambda , lambda+2.0*mu, lambda , 0.0 , 0.0 , 0.0 },
{lambda , lambda , lambda+2.0*mu, 0.0 , 0.0 , 0.0 },
{ 0.0 , 0.0 , 0.0 , 2.0*mu, 0.0 , 0.0 },
{ 0.0 , 0.0 , 0.0 , 0.0 , 2.0*mu, 0.0 },
{ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 2.0*mu}
};
}
//-------------------------------------------------------------------------------
//--- Definition of orthotropic elasticity tensor (in voigt notation)
// - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
// - Output: compliance Matrix
Dune::FieldMatrix<double,6,6> orthotropic(const VectorRT& framevector1,
const VectorRT& framevector2,
const VectorRT& framevector3,
const Dune::FieldVector<double,9>& p //elastic moduli {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
)
{
// (see https://en.wikipedia.org/wiki/Orthotropic_material)
auto E_1 = p[0];
auto E_2 = p[1];
auto E_3 = p[2];
auto G_12 = p[3];
auto G_23 = p[4];
auto G_31 = p[5];
auto nu_12 = p[6];
auto nu_13 = p[7];
auto nu_23 = p[8];
//other three poisson ratios are not independent
auto nu_21 = (p[6]/p[0])*p[1];
auto nu_31 = (p[7]/p[0])*p[2];
auto nu_32 = (p[8]/p[1])*p[2];
//compliance matrix
Dune::FieldMatrix<double,6,6> C = {{1.0/E_1 , -nu_21/E_2 , -nu_31/E_3 , 0.0 , 0.0 , 0.0 },
{-nu_12/E_1 , 1.0/E_2 , -nu_32/E_3 , 0.0 , 0.0 , 0.0 },
{-nu_13/E_1 , -nu_23/E_2 , 1.0/E_3 , 0.0 , 0.0 , 0.0 },
{ 0.0 , 0.0 , 0.0 , 1.0/G_23, 0.0 , 0.0 },
{ 0.0 , 0.0 , 0.0 , 0.0 , 1.0/G_31, 0.0 },
{ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 1.0/G_12}
};
// printmatrix(std::cout, C, "C", "--");
//--- return elasticity tensor
C.invert();
// printmatrix(std::cout, C, "C after .invert()", "--");
return C;
}
//-------------------------------------------------------------------------------
//--- Definition of transversely isotropic elasticity tensor (in voigt notation)
// - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
// - Output: compliance Matrix
Dune::FieldMatrix<double,6,6> transversely_isotropic(const VectorRT& framevector1,
const VectorRT& framevector2,
const VectorRT& framevector3,
const Dune::FieldVector<double,5>& p //elastic moduli {E1,E2,G12,nu12,nu23}
)
{
// (see https://en.wikipedia.org/wiki/Poisson%27s_ratio)
auto E_1 = p[0];
auto E_2 = p[1];
auto E_3 = E_2;
auto nu_12 = p[3];
auto nu_13 = nu_12;
auto nu_21 = (nu_12/E_1)*E_2;
auto nu_31 = nu_21;
auto nu_23 = p[4];
auto nu_32 = nu_23;
auto G_12 = p[2];
auto G_23 = E_2/(2.0*(1+nu_23));
auto G_31 = G_23;
//other three poisson ratios are not independent
//---compliance matrix
Dune::FieldMatrix<double,6,6> C = {{1.0/E_1 , -nu_21/E_2 , -nu_31/E_3 , 0.0 , 0.0 , 0.0 },
{-nu_12/E_1 , 1.0/E_2 , -nu_32/E_3 , 0.0 , 0.0 , 0.0 },
{-nu_13/E_1 , -nu_23/E_2 , 1.0/E_3 , 0.0 , 0.0 , 0.0 },
{ 0.0 , 0.0 , 0.0 , 1.0/G_23, 0.0 , 0.0 },
{ 0.0 , 0.0 , 0.0 , 0.0 , 1.0/G_31, 0.0 },
{ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 1.0/G_12}
};
// printmatrix(std::cout, C, "C", "--");
//--- return elasticity tensor
C.invert();
// printmatrix(std::cout, C, "C after .invert()", "--");
return C;
}
//-------------------------------------------------------------------------------
// --- Definition of transversely isotropic elasticity tensor (in voigt notation)
// - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
// - Output: inverse compliance Matrix
Dune::FieldMatrix<double,6,6> general_anisotropic(const VectorRT& framevector1,
const VectorRT& framevector2,
const VectorRT& framevector3,
const Dune::FieldMatrix<double,6,6>& C //compliance matrix
)
{
//--- return elasticity tensor (in Voigt notation)
auto tmp = C;
tmp.invert();
// printmatrix(std::cout, C, "C after .invert()", "--");
return tmp;
}
//-------------------------------------------------------------------------------
// --- Helpers
static Dune::FieldVector<double,6> MatrixToVoigt(const MatrixRT& G)
{
// return {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
return {G[0][0], G[1][1], G[2][2], 2.0*G[1][2], 2.0*G[0][2], 2.0*G[0][1]};
}
static MatrixRT VoigtToMatrix(const Dune::FieldVector<double,6>& v)
{
// return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
return {{v[0], v[5]/2.0, v[4]/2.0}, {v[5]/2.0, v[1], v[3]/2.0}, {v[4]/2.0, v[3]/2.0, v[2]}};
}
MatrixRT applyElasticityTensor(const MatrixRT& G, const int& phase) const
{
// FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], sqrt(2.0)*G[1][2], sqrt(2.0)*G[0][2], sqrt(2.0)*G[0][1]};
// FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
Dune::FieldVector<double,6> G_tmp = MatrixToVoigt(G);
Dune::FieldVector<double,6> out(0);
// printvector(std::cout, G_tmp, "G_tmp", "--");
if (phase == 1)
{
// printmatrix(std::cout, L1_, "L1_", "--");
L1_.mv(G_tmp,out);
}
else if (phase == 2)
{
// printmatrix(std::cout, L2_, "L2_", "--");
L2_.mv(G_tmp,out);
}
else
{
// printmatrix(std::cout, L3_, "L3_", "--");
L3_.mv(G_tmp,out);
}
// printvector(std::cout, out, "out", "--");
return VoigtToMatrix(out);
// return image as Matrix again
// return {{out[0], (1.0/sqrt(2.0))*out[5], (1.0/sqrt(2.0))*out[4]}, {(1.0/sqrt(2.0))*out[5], out[1], (1.0/sqrt(2.0))*out[3]}, {(1.0/sqrt(2.0))*out[4], (1.0/sqrt(2.0))*out[3], out[2]}};
// return {{out[0], (1.0/2.0)*out[5], (1.0/2.0)*out[4]}, {(1.0/2.0)*out[5], out[1], (1.0/2.0)*out[3]}, {(1.0/2.0)*out[4], (1.0/2.0)*out[3], out[2]}};
// return {{out[0], out[5], out[4]}, {out[5], out[1], out[3]}, {out[4], out[3], out[2]}};
}
////////////////////////////////////////////////////////////////////////////////////////
MatrixRT prestrain(const int& phase,const Domain& x) const
{
if (phase == 1)
return prestrain1_(x);
else if (phase == 2)
return prestrain2_(x);
else if (phase ==3)
return prestrain3_(x);
else
DUNE_THROW(Dune::Exception, "Phase number not implemented.");
}
////////////////////////////////////////////////////////
//--- Apply elasticityTensor_ to input Matrix G at material phase
// input: Matrix G, material-phase
MatrixRT ElasticityTensor(const MatrixRT& G, const int& phase) const //deprecated
{
return elasticityTensor_(G,phase);
}
int IndicatorFunction(const Domain& x) const
{
return indicatorFunction_(x);
}
///////////////////////////////
// Getter
///////////////////////////////
Dune::ParameterTree getParameterSet() const {return parameterSet_;}
// Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
// Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
MatrixPhase getPrestrainFunction() const {return prestrain_;}
auto getLocalIndicatorFunction() const {return localIndicatorFunction_;} //get as localFunction
auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;} //get as GridViewFunction
auto getIndicatorFunction() const {return indicatorFunction_;}
//getLameParameters()? .. Throw Exception if isotropic_ = false!
///////////////////////////////
// VTK - Writer
///////////////////////////////
// -----------------------------------------------------------------
// --- write material (grid)functions to VTK
void writeVTKMaterialFunctions(std::array<int, dim> nElements, const int level)
{
std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
using VTKGridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
// VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
using GridViewVTK = typename VTKGridType::LeafGridView;
const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
std::vector<double> indicatorFunction_CoeffP0;
Dune::Functions::interpolate(scalarP0FeBasis, indicatorFunction_CoeffP0, indicatorFunction_);
auto indicatorFunction_P0 = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, indicatorFunction_CoeffP0);
std::vector<double> indicatorFunction_CoeffP1;
Dune::Functions::interpolate(scalarP1FeBasis, indicatorFunction_CoeffP1, indicatorFunction_);
auto indicatorFunction_P1 = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, indicatorFunction_CoeffP1);
Dune::VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
MaterialVtkWriter.addCellData(
indicatorFunction_P0,
Dune::VTK::FieldInfo("indicatorFunction_P0", Dune::VTK::FieldInfo::Type::scalar, 1));
MaterialVtkWriter.addVertexData(
indicatorFunction_P1,
Dune::VTK::FieldInfo("indicatorFunction_P1", Dune::VTK::FieldInfo::Type::scalar, 1));
MaterialVtkWriter.write(outputPath + "/MaterialFunctions-level"+ std::to_string(level) );
std::cout << "wrote data to file:" + outputPath +"/MaterialFunctions-level" + std::to_string(level) << std::endl;
return;
};
// -----------------------------------------------------------------
// --- write prestrain (grid)functions to VTK
// void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
// {
// std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
// using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
// //VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
// VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
// using GridViewVTK = VTKGridType::LeafGridView;
// const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
// auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
// auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
// std::vector<double> prestrainFunction_CoeffP0;
// Functions::interpolate(scalarP0FeBasis, prestrainFunction_CoeffP0, prestrain_);
// auto prestrainFunction_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, prestrainFunction_CoeffP0);
// std::vector<double> prestrainFunction_CoeffP1;
// Functions::interpolate(scalarP1FeBasis, prestrainFunction_CoeffP1, prestrain_);
// auto prestrainFunction_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, prestrainFunction_CoeffP1);
// VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
// MaterialVtkWriter.addVertexData(
// prestrainFunction_P0,
// VTK::FieldInfo("prestrainFunction_P0", VTK::FieldInfo::Type::scalar, 1));
// MaterialVtkWriter.addCellData(
// prestrainFunction_P1,
// VTK::FieldInfo("prestrainFunction_P1", VTK::FieldInfo::Type::scalar, 1));
// MaterialVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
// std::cout << "wrote data to file:" + outputPath + "/PrestrainFunctions-level" + std::to_string(level) << std::endl;
// return;
// };
// void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
// {
// using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
// // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
// // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{40,40,40});
// VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
// using GridViewVTK = VTKGridType::LeafGridView;
// const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
// FTKfillerContainer<dim> VTKFiller;
// VTKFiller.vtkPrestrainNorm(gridView_VTK, B_Term, "PrestrainBNorm");
// // WORKS Too
// VTKFiller.vtkProblemCell(gridView_VTK, B_Term, muLocal,"VTKProblemCell");;
// // TEST
// auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
// auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
// std::vector<double> B_CoeffP0;
// Functions::interpolate(scalarP0FeBasis, B_CoeffP0, B_Term);
// auto B_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, B_CoeffP0);
// VTKWriter<GridView> PrestrainVtkWriter(gridView_VTK);
// PrestrainVtkWriter.addCellData(
// B_DGBF_P0,
// VTK::FieldInfo("B_P0", VTK::FieldInfo::Type::scalar, 1));
// PrestrainVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
// std::cout << "wrote data to file:" + outputPath +"/PrestrainFunctions-level" + std::to_string(level) << std::endl;
// // };
// MatrixRT ElasticityTensorGlobal(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
// {
// return MAT(G,x);
// }
// //--- apply elasticityTensor_ to input Matrix G at position x
// MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
// {
// // Python::Module module = Python::import("material");
// // auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
// // return elasticityTensor_(G,indicatorFunctionTest(x));
// // auto IFunction = localFunction(indicatorFunction_);
// // auto IFunction = this->getIndicatorFunction();
// // auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// // auto tmpLocal = LocalF
// // return elasticityTensor_(G,IFunction(x));
// // return elasticityTensor_(G,localIndicatorFunction_(x));
// return elasticityTensor_(G,indicatorFunction_(x));
// // return elasticityTensor_(G,indicatorFunctionGVF_(x));
// // int phase = indicatorFunction_(x);
// // auto tmp = elasticityTensor_(G,phase);
// // auto tmp = material_1(G,indicatorFunction_(x));
// // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
// // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
// }
// template<class Domain>
// static MatrixRT MAT(const MatrixRT& G, const Domain& x) //unfortunately not possible as a GridViewFunction?
// {
// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
// double theta=0.25;
// if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
// else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
// else
// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
// }
}; // end class
#endif