From 696882a7e6bac19d577232fe0a791fa6529f3304 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 10 Jul 2023 15:16:23 +0200
Subject: [PATCH] Precompute the deformation gradients for each shape function
& component
This saves about 20% run time.
---
dune/microstructure/CorrectorComputer.hh | 36 +++++++++++-------------
1 file changed, 16 insertions(+), 20 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 3c78ef0a..8a6f9c93 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -303,6 +303,20 @@ public:
for (size_t i=0; i< jacobians.size(); i++)
jacobians[i] = jacobians[i] * geometryJacobianInverse;
+ // Compute the scaled deformation gradient for each shape function
+ std::vector<std::array<MatrixRT,dimworld> > deformationGradient(nSf);
+
+ for (size_t i=0; i < nSf; i++)
+ {
+ for (size_t k=0; k< dimworld; k++)
+ {
+ MatrixRT defGradientV(0);
+ defGradientV[k] = jacobians[i][0];
+
+ deformationGradient[i][k] = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
+ }
+ }
+
// Compute the material phase that the current quadrature point is in
const auto phase = localIndicatorFunction(quadPos);
@@ -310,29 +324,11 @@ public:
for (size_t j=0; j < nSf; j++ )
{
size_t row = localView.tree().child(l).localIndex(j);
- // (scaled) Deformation gradient of the test basis function
- MatrixRT defGradientV(0);
- defGradientV[l][0] = jacobians[j][0][0]; // Y
- defGradientV[l][1] = jacobians[j][0][1]; //X2
- // defGradientV[l][2] = (1.0/gamma)*jacobians[j][0][2]; //X3
- defGradientV[l][2] = jacobians[j][0][2]; //X3
- defGradientV = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
// "phi*phi"-part
for (size_t k=0; k < dimworld; k++)
for (size_t i=0; i < nSf; i++)
{
- // (scaled) Deformation gradient of the ansatz basis function
- MatrixRT defGradientU(0);
- defGradientU[k][0] = jacobians[i][0][0]; // Y
- defGradientU[k][1] = jacobians[i][0][1]; //X2
- // defGradientU[k][2] = (1.0/gamma)*jacobians[i][0][2]; //X3
- defGradientU[k][2] = jacobians[i][0][2]; //X3
- // printmatrix(std::cout, defGradientU , "defGradientU", "--");
- defGradientU = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientU));
-
-
-
// auto etmp = material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos));
// auto etmp = elasticityTensor(defGradientU,element.geometry().global(quadPos));
// auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos);
@@ -375,7 +371,7 @@ public:
// printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
- double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,phase),sym(defGradientV));
+ double energyDensity= scalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),sym(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));
@@ -416,7 +412,7 @@ public:
- double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(defGradientV));
+ double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(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;
--
GitLab