From f62c44f5c9253f4f3f59fe598725e0f6c2243913 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sun, 1 Oct 2023 18:41:02 +0200
Subject: [PATCH] Add option to rotate material coordinate frame
---
dune/microstructure/matrix_operations.hh | 58 +++
dune/microstructure/prestrainedMaterial.hh | 494 +++++++++++++--------
materials/material_orthotropic.py | 14 +-
3 files changed, 378 insertions(+), 188 deletions(-)
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index fc1670e5..3eca389f 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -105,6 +105,64 @@ namespace MatrixOperations {
return sum;
}
+
+ /**
+ * @brief Determines rotation matrix based on an axis and an angle
+ *
+ * @param axis
+ * @param angle
+ * @return MatrixRT
+ */
+ static MatrixRT rotationMatrix(int axis, double angle){
+
+ switch (axis)
+ {
+ case 0:
+ {
+ return {{1.0, 0, 0 },
+ { 0, cos(angle), -1.0*sin(angle)},
+ { 0, sin(angle), cos(angle)}};
+ }
+ case 1:
+ {
+ return {{ cos(angle), 0, sin(angle)},
+ { 0, 1.0, 0 },
+ {-1.0*sin(angle), 0, cos(angle)}};
+ }
+ case 2:
+ {
+ return {{cos(angle), -1.0*sin(angle), 0},
+ {sin(angle), cos(angle), 0},
+ { 0, 0, 1.0}};
+ }
+ default:
+ DUNE_THROW(Dune::Exception, " axis not feasible. rotationMatrix is only implemented for 3x3-matrices.");
+ }
+ }
+
+ /**
+ * @brief 6x6 matrix that transforms the strain tensor. This matrix is used to
+ * transform the compliance matrix (given in Voigt notation) into another frame.
+ * see 'https://en.wikipedia.org/wiki/Orthotropic_material#Condition_for_material_symmetry_2'
+ * for details.
+ *
+ * @param axis
+ * @param angle
+ * @return MatrixRT
+ */
+ static Dune::FieldMatrix<double,6,6> rotationMatrixCompliance(int axis, double angle){
+
+ MatrixRT R = rotationMatrix(axis,angle);
+
+ return {{ R[0][0]*R[0][0], R[0][1]*R[0][1], R[0][2]*R[0][2], R[0][1]*R[0][2], R[0][0]*R[0][2], R[0][0]*R[0][1]},
+ { R[1][0]*R[1][0], R[1][1]*R[1][1], R[1][2]*R[1][2], R[1][1]*R[1][2], R[1][0]*R[1][2], R[1][0]*R[1][1]},
+ { R[2][0]*R[2][0], R[2][1]*R[2][1], R[2][2]*R[2][2], R[2][1]*R[2][2], R[2][0]*R[2][2], R[2][0]*R[2][1]},
+ {2.0*R[1][0]*R[2][0], 2.0*R[1][1]*R[2][1], 2.0*R[1][2]*R[2][2], R[1][1]*R[2][2]+R[1][2]*R[2][1], R[1][0]*R[2][2]+R[1][2]*R[2][0], R[1][0]*R[2][1]+R[1][1]*R[2][0]},
+ {2.0*R[0][0]*R[2][0], 2.0*R[0][1]*R[2][1], 2.0*R[0][2]*R[2][2], R[0][1]*R[2][2]+R[0][2]*R[2][1], R[0][0]*R[2][2]+R[0][2]*R[2][0], R[0][0]*R[2][1]+R[0][1]*R[2][0]},
+ {2.0*R[0][0]*R[0][0], 2.0*R[0][1]*R[1][1], 2.0*R[0][2]*R[1][2], R[0][1]*R[1][2]+R[0][2]*R[1][1], R[0][0]*R[1][2]+R[0][2]*R[1][0], R[0][0]*R[1][1]+R[0][1]*R[1][0]}
+ };
+ }
+
// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){ // ?? Check with Robert
// E1= sym(E1);
// E2 = sym(E2);
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 60ddc691..8651815a 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -185,10 +185,6 @@ public:
//-- 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};
@@ -262,172 +258,213 @@ public:
// auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
}
+/**
+ * @brief TODO: Code allows a spatially varying prestrain. However most of the time we only
+ * use a constant prestrain.
+ * - Add switch (constant/non-constant) prestrain?
+ *
+ * @param module python module used
+ * @param phase material phase
+ * @return Func2Tensor
+ */
+Func2Tensor setupPrestrainPhase(Python::Module module, int phase){
+
+ auto prestrain = Python::make_function<MatrixRT>(module.get("prestrain_phase" + std::to_string(phase)));
+
+ int axis = 0;
+ double angle = 0;
+ try
+ {
+ module.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
+ module.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
+ }
+ catch(Dune::Exception)
+ {
+ //default frame is used.
+ }
+ /**
+ * @brief (optional) Tramsform prestrain to the correct frame.
+ */
+ if(abs(angle) > 1e-12){
+ auto R = rotationMatrix(axis,angle);
+ Func2Tensor f = [prestrain,R] (const Domain& x) { return (prestrain(x).leftmultiply(R)).rightmultiply(R.transposed()); };
+
+ return f;
+ }
+ else
+ return prestrain;
+}
-
-Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
+// Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
+Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
{
- std::string materialParameters_string;
- std::cout << "Setup phase "<< phase << " as " << phaseType << " material." << std::endl;
+ std::string phaseType;
+ module.get("phase" + std::to_string(phase) + "_type").toC<std::string>(phaseType);
- switch (phase)
- {
- case 1:
- {
- materialParameters_string = "materialParameters_phase1";
- break;
- }
- case 2:
- {
- materialParameters_string = "materialParameters_phase2";
- break;
- }
- case 3:
- {
- materialParameters_string= "materialParameters_phase3";
- break;
- }
- case 4:
- {
- materialParameters_string= "materialParameters_phase4";
- break;
- }
- default:
- DUNE_THROW(Dune::Exception, "Phase number not implemented.");
- break;
- }
+ std::cout << "Setup phase "<< phase << " as " << phaseType << " material." << std::endl;
- if(phaseType == "isotropic") //TODO SetupPHASE (phase_type)
+ if(phaseType == "isotropic")
{
Dune::FieldVector<double,2> materialParameters(0);
- module.get(materialParameters_string).toC<Dune::FieldVector<double,2>>(materialParameters);
+ module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,2>>(materialParameters);
std::cout << "Lame-Parameters (mu,lambda): (" << materialParameters[0] << "," << materialParameters[1] << ") " << std::endl;
return isotropic(materialParameters[0],materialParameters[1]);
}
- if(phaseType == "orthotropic") //TODO SetupPHASE (phase_type)
+
+ /**
+ * @brief For other material classes. The stiffness-/compliance matrix depends on the
+ * coordinate frame.
+ *
+ * The default frame vectors:
+ */
+ // VectorRT frameVector1 = {1.0,0.0,0.0};
+ // VectorRT frameVector2 = {0.0,1.0,0.0};
+ // VectorRT frameVector3 = {0.0,0.0,1.0};
+ /**
+ * @brief For other material classes. The stiffness-/compliance matrix depends on the
+ * coordinate frame.
+ *
+ * (optional) rotate material coordinate frame by an 'angle' around 'axis'
+ */
+ int axis = 0;
+ double angle = 0;
+
+ try
+ {
+ // module.get("phase" + std::to_string(phase) + "_FrameVector1").toC<VectorRT>(frameVector1);
+ // module.get("phase" + std::to_string(phase) + "_FrameVector2").toC<VectorRT>(frameVector2);
+ // module.get("phase" + std::to_string(phase) + "_FrameVector3").toC<VectorRT>(frameVector3);
+ // printvector(std::cout, frameVector1, "frameVector1: ", "--");
+ // printvector(std::cout, frameVector2, "frameVector2: ", "--");
+ // printvector(std::cout, frameVector3, "frameVector3: ", "--");
+ module.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
+ module.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
+ }
+ catch(Dune::Exception)
+ {
+ // std::cout << "(Could not read frame vectors for material phase) " << phase << ": - Canonical Frame (default) is used." << std::endl;
+ std::cout << "(Could not read rotaton axis & angle for material phase) " << phase << ": - Canonical Frame (default) is used." << std::endl;
+ }
+
+
+ if(phaseType == "orthotropic")
{
- // 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);
+ module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,9>>(materialParameters);
+ // return orthotropic(frameVector1,frameVector2,frameVector3,materialParameters);
+ return orthotropic(axis,angle,materialParameters);
}
- if(phaseType == "transversely_isotropic") //TODO SetupPHASE (phase_type)
+ if(phaseType == "transversely_isotropic")
{
- // 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);
+ module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,5>>(materialParameters);
+ // return transversely_isotropic(frameVector1,frameVector2,frameVector3,materialParameters);
+ return transversely_isotropic(axis,angle,materialParameters);
}
- if(phaseType == "general_anisotropic") //TODO SetupPHASE (phase_type)
+ if(phaseType == "general_anisotropic")
{
- // 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);
+ module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldMatrix<double,6,6>>(materialParameters);
+ // return general_anisotropic(frameVector1,frameVector2,frameVector3,materialParameters);
+ return general_anisotropic(axis,angle,materialParameters);
}
else
- DUNE_THROW(Dune::Exception, " Unknown material class for phase ");
+ DUNE_THROW(Dune::Exception, " Unknown material class for phase " + std::to_string(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!?
+ /**
+ * @brief Setup method that determines the correct stiffness-matrix (constitutive law)
+ * for each material phase separately.
+ *
+ * @param name name of the material file being used.
+ */
+ 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;
+ int phases;
+ module.get("Phases").toC<int>(phases);
+ std::cout << "---- Starting material setup... (Number of Phases: "<< phases << ") " << std::endl;
-
- switch (Phases)
- {
- case 1:
+ switch (phases)
{
- 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"));
+ case 1:
+ {
+ // std::string phase1_type;
+ // module.get("phase1_type").toC<std::string>(phase1_type);
- 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"));
-
- 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(phase1_type, module,1);
- L2_ = setupPhase(phase2_type, module,2);
- L3_ = setupPhase(phase3_type, module,3);
- L4_ = setupPhase(phase4_type, 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"));
-
- break;
- }
- default:
- break;
+ 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! ");
}
std::cout << "Material setup done." << std::endl;
}
@@ -440,26 +477,41 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
//--- Definition of isotropic elasticity tensor (stiffness matrix 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 , mu , 0.0 , 0.0 },
- { 0.0 , 0.0 , 0.0 , 0.0 , mu , 0.0 },
- { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , mu }
- };
+ 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 , mu , 0.0 , 0.0 },
+ { 0.0 , 0.0 , 0.0 , 0.0 , mu , 0.0 },
+ { 0.0 , 0.0 , 0.0 , 0.0 , 0.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}
- )
+ // 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}
+ // )
+
+ /**
+ * @brief Definition of orthotropic elasticity/stiffness tensor (in Voigt notation).
+ * The convention in literature is to provide the elastic constants in form of the 'compliance matrix'
+ * In order to compute stresses from strain, we have to invert the 'compliance matrix' to obtain the 'stiffness matrix'
+ * (see https://en.wikipedia.org/wiki/Orthotropic_material for more details)
+ *
+ * @param axis rotate material frame around axis
+ * @param angle rotation angle
+ * @param p elastic constants: (E1,E2,E3,G12,G23,G31,nu12,nu13,nu23)
+ * @return Stiffness matrix
+ */
+ Dune::FieldMatrix<double,6,6> orthotropic(const int& axis,
+ const double& angle,
+ const Dune::FieldVector<double,9>& p
+ )
{
- // (see https://en.wikipedia.org/wiki/Orthotropic_material)
auto E_1 = p[0];
auto E_2 = p[1];
auto E_3 = p[2];
@@ -474,34 +526,58 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
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}
- };
+ std::string axis_string;
+ switch (axis)
+ {
+ case 0: axis_string = "x - axis";
+ case 1: axis_string = "y - axis";
+ case 2: axis_string = "z - axis";
+ default: break;
+ }
- printmatrix(std::cout, C, "(Orthotropic) Compliance matrix used:", "--");
- // printmatrix(std::cout, C, "C", "--");
+ // Compliance matrix
+ Dune::FieldMatrix<double,6,6> Compliance = {{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}
+ };
+
+
+ /**
+ * @brief (optional) Tramsform compliance matrix to the correct frame.
+ */
+ if(abs(angle) > 1e-12){
+ // std::cout << "material frame rotated around axis: " << axis << " by an angle of: " << angle << std::endl;
+ std::cout << "material frame rotated around axis: " << axis_string << " by an angle of: " << angle << std::endl;
+ Dune::FieldMatrix<double,6,6> C_rot = rotationMatrixCompliance(axis,angle);
+ (Compliance.leftmultiply(C_rot)).rightmultiply(C_rot.transposed());
+ }
- //--- return elasticity tensor
- C.invert();
- // printmatrix(std::cout, C, "C after .invert()", "--");
- return C;
+
+ printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
+
+ // invert to obtain elasticity/stiffness-Tensor
+ Compliance.invert();
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
+ return Compliance;
}
//-------------------------------------------------------------------------------
+ // 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}
+ // )
//--- 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}
- )
+ Dune::FieldMatrix<double,6,6> transversely_isotropic(const int& axis,
+ const double& angle,
+ const Dune::FieldVector<double,5>& p
+ )
{
// (see https://en.wikipedia.org/wiki/Poisson%27s_ratio)
auto E_1 = p[0];
@@ -517,22 +593,42 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
auto G_23 = E_2/(2.0*(1+nu_23));
auto G_31 = G_23;
//other three poisson ratios are not independent
+
+ std::string axis_string;
+ switch (axis)
+ {
+ case 0: axis_string = "x - axis";
+ case 1: axis_string = "y - axis";
+ case 2: axis_string = "z - axis";
+ default: break;
+ }
//---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}
- };
+ Dune::FieldMatrix<double,6,6> Compliance = {{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}
+ };
+
+ /**
+ * @brief (optional) Tramsform compliance matrix to the correct frame.
+ */
+ if(abs(angle) > 1e-12){
+ // std::cout << "material frame rotated around axis: " << axis << " by an angle of: " << angle << std::endl;
+ std::cout << "material frame rotated around axis: " << axis_string << " by an angle of: " << angle << std::endl;
+ Dune::FieldMatrix<double,6,6> C_rot = rotationMatrixCompliance(axis,angle);
+ (Compliance.leftmultiply(C_rot)).rightmultiply(C_rot.transposed());
+ }
+ printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
// printmatrix(std::cout, C, "C", "--");
//--- return elasticity tensor
- C.invert();
+ Compliance.invert();
// printmatrix(std::cout, C, "C after .invert()", "--");
- return C;
+ return Compliance;
}
//-------------------------------------------------------------------------------
@@ -540,17 +636,41 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
// --- 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
- )
+ // 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
+ // )
+ Dune::FieldMatrix<double,6,6> general_anisotropic(const int& axis,
+ const double& angle,
+ Dune::FieldMatrix<double,6,6>& Compliance
+ )
{
+
+ std::string axis_string;
+ switch (axis)
+ {
+ case 0: axis_string = "x - axis";
+ case 1: axis_string = "y - axis";
+ case 2: axis_string = "z - axis";
+ default: break;
+ }
+
+ /**
+ * @brief (optional) Tramsform compliance matrix to the correct frame.
+ */
+ if(abs(angle) > 1e-12){
+ // std::cout << "material frame rotated around axis: " << axis << " by an angle of: " << angle << std::endl;
+ std::cout << "material frame rotated around axis: " << axis_string << " by an angle of: " << angle << std::endl;
+ Dune::FieldMatrix<double,6,6> C_rot = rotationMatrixCompliance(axis,angle);
+ (Compliance.leftmultiply(C_rot)).rightmultiply(C_rot.transposed());
+ }
+ printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
+
//--- return elasticity tensor (in Voigt notation)
- auto tmp = C;
- tmp.invert();
+ Compliance.invert();
// printmatrix(std::cout, C, "C after .invert()", "--");
- return tmp;
+ return Compliance;
}
//-------------------------------------------------------------------------------
diff --git a/materials/material_orthotropic.py b/materials/material_orthotropic.py
index ffaf3901..3f9ab71e 100644
--- a/materials/material_orthotropic.py
+++ b/materials/material_orthotropic.py
@@ -41,7 +41,19 @@ materialParameters_phase1 = [11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37] # w
#- PHASE 2
phase2_type="orthotropic"
-materialParameters_phase2 = [10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31] # Norway spruce parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+# materialParameters_phase2 = [10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31] # Norway spruce parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+materialParameters_phase2 = [11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37]
+
+# Pass a set of FrameVectors to transform material properties to this Frame
+# phase2_FrameVector1 = [1, 0 ,0]
+# phase2_FrameVector2 = [0, 5, 0]
+# phase2_FrameVector3 = [0, 0, 1]
+
+phase2_axis = 2
+phase2_angle = np.pi/2.0
+# phase2_angle = 2*np.pi/12
+
+
#- PHASE 3
phase3_type="isotropic"
--
GitLab