From c4ac0bd63e46bb8eced07c28932c1f773ee40655 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 2 Oct 2023 11:21:46 +0200
Subject: [PATCH] Store deformation gradient and stress as Voigt vectors
This avoids many transformations from symmetric matrices to
Voigt vectors and back. In my (limited) testing, this reduces
the time to assemble the stiffness matrix by about 25%.
This patch also introduces a custom scalar product method for
Voigt vectors, which reproduces the Frobenius scalar product
in matrix space. That way, the potentially confusing distinction
between stress-like Voigt vectors and strain-like Voigt vectors
can be avoided.
---
dune/microstructure/CorrectorComputer.hh | 38 +++++++++----------
.../EffectiveQuantitiesComputer.hh | 6 +--
dune/microstructure/prestrainedMaterial.hh | 10 ++---
dune/microstructure/voigthelper.hh | 24 ++----------
4 files changed, 29 insertions(+), 49 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index d8155240..8d8c53e7 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 b6f3e8af..c6e358d4 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 16e43bff..4574cafa 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 cf5063a6..24984e84 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
--
GitLab