From b6a2ff90074e59bb6bdf4d9873395ef72a6a5994 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Fri, 29 May 2015 20:59:10 +0200
Subject: [PATCH] exploit symmetries in Local integrators of mass and laplace
 matrices without assignment operation

---
 include/deal.II/integrators/l2.h      | 45 ++++++++++++++++-----------
 include/deal.II/integrators/laplace.h | 24 +++++++-------
 2 files changed, 40 insertions(+), 29 deletions(-)

diff --git a/include/deal.II/integrators/l2.h b/include/deal.II/integrators/l2.h
index baadc9cbab..41672153ee 100644
--- a/include/deal.II/integrators/l2.h
+++ b/include/deal.II/integrators/l2.h
@@ -63,17 +63,21 @@ namespace LocalIntegrators
           const double dx = fe.JxW(k) * factor;
 
           for (unsigned int i=0; i<n_dofs; ++i)
-            for (unsigned int j=i; j<n_dofs; ++j)
-              for (unsigned int d=0; d<n_components; ++d)
-                M(i,j) += dx
-                          * fe.shape_value_component(j,k,d)
+            for (unsigned int d=0; d<n_components; ++d)
+              {
+                M(i,i) += dx
+                          * fe.shape_value_component(i,k,d)
                           * fe.shape_value_component(i,k,d);
+                for (unsigned int j=i+1; j<n_dofs; ++j)
+                  {
+                    const double Mij =  dx
+                                        * fe.shape_value_component(j,k,d)
+                                        * fe.shape_value_component(i,k,d);
+                    M(i,j) += Mij;
+                    M(j,i) += Mij;
+                  }
+              }
         }
-
-      // exploit symmetry
-      for (unsigned int i=0; i<n_dofs; ++i)
-        for (unsigned int j=i+1; j<n_dofs; ++j)
-          M(j,i) = M(i,j);
     }
 
     /**
@@ -108,17 +112,22 @@ namespace LocalIntegrators
           const double dx = fe.JxW(k) * weights[k];
 
           for (unsigned int i=0; i<n_dofs; ++i)
-            for (unsigned int j=i; j<n_dofs; ++j)
-              for (unsigned int d=0; d<n_components; ++d)
-                M(i,j) += dx
-                          * fe.shape_value_component(j,k,d)
+            for (unsigned int d=0; d<n_components; ++d)
+              {
+                M(i,i) += dx
+                          * fe.shape_value_component(i,k,d)
                           * fe.shape_value_component(i,k,d);
-        }
 
-      //exploit symmetry
-      for (unsigned int i=0; i<n_dofs; ++i)
-        for (unsigned int j=i+1; j<n_dofs; ++j)
-          M(j,i) = M(i,j);
+                for (unsigned int j=i+1; j<n_dofs; ++j)
+                  {
+                    const double Mij = dx
+                                       * fe.shape_value_component(j,k,d)
+                                       * fe.shape_value_component(i,k,d);
+                    M(i,j) += Mij;
+                    M(j,i) += Mij;
+                  }
+              }
+        }
     }
 
     /**
diff --git a/include/deal.II/integrators/laplace.h b/include/deal.II/integrators/laplace.h
index bf294d6e3a..8c74227f87 100644
--- a/include/deal.II/integrators/laplace.h
+++ b/include/deal.II/integrators/laplace.h
@@ -61,19 +61,21 @@ namespace LocalIntegrators
         {
           const double dx = fe.JxW(k) * factor;
           for (unsigned int i=0; i<n_dofs; ++i)
-            {
-              for (unsigned int j=i; j<n_dofs; ++j)
-                for (unsigned int d=0; d<n_components; ++d)
-                  M(i,j) += dx *
-                            (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
-            }
+            for (unsigned int d=0; d<n_components; ++d)
+              {
+                M(i,i) += dx *
+                          (fe.shape_grad_component(i,k,d) * fe.shape_grad_component(i,k,d));
+
+                for (unsigned int j=i+1; j<n_dofs; ++j)
+                  {
+                    const double Mij = dx *
+                                       (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
+                    M(i,j) += Mij;
+                    M(j,i) += Mij;
+                  }
+              }
         }
 
-      // exploit symmetry
-      for (unsigned int i=0; i<n_dofs; ++i)
-        for (unsigned int j=i+1; j<n_dofs; ++j)
-          M(j,i) = M(i,j);
-
     }
 
     /**
-- 
2.39.5