Fix Voigt scaling issue

Voigt-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 gets scaled by a factor of 2 in the last three columns
  - The 'voigtScalarProduct' scales the last three products by a factor of 2

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);

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;
   }
   //-------------------------------------------------------------------------------

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;
 }