Store the basis of matrix space in Voigt notation

This again saves a few conversions.

dune/microstructure/CorrectorComputer.hh

@@ -67,11 +67,9 @@ protected:
   const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
 
   // ---- Basis for R_sym^{2x2}
-  constexpr static double sqrtOf2 = 1.414213562373095;  // The sqrt method is not constexpr
-  constexpr static MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
-  constexpr static MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
-  constexpr static MatrixRT G3_ {{0.0, 1.0/sqrtOf2, 0.0}, {1.0/sqrtOf2, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-  constexpr static std::array<MatrixRT,3 > matrixBasis_ = {G1_, G2_, G3_};
+  // Could really be constexpr static, but then we would need to write the basis here directly in Voigt notation.
+  // However, I suspect that our Voigt representation may still change in the future.
+  const std::array<VoigtVector<double,3>,3 > matrixBasis_;
 
   Func2Tensor x3G_1_ = [] (const Domain& x) {
                             return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};    //TODO könnte hier sign übergeben?
@@ -104,6 +102,9 @@ public:
             gamma_(gamma),
             log_(log),
             parameterSet_(parameterSet),
+            matrixBasis_(std::array<VoigtVector<double,3>,3>{matrixToVoigt(Dune::FieldMatrix<double,3,3>({{1, 0, 0}, {0, 0, 0}, {0, 0, 0}})),
+                                                             matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
+                                                             matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))}),
             phiOffset_(basis.size())
     {
 
@@ -157,7 +158,7 @@ public:
 
     
     // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
-    auto getMatrixBasiscontainer(){return make_shared<std::array<MatrixRT,3 >>(matrixBasis_);}
+    auto getMatrixBasiscontainer(){return make_shared<std::array<VoigtVector<double,3>,3 >>(matrixBasis_);}
     // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
     auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
 
@@ -408,7 +409,7 @@ public:
 
 
 
-              double energyDensityGphi = voigtScalarProduct(material_.applyElasticityTensor(matrixToVoigt(matrixBasis_[m]),phase),deformationGradient[j][l]);
+              double energyDensityGphi = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),deformationGradient[j][l]);
               // double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
               // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
               // 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;
@@ -457,7 +458,7 @@ public:
 
 
 
-          double energyDensityGG = voigtScalarProduct(material_.applyElasticityTensor(matrixToVoigt(matrixBasis_[m]),phase),matrixToVoigt(matrixBasis_[n]));
+          double energyDensityGG = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),matrixBasis_[n]);
           // double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
           // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
           // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
@@ -531,14 +532,6 @@ public:
     // LocalBasis-Offset
     const int localPhiOffset = localView.size();
 
-    ///////////////////////////////////////////////
-    // Basis for R_sym^{2x2}
-    //////////////////////////////////////////////
-    MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-    MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}};
-    MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
-    std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-
   //   int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;  // TEST
   //   int orderQR = 0;
   //   int orderQR = 1;
@@ -620,7 +613,7 @@ public:
         // double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
 
 
-        double energyDensityfG= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixToVoigt(basisContainer[m]));
+        double energyDensityfG= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixBasis_[m]);
         // double energyDensityfG= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
         // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
         // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
@@ -1101,9 +1094,9 @@ public:
 
       for(size_t i=0; i<3; i++)
       {
-          mContainer[0] += m_1_[i]*matrixBasis_[i];
-          mContainer[1] += m_2_[i]*matrixBasis_[i];
-          mContainer[2] += m_3_[i]*matrixBasis_[i];
+          mContainer[0] += m_1_[i]*voigtToMatrix(matrixBasis_[i]);
+          mContainer[1] += m_2_[i]*voigtToMatrix(matrixBasis_[i]);
+          mContainer[2] += m_3_[i]*voigtToMatrix(matrixBasis_[i]);
       }
 
       std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;