diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 8d8c53e716cfd48de33612313846ee798ae1a342..3309954055f07413fdbf047825fca443e7903f14 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/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;