From 897f9939e4c84203d64d19f5964e946f712f5b3b Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 2 Oct 2023 12:33:57 +0200
Subject: [PATCH] Compute the phases at the quad points outside of the element
 assembler

This will be the cache key.
---
 dune/microstructure/CorrectorComputer.hh | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index c1172775..eebca42d 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -244,9 +244,11 @@ public:
   /** \brief Compute the stiffness matrix for one element
    *
    * \param quadRule The quadrature rule to be used
+   * \param phaseAtQuadPoint The material phase (i.e., the type of material) at each quadrature point
    */
   void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
                                      const Dune::QuadratureRule<double,dim>& quadRule,
+                                     const std::vector<int>& phaseAtQuadPoint,
                                     ElementMatrixCT& elementMatrix
                                     )
   {
@@ -262,14 +264,6 @@ public:
 
 
     // auto elasticityTensor = material_.getElasticityTensor();
-    // auto indicatorFunction = material_.getIndicatorFunction();
-    auto localIndicatorFunction = material_.getLocalIndicatorFunction();
-    localIndicatorFunction.bind(element);
-    // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
-    // auto indicatorFunction = localFunction(indicatorFunctionGVF);
-    // indicatorFunction.bind(element);
-    // Func2int indicatorFunction = *material_.getIndicatorFunction();
-
 
     // LocalBasis-Offset
     const int localPhiOffset = localView.size();
@@ -308,7 +302,7 @@ public:
       }
 
       // Compute the material phase that the current quadrature point is in
-      const auto phase = localIndicatorFunction(quadPos);
+      const auto phase = phaseAtQuadPoint[QPcounter-1];
 
       for (size_t l=0; l< dimworld; l++)
       for (size_t j=0; j < nSf; j++ )
@@ -658,9 +652,18 @@ public:
       int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
       const auto& quadRule = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
 
+      // Determine the material phase at each quadrature point
+      auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+      localIndicatorFunction.bind(element);
+
+      std::vector<int> phaseAtQuadPoint(quadRule.size());
+
+      for (std::size_t i=0; i<quadRule.size(); ++i)
+        phaseAtQuadPoint[i] = localIndicatorFunction(quadRule[i].position());
+
       ElementMatrixCT elementMatrix;
       // computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
-      computeElementStiffnessMatrix(localView, quadRule, elementMatrix);
+      computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
       
       
   //     std::cout << "local assembly time:" << Time.elapsed() << std::endl;
-- 
GitLab