diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 1796b1b2ffcf1686ba58b552c376b7d38331f20f..59224003ec84df61f6de1c1c19be766679d62981 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -143,6 +143,22 @@ namespace MatrixOperations {
};
}
+
+ /**
+ * @brief Scale last three columns of stiffness matrix by a factor of 2.
+ * Inserting this facotr in the stiffness matrix allows to use the
+ * same Matrix-to-Vector mapping for both the stress and the strain.
+ *
+ * @param S Stiffness matrix
+ */
+ static void scaleStiffnessMatrix(Dune::FieldMatrix<double,6,6>& S){
+ for(size_t i = 0; i<6; i++)
+ for(size_t j = 3; j<6; j++)
+ {
+ S[i][j] = 2.0*S[i][j];
+ }
+ }
+
// 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 4574cafab436b2cea03ee4e98b46ff13d9e5672d..deae82234e1553f7c3c0646f9e11b854e3b4fe01 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -395,15 +395,15 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
// MATERIAL CLASSES DEFINITIONS:
////////////////////////////////////////////////////////
- //--- Definition of isotropic elasticity tensor (stiffness matrix in voigt notation)
+ /** \brief Isotropic stiffness matrix in Voigt notation but scaled by a factor of 2 in the last three columns. */
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 , 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 }
};
}
//-------------------------------------------------------------------------------
@@ -476,12 +476,14 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
(Compliance.leftmultiply(C_rot)).rightmultiply(C_rot.transposed());
}
-
printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
- // invert to obtain elasticity/stiffness-Tensor
+ // invert Compliance matrix to obtain elasticity/stiffness matrix
Compliance.invert();
// printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
+ scaleStiffnessMatrix(Compliance);
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor scaled", "--");
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
return Compliance;
}
//-------------------------------------------------------------------------------
@@ -548,7 +550,10 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
//--- return elasticity tensor
Compliance.invert();
- // printmatrix(std::cout, C, "C after .invert()", "--");
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
+
+ scaleStiffnessMatrix(Compliance);
+
return Compliance;
}
//-------------------------------------------------------------------------------
@@ -588,9 +593,15 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
}
printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
- //--- return elasticity tensor (in Voigt notation)
+
+
+
+ //--- return stiffness/elasticity tensor (in Voigt notation)
Compliance.invert();
- // printmatrix(std::cout, C, "C after .invert()", "--");
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
+
+ scaleStiffnessMatrix(Compliance);
+
return Compliance;
}
//-------------------------------------------------------------------------------
diff --git a/dune/microstructure/voigthelper.hh b/dune/microstructure/voigthelper.hh
index 05668a29891b7106b7420242e168964c753e1ec4..5058f710ec2b21a535b8da7f51aec7082d498646 100644
--- a/dune/microstructure/voigthelper.hh
+++ b/dune/microstructure/voigthelper.hh
@@ -16,16 +16,26 @@ template<typename T, int dim>
using VoigtMatrix = Dune::FieldMatrix<T,dim*(dim+1)/2,dim*(dim+1)/2>;
+/**
+ * @brief Voigt-notation(https://de.wikipedia.org/wiki/Voigtsche_Notation) distinguishes
+ * in the transformation from Matrix to Vector between stresses and strains. The transformation
+ * for strains features an additional factor 2 for the non-diagonal entries. In order to avoid
+ * the use of different data structures for both stresses & strains we use the same Matrix-to-Vector
+ * mapping ('matrixToVoigt') and incorporate the factors in suitable places. namely:
+ * - The Stiffness matrix of the constitutive relation get scaled by a factor of 2 in the last three columns
+ * - The 'voigtScalarProduct' scales the lase three products by a factor of 2
+ */
+
template<typename T>
VoigtVector<T,3> matrixToVoigt(const Dune::FieldMatrix<T,3,3>& matrix)
{
- return {matrix[0][0], matrix[1][1], matrix[2][2], 2.0*matrix[1][2], 2.0*matrix[0][2], 2.0*matrix[0][1]};
+ return {matrix[0][0], matrix[1][1], matrix[2][2], matrix[1][2], matrix[0][2], matrix[0][1]};
}
template<typename T>
Dune::FieldMatrix<T,3,3> voigtToMatrix(const VoigtVector<T,3>& v)
{
- 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]}};
+ return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
}
/** \brief Interpret Voigt vectors as symmetric matrices and compute their Frobenius scalar voigtScalarProduct
@@ -33,7 +43,7 @@ Dune::FieldMatrix<T,3,3> voigtToMatrix(const VoigtVector<T,3>& v)
template<typename T>
T voigtScalarProduct(const VoigtVector<T,3>& a, const VoigtVector<T,3>& b)
{
- return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + 2*a[3]*b[3] + 2*a[4]*b[4] + 2*a[5]*b[5];
+ return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + 2.0*a[3]*b[3] + 2.0*a[4]*b[4] + 2.0*a[5]*b[5];
}
/** \brief Return the symmetric part of a matrix, in Voigt notation */
@@ -44,9 +54,9 @@ VoigtVector<T,3> symVoigt(const Dune::FieldMatrix<T,3,3>& matrix)
matrix[1][1],
matrix[2][2],
// The 0.5 is part of the sym operator, the 2.0 is from the conversion to Voigt notation
- 2 * 0.5 * (matrix[1][2] + matrix[2][1]),
- 2 * 0.5 * (matrix[0][2] + matrix[2][0]),
- 2 * 0.5 * (matrix[0][1] + matrix[1][0])};
+ 0.5 * (matrix[1][2] + matrix[2][1]),
+ 0.5 * (matrix[0][2] + matrix[2][0]),
+ 0.5 * (matrix[0][1] + matrix[1][0])};
return result;
}