Add method return the symmetric part of a matrix in Voigt notation

That's cheaper than first computing the symmetric part as a matrix,
and then converting that to Voigt notation.

dune/microstructure/CorrectorComputer.hh

@@ -310,7 +310,7 @@ public:
           MatrixRT defGradientV(0);
           defGradientV[k] = jacobians[i][0];
 
-          deformationGradient[i][k] = matrixToVoigt(sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV)));
+          deformationGradient[i][k] = symVoigt(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
         }
       }
 
@@ -581,7 +581,7 @@ public:
   //         tmpDefGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2];                     //
           tmpDefGradientV[k][2] = jacobians[i][0][2];                     // X3
           
-          VoigtVector<double,3> defGradientV = matrixToVoigt(sym(crossSectionDirectionScaling((1.0/gamma_),tmpDefGradientV)));
+          VoigtVector<double,3> defGradientV = symVoigt(crossSectionDirectionScaling((1.0/gamma_),tmpDefGradientV));
 
           // double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
 
@@ -1373,7 +1373,7 @@ public:
           
           auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + mContainer[b];
 
-          const auto strain1 = matrixToVoigt(sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi1(quadPos))));
+          const auto strain1 = symVoigt(crossSectionDirectionScaling(1.0/gamma_, DerPhi1(quadPos)));
           
         
           auto G = matrixFieldG(quadPos);
@@ -1381,7 +1381,7 @@ public:
           
           auto tmp = matrixToVoigt(G + mContainer[a]) + strain1;
 
-          double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),matrixToVoigt(sym(Chi)));
+          double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),symVoigt(Chi));
           // double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
           // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
 

dune/microstructure/voigthelper.hh

@@ -36,4 +36,19 @@ 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];
 }
 
+/** \brief Return the symmetric part of a matrix, in Voigt notation */
+template<typename T>
+VoigtVector<T,3> symVoigt(const Dune::FieldMatrix<T,3,3>& matrix)
+{
+  VoigtVector<T,3> result = {matrix[0][0],
+                             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])};
+
+  return result;
+}
+
 #endif   // DUNE_MICROSTRUCTURE_VOIGTHELPER_HH