From 8642a2e42a0e475df0e6749632782fc87a248b39 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 29 Sep 2023 15:09:52 +0200
Subject: [PATCH] Define basis of matrix space only once

---
 dune/microstructure/CorrectorComputer.hh | 21 +++++++--------------
 1 file changed, 7 insertions(+), 14 deletions(-)

diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 576d2afa..b64f5a5e 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -64,10 +64,11 @@ protected:
   const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
 
   // ---- Basis for R_sym^{2x2}
-  MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
-  MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
-  MatrixRT G3_ {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
-  std::array<MatrixRT,3 > MatrixBasisContainer_ = {G1_, G2_, G3_};
+  constexpr static double sqrtOf2 = 1.414213562373095;  // The sqrt method is not constexpr
+  constexpr static MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+  constexpr static MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+  constexpr static MatrixRT G3_ {{0.0, 1.0/sqrtOf2, 0.0}, {1.0/sqrtOf2, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+  constexpr static std::array<MatrixRT,3 > MatrixBasisContainer_ = {G1_, G2_, G3_};
 
   Func2Tensor x3G_1_ = [] (const Domain& x) {
                             return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};    //TODO könnte hier sign übergeben?
@@ -269,14 +270,6 @@ public:
   //   std::cout << "localView.size(): " << localView.size() << std::endl;
   //   std::cout << "nSf : " << nSf << std::endl;
 
-    ///////////////////////////////////////////////
-    // Basis for R_sym^{2x2}            // wird nicht als Funktion benötigt da konstant...
-    //////////////////////////////////////////////
-    MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-    MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}};
-    MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
-    std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-
   //   int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;  // TEST
     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
   //   int orderQR = 0;
@@ -412,7 +405,7 @@ public:
 
 
 
-              double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),deformationGradient[j][l]);
+              double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(MatrixBasisContainer_[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;
@@ -461,7 +454,7 @@ public:
 
 
 
-          double energyDensityGG = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),basisContainer[n]);
+          double energyDensityGG = scalarProduct(material_.applyElasticityTensor(MatrixBasisContainer_[m],phase),MatrixBasisContainer_[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]));
-- 
GitLab