Newer
Older
#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;
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
// // ----------------------------------------------------------------------------------
// 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>
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 VoigtMatrix = FieldMatrix< double, 6, 6>;
const Dune::ParameterTree& parameterSet_;
// MatrixFunc L1_;
// MatrixFunc L2_;
// MatrixFunc L3_;
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_;
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});
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
// 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"));
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
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;
}
///////////////////////////////////////////////////////////////////////
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
// //---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}
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
// printmatrix(std::cout, C, "C after .invert()", "--");
return tmp;
}
//-------------------------------------------------------------------------------
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", "--");
// 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_", "--");
// 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);
DUNE_THROW(Dune::Exception, "Phase number not implemented.");
////////////////////////////////////////////////////////
//--- Apply elasticityTensor_ to input Matrix G at material phase
MatrixRT ElasticityTensor(const MatrixRT& G, const int& phase) const //deprecated
}
int IndicatorFunction(const Domain& x) const
{
return indicatorFunction_(x);
}
///////////////////////////////
///////////////////////////////
Dune::ParameterTree getParameterSet() const {return parameterSet_;}
// Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
// Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
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);
Dune::VTK::FieldInfo("indicatorFunction_P0", Dune::VTK::FieldInfo::Type::scalar, 1));
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;
};
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
// -----------------------------------------------------------------
// --- 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;
// // };
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
// 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();
// }