exploit symmetries in Local integrators of mass and laplace matrices without assignme...

diff --git a/include/deal.II/integrators/l2.h b/include/deal.II/integrators/l2.h
index baadc9cbabb2c4bac07c7ba7bcfcad5060753af0..41672153ee28e10989e5c4ab84d736410e03aefa 100644 (file)
--- 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 bf294d6e3a2e02f10e2d0d93dc0d98fbc0269f62..8c74227f87a6462756f1bfcd6763d8de86bf82b5 100644 (file)
--- 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);
-
     }
 
     /**