move m*phi assembly out of for-loop

inputs/cellsolver.parset

@@ -68,8 +68,8 @@ alpha = 2.0
 #--- Lame-Parameters
 mu1=1.0
 
-lambda1=0.0
-#lambda1=1.0
+#lambda1=0.0
+lambda1=1.0
 
 
 # ---volume fraction  (default value = 1.0/4.0)

src/Cell-Problem.cc

@@ -194,12 +194,18 @@ void computeElementStiffnessMatrix(const LocalView& localView,
     for (size_t i=0; i<gradients.size(); i++)
       jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
 
-    for (size_t k=0; k < dimWorld; k++)
-      for (size_t i=0; i < nSf; i++)
+    for (size_t l=0; l< dimWorld; l++)
+      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] = gradients[j][0];                       // Y
+        defGradientV[l][1] = gradients[j][1];                       //X2
+        defGradientV[l][2] = (1.0/gamma)*gradients[j][2];           //X3
         // "phi*phi"-part
-        for (size_t l=0; l< dimWorld; l++)
-          for (size_t j=0; j < nSf; j++ )
+        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);
@@ -207,29 +213,24 @@ void computeElementStiffnessMatrix(const LocalView& localView,
             defGradientU[k][1] = gradients[i][1];                       //X2
             defGradientU[k][2] = (1.0/gamma)*gradients[i][2];           //X3
             // printmatrix(std::cout, defGradientU , "defGradientU", "--");
-            // (scaled)  Deformation gradient of the test basis function
-            MatrixRT defGradientV(0);
-            defGradientV[l][0] = gradients[j][0];                       // Y
-            defGradientV[l][1] = gradients[j][1];                       //X2
-            defGradientV[l][2] = (1.0/gamma)*gradients[j][2];           //X3
-              
+
             double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
 //             double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);  // also works..
 
             size_t col = localView.tree().child(k).localIndex(i);                       
-            size_t row = localView.tree().child(l).localIndex(j);
+            
             elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
+           }
             
-            // "m*phi"  & "phi*m" - part
-            for( size_t m=0; m<3; m++)
-            {
-                double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
-
-                auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
-                elementMatrix[row][localPhiOffset+m] += value;
-                elementMatrix[localPhiOffset+m][row] += value;
+           // "m*phi"  & "phi*m" - part
+           for( size_t m=0; m<3; m++)
+           {
+               double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+               auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+               elementMatrix[row][localPhiOffset+m] += value;
+               elementMatrix[localPhiOffset+m][row] += value;
             }
-          }
+          
       }
     // "m*m"-part
     for(size_t m=0; m<3; m++)
@@ -418,10 +419,13 @@ void assembleCellStiffness(const Basis& basis,
     muLocal.bind(element);
     lambdaLocal.bind(element);
     const int localPhiOffset = localView.size();
-    
+//     Dune::Timer Time;
     //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
     Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
     computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
+    
+//     std::cout << "local assembly time:" << Time.elapsed() << std::endl;
+    
     //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
     //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
     //std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
@@ -848,6 +852,8 @@ double computeL2Norm(const Basis& basis,
 int main(int argc, char *argv[])
 {
   MPIHelper::instance(argc, argv);
+  
+  Dune::Timer globalTimer;
 
   ParameterTree parameterSet;
   if (argc < 2)
@@ -1875,4 +1881,6 @@ int main(int argc, char *argv[])
     log       << std::string(tableWidth*7 + 3*7, '-') << "\n";  
 
     log.close();
+    
+//     std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
 }