From b371988e3bcdbe0d18c1ee0737dba5fa7203304d Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sun, 1 Oct 2023 19:38:03 +0200
Subject: [PATCH] Do not hard-wire the number of phases
There is no need for this. At least, I don't see any.
---
dune/microstructure/prestrainedMaterial.hh | 197 ++++-----------------
1 file changed, 33 insertions(+), 164 deletions(-)
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index e2dcfc70..a35aa2fe 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -104,20 +104,10 @@ protected:
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_;
+ std::vector<VoigtMatrix<double,3> > L_;
- Func2Tensor prestrain1_;
- Func2Tensor prestrain2_;
- Func2Tensor prestrain3_;
- Func2Tensor prestrain4_;
+ std::vector<Func2Tensor> prestrain_;
Python::Module module_;
@@ -127,7 +117,6 @@ protected:
// Func2Tensor elasticityTensor_;
// Func2TensorParam elasticityTensor_;
Func2TensorPhase elasticityTensor_;
- MatrixPhase prestrain_;
Func2int indicatorFunction_;
@@ -186,9 +175,9 @@ public:
- // L1_ = orthotropic(e1,e2,e3,walnut_parameters);
- // L2_ = orthotropic(e1,e2,e3,Nspruce_parameters);
- // L3_ = orthotropic(e1,e2,e3,Nspruce_parameters);
+ // L_[0] = orthotropic(e1,e2,e3,walnut_parameters);
+ // L_[1] = orthotropic(e1,e2,e3,Nspruce_parameters);
+ // L_[2] = orthotropic(e1,e2,e3,Nspruce_parameters);
setup(materialFunctionName_);
// setup("material");
@@ -198,18 +187,18 @@ public:
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});
+ // L_[0] = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
+ // L_[1] = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63, 0.37});
+ // L_[2] = 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]);
+ // L_[0] = isotropic(mu[0],lambda[0]);
+ // L_[1] = isotropic(mu[1],lambda[1]);
+ // L_[2] = isotropic(mu[2],lambda[2]);
- // L1_ = {{lambda[0]+2.0*mu[0], lambda[0], lambda[0], 0.0 , 0.0, 0.0},
+ // L_[0] = {{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},
@@ -217,7 +206,7 @@ public:
// { 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},
+ // L_[1] = {{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},
@@ -225,7 +214,7 @@ public:
// { 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},
+ // L_[2] = {{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},
@@ -235,11 +224,11 @@ public:
// 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_)", "--");
+ // printmatrix(std::cout, L_[0], "L_[0]", "--");
+ // L_[0].leftmultiply(D_inv);
+ // printmatrix(std::cout, L_[0], "L_[0].leftmultiply(D_inv)", "--");
+ // L_[0].rightmultiply(D_);
+ // printmatrix(std::cout, L_[0], "L_[0].rightmultiply(D_)", "--");
// Python::Module module = Python::import(materialFunctionName_);
@@ -249,9 +238,9 @@ public:
// 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"));
+ // L_[0] = Python::make_function<MatrixRT>(module.get("L1"));
+ // L_[1] = Python::make_function<MatrixRT>(module.get("L2"));
+ // L_[2] = Python::make_function<MatrixRT>(module.get("L3"));
// bool isotropic_ = true; // read from module File
// auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
}
@@ -389,81 +378,15 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
module.get("Phases").toC<int>(phases);
std::cout << "---- Starting material setup... (Number of Phases: "<< phases << ") " << std::endl;
- switch (phases)
- {
- case 1:
- {
- // std::string phase1_type;
- // module.get("phase1_type").toC<std::string>(phase1_type);
+ L_.resize(phases);
+ prestrain_.resize(phases);
- L1_ = setupPhase(module,1);
- // prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
- prestrain1_ = setupPrestrainPhase(module,1);
- 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(module,1);
- L2_ = setupPhase(module,2);
- // prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
- // prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
- prestrain1_ = setupPrestrainPhase(module,1);
- prestrain2_ = setupPrestrainPhase(module,2);
- 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(module,1);
- L2_ = setupPhase(module,2);
- L3_ = setupPhase(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"));
- prestrain1_ = setupPrestrainPhase(module,1);
- prestrain2_ = setupPrestrainPhase(module,2);
- prestrain3_ = setupPrestrainPhase(module,3);
- break;
- }
- case 4:
- {
- // 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);
- // std::string phase4_type;
- // module.get("phase4_type").toC<std::string>(phase4_type);
-
- L1_ = setupPhase(module,1);
- L2_ = setupPhase(module,2);
- L3_ = setupPhase(module,3);
- L4_ = setupPhase(module,4);
- // 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"));
- // prestrain4_ = Python::make_function<MatrixRT>(module.get("prestrain_phase4"));
- prestrain1_ = setupPrestrainPhase(module,1);
- prestrain2_ = setupPrestrainPhase(module,2);
- prestrain3_ = setupPrestrainPhase(module,3);
- prestrain4_ = setupPrestrainPhase(module,4);
- break;
- }
- default:
- DUNE_THROW(Dune::Exception, " No more than 4 material phases are implemented yet! ");
+ for (int i=0; i<phases; i++)
+ {
+ L_[i] = setupPhase(module,i+1);
+ prestrain_[i] = setupPrestrainPhase(module, i+1);
}
+
std::cout << "Material setup done." << std::endl;
}
@@ -679,71 +602,17 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
VoigtVector<double,3> G_tmp = strainToVoigt(G);
VoigtVector<double,3> out(0);
- switch (phase)
- {
- case 1:
- {
- // printmatrix(std::cout, L1_, "L1_", "--");
- L1_.mv(G_tmp,out);
- break;
- }
- case 2:
- {
- // printmatrix(std::cout, L2_, "L2_", "--");
- L2_.mv(G_tmp,out);
- break;
- }
- case 3:
- {
- // printmatrix(std::cout, L3_, "L3_", "--");
- L3_.mv(G_tmp,out);
- break;
- }
- case 4:
- {
- // printmatrix(std::cout, L4_, "L4_", "--");
- L4_.mv(G_tmp,out);
- break;
- }
- default:
- DUNE_THROW(Dune::Exception, "Phase number not implemented.");
- break;
- }
+ // printmatrix(std::cout, L_[0], "L_[0]", "--");
+ L_[phase-1].mv(G_tmp,out);
return voigtToStress(out);
-
- }
+ }
////////////////////////////////////////////////////////////////////////////////////////
MatrixRT prestrain(const int& phase,const Domain& x) const
{
- switch (phase)
- {
- case 1:
- {
- return prestrain1_(x);
- break;
- }
- case 2:
- {
- return prestrain2_(x);
- break;
- }
- case 3:
- {
- return prestrain3_(x);
- break;
- }
- case 4:
- {
- return prestrain3_(x);
- break;
- }
- default:
- DUNE_THROW(Dune::Exception, "Phase number not implemented.");
- break;
- }
+ return prestrain_[phase-1](x);
}
--
GitLab