From e70be8b72f6664ca575864efdae6ad08bccbb32e Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 31 Aug 2022 12:02:33 +0200
Subject: [PATCH] Add setup method for: isotropic, orthotropic,
transversely_isotropic
---
dune/microstructure/CorrectorComputer.hh | 10 +-
dune/microstructure/prestrainedMaterial.hh | 212 ++++++++++++++++++---
inputs/cellsolver.parset | 2 +-
3 files changed, 188 insertions(+), 36 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index f5a7e7ef..d9edddfd 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -371,10 +371,12 @@ public:
// std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
// }
- std::cout << "--------------- \n";
- printmatrix(std::cout, defGradientU, "defGradientU", "--");
- printmatrix(std::cout, material_.applyL(defGradientU,localIndicatorFunction(quadPos)), "material_.applyL(defGradientU,localIndicatorFunction(quadPos))", "--");
- printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
+
+
+ // std::cout << "--------------- \n";
+ // printmatrix(std::cout, defGradientU, "defGradientU", "--");
+ // printmatrix(std::cout, material_.applyL(defGradientU,localIndicatorFunction(quadPos)), "material_.applyL(defGradientU,localIndicatorFunction(quadPos))", "--");
+ // printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
double energyDensity= scalarProduct(material_.applyL(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 5bfb83ed..4c9e489d 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -166,29 +166,61 @@ public:
const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
- 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]}
- };
+
+ //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);
+
+
+ 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, L1_, "L1_", "--");
@@ -247,13 +279,131 @@ public:
//////////// TESTING //////////////////////////////////////
+
+ static 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]};
+ }
+
+ static MatrixRT VoigtToMatrix(const FieldVector<double,6>& v)
+ {
+ return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
+ }
+
+
+
+ //--- Definition of isotropic elasticity tensor (in voigt notation)
+ 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
+ FieldMatrix<double,6,6> orthotropic(const VectorRT& framevector1,
+ const VectorRT& framevector2,
+ const VectorRT& framevector3,
+ const 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
+ 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
+ FieldMatrix<double,6,6> transversely_isotropic(const VectorRT& framevector1,
+ const VectorRT& framevector2,
+ const VectorRT& framevector3,
+ const 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
+ 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;
+ }
+
+
+
+
+
+
+
+
+
MatrixRT applyL(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]};
+ // FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
+ FieldVector<double,6> G_tmp = MatrixToVoigt(G);
FieldVector<double,6> out(0);
- printvector(std::cout, G_tmp, "G_tmp", "--");
+ // printvector(std::cout, G_tmp, "G_tmp", "--");
//
// printmatrix(std::cout, L2_, "L2_", "--");
@@ -261,27 +411,27 @@ public:
if (phase == 1)
{
- printmatrix(std::cout, L1_, "L1_", "--");
+ // printmatrix(std::cout, L1_, "L1_", "--");
L1_.mv(G_tmp,out);
}
else if (phase == 2)
{
- printmatrix(std::cout, L2_, "L2_", "--");
+ // printmatrix(std::cout, L2_, "L2_", "--");
L2_.mv(G_tmp,out);
}
else
{
- printmatrix(std::cout, L3_, "L3_", "--");
+ // printmatrix(std::cout, L3_, "L3_", "--");
L3_.mv(G_tmp,out);
}
- printvector(std::cout, out, "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]}};
+ // return {{out[0], out[5], out[4]}, {out[5], out[1], out[3]}, {out[4], out[3], out[2]}};
}
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 163ae6c4..88cfed92 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -13,7 +13,7 @@ outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
# --- DEBUG (Output) Option:
-#print_debug = true #(default=false)
+print_debug = true #(default=false)
#############################################
--
GitLab