diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index d81552408b56790f2e5a3ae0a1584916f22360b8..8d8c53e716cfd48de33612313846ee798ae1a342 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -12,6 +12,8 @@
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
#include <dune/solvers/solvers/umfpacksolver.hh>
+#include <dune/microstructure/voigthelper.hh>
+
using namespace MatrixOperations;
using std::shared_ptr;
@@ -298,7 +300,7 @@ public:
jacobians[i] = jacobians[i] * geometryJacobianInverse;
// Compute the scaled deformation gradient for each shape function
- std::vector<std::array<MatrixRT,dimworld> > deformationGradient(nSf);
+ std::vector<std::array<VoigtVector<double,dim>,dimworld> > deformationGradient(nSf);
for (size_t i=0; i < nSf; i++)
{
@@ -307,7 +309,7 @@ public:
MatrixRT defGradientV(0);
defGradientV[k] = jacobians[i][0];
- deformationGradient[i][k] = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
+ deformationGradient[i][k] = matrixToVoigt(sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV)));
}
}
@@ -365,7 +367,7 @@ public:
// printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
- double energyDensity= scalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
+ double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
// double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
@@ -406,7 +408,7 @@ public:
- double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),deformationGradient[j][l]);
+ double energyDensityGphi = voigtScalarProduct(material_.applyElasticityTensor(matrixToVoigt(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;
@@ -455,7 +457,7 @@ public:
- double energyDensityGG = scalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),matrixBasis_[n]);
+ double energyDensityGG = voigtScalarProduct(material_.applyElasticityTensor(matrixToVoigt(matrixBasis_[m]),phase),matrixToVoigt(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]));
@@ -580,19 +582,19 @@ public:
for (size_t k=0; k < dimworld; k++)
{
// Deformation gradient of the ansatz basis function
- MatrixRT defGradientV(0);
- defGradientV[k][0] = jacobians[i][0][0]; // Y
- defGradientV[k][1] = jacobians[i][0][1]; // X2
- // defGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2]; //
- defGradientV[k][2] = jacobians[i][0][2]; // X3
+ MatrixRT tmpDefGradientV(0);
+ tmpDefGradientV[k][0] = jacobians[i][0][0]; // Y
+ tmpDefGradientV[k][1] = jacobians[i][0][1]; // X2
+ // tmpDefGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2]; //
+ tmpDefGradientV[k][2] = jacobians[i][0][2]; // X3
- defGradientV = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
+ VoigtVector<double,3> defGradientV = matrixToVoigt(sym(crossSectionDirectionScaling((1.0/gamma_),tmpDefGradientV)));
// double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
- double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
+ double energyDensity= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),defGradientV);
// double energyDensity= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
@@ -618,7 +620,7 @@ public:
// double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
- double energyDensityfG= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
+ double energyDensityfG= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixToVoigt(basisContainer[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]));
@@ -1378,19 +1380,15 @@ public:
auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + mContainer[b];
- auto strain1 = DerPhi1(quadPos);
-
-
- strain1 = crossSectionDirectionScaling(1.0/gamma_, strain1);
- strain1 = sym(strain1);
+ const auto strain1 = matrixToVoigt(sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi1(quadPos))));
auto G = matrixFieldG(quadPos);
// auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
- auto tmp = G + mContainer[a] + strain1;
+ auto tmp = matrixToVoigt(G + mContainer[a]) + strain1;
- double energyDensity= scalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
+ double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),matrixToVoigt(sym(Chi)));
// double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index b6f3e8aff123858541694bd5eaa089bedcc79ac7..c6e358d4c796396fba5d84d5371d5a975eda3457 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -170,10 +170,10 @@ public:
// auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
// auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
- auto X1 = G1 + Chi1;
+ auto X1 = matrixToVoigt(G1 + Chi1);
// auto X2 = G2 + Chi2;
- double energyDensity= scalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
+ double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)), matrixToVoigt(sym(G2)));
// double energyDensity= scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
@@ -188,7 +188,7 @@ public:
// printmatrix(std::cout, prestrainPhaseValue_old, "prestrainPhaseValue_old", "--");
// printmatrix(std::cout, prestrainPhaseValue, "prestrainPhaseValue", "--");
// }
- auto value = scalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue));
+ auto value = voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),matrixToVoigt(sym(prestrainPhaseValue)));
elementPrestrain += value * quadPoint.weight() * integrationElement;
// elementPrestrain += scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue)) * quadPoint.weight() * integrationElement;
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 16e43bff33a63940e9cd500300c8405e5410a7df..4574cafab436b2cea03ee4e98b46ff13d9e5672d 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -596,17 +596,15 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
//-------------------------------------------------------------------------------
- MatrixRT applyElasticityTensor(const MatrixRT& G, const int& phase) const
+ VoigtVector<double,3> applyElasticityTensor(const VoigtVector<double,3>& G, const int& phase) const
{
-
- VoigtVector<double,3> G_tmp = strainToVoigt(G);
VoigtVector<double,3> out(0);
// printmatrix(std::cout, L_[0], "L_[0]", "--");
- L_[phase-1].mv(G_tmp,out);
+ L_[phase-1].mv(G,out);
- return voigtToStress(out);
- }
+ return out;
+ }
////////////////////////////////////////////////////////////////////////////////////////
diff --git a/dune/microstructure/voigthelper.hh b/dune/microstructure/voigthelper.hh
index cf5063a64c151cc9a9a2ee24f28ca585a10b0cdd..24984e84f870ad001924a89421198dbeaba1558e 100644
--- a/dune/microstructure/voigthelper.hh
+++ b/dune/microstructure/voigthelper.hh
@@ -28,28 +28,12 @@ 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]}};
}
+/** \brief Interpret Voigt vectors as symmetric matrices and compute their Frobenius scalar voigtScalarProduct
+ */
template<typename T>
-VoigtVector<T,3> strainToVoigt(const Dune::FieldMatrix<T,3,3>& G)
-{
- return {G[0][0], G[1][1], G[2][2], 2.0*G[1][2], 2.0*G[0][2], 2.0*G[0][1]};
-}
-
-template<typename T>
-VoigtVector<T,3> stressToVoigt(const Dune::FieldMatrix<T,3,3>& G)
-{
- return {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
-}
-
-template<typename T>
-Dune::FieldMatrix<T,3,3> voigtToStress(const VoigtVector<T,3>& v)
-{
- return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
-}
-
-template<typename T>
-Dune::FieldMatrix<T,3,3> voigtToStrain(const Dune::FieldVector<T,6>& v)
+T voigtScalarProduct(const VoigtVector<T,3>& a, const VoigtVector<T,3>& b)
{
- 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 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];
}
#endif // DUNE_MICROSTRUCTURE_VOIGTHELPER_HH