From 9be7b2cda6810b0ddac662774eec03ce66dc5122 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 10 Jul 2023 14:34:52 +0200
Subject: [PATCH] Compute phase only once per quadrature point
Computing the material phase at a particular point may involve
a call to a Python function. That is usually quite costly.
Therefore, make sure that this call is really made only once
per quadrature point.
---
dune/microstructure/CorrectorComputer.hh | 11 +++++++----
1 file changed, 7 insertions(+), 4 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index cdc6dd1a..f0b97b97 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -309,6 +309,9 @@ public:
for (size_t i=0; i<gradients.size(); i++)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+ // Compute the material phase that the current quadrature point is in
+ const auto phase = localIndicatorFunction(quadPos);
+
for (size_t l=0; l< dimworld; l++)
for (size_t j=0; j < nSf; j++ )
{
@@ -378,7 +381,7 @@ public:
// printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
- double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
+ double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,phase),sym(defGradientV));
// 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));
@@ -419,7 +422,7 @@ public:
- double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
+ double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(defGradientV));
// 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;
@@ -468,7 +471,7 @@ public:
- double energyDensityGG = scalarProduct(material_.applyElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
+ double energyDensityGG = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),sym(basisContainer[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]));
@@ -1473,4 +1476,4 @@ public:
}; // end class
-#endif
\ No newline at end of file
+#endif
--
GitLab