From: Denis Davydov <davydden@gmail.com>
Date: Fri, 29 May 2015 18:59:10 +0000 (+0200)
Subject: exploit symmetries in Local integrators of mass and laplace matrices without assignme... 
X-Git-Tag: v8.3.0-rc1~139^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b6a2ff90074e59bb6bdf4d9873395ef72a6a5994;p=dealii.git

exploit symmetries in Local integrators of mass and laplace matrices without assignment operation
---

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);
-
     }
 
     /**