accumulate sums over n_components

diff --git a/include/deal.II/integrators/l2.h b/include/deal.II/integrators/l2.h
index 41672153ee28e10989e5c4ab84d736410e03aefa..82071cc40a6e65bfc8458237859c1e4164e1d642 100644 (file)
--- a/include/deal.II/integrators/l2.h
+++ b/include/deal.II/integrators/l2.h
@@ -61,22 +61,28 @@ namespace LocalIntegrators
       for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
         {
           const double dx = fe.JxW(k) * factor;
-
           for (unsigned int i=0; i<n_dofs; ++i)
-            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;
-                  }
-              }
+            {
+              double Mii = 0.0;
+              for (unsigned int d=0; d<n_components; ++d)
+                Mii += dx
+                       * fe.shape_value_component(i,k,d)
+                       * fe.shape_value_component(i,k,d);
+
+              M(i,i) += Mii;
+
+              for (unsigned int j=i+1; j<n_dofs; ++j)
+                {
+                  double Mij = 0.0;
+                  for (unsigned int d=0; d<n_components; ++d)
+                    Mij += dx
+                           * fe.shape_value_component(j,k,d)
+                           * fe.shape_value_component(i,k,d);
+
+                  M(i,j) += Mij;
+                  M(j,i) += Mij;
+                }
+            }
         }
     }
 
@@ -110,23 +116,28 @@ namespace LocalIntegrators
       for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
         {
           const double dx = fe.JxW(k) * weights[k];
-
           for (unsigned int i=0; i<n_dofs; ++i)
-            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;
-                  }
-              }
+            {
+              double Mii = 0.0;
+              for (unsigned int d=0; d<n_components; ++d)
+                Mii += dx
+                       * fe.shape_value_component(i,k,d)
+                       * fe.shape_value_component(i,k,d);
+
+              M(i,i) += Mii;
+
+              for (unsigned int j=i+1; j<n_dofs; ++j)
+                {
+                  double Mij = 0.0;
+                  for (unsigned int d=0; d<n_components; ++d)
+                    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 8c74227f87a6462756f1bfcd6763d8de86bf82b5..a802e13a1697235c13977d96a195eeff2a2ed7c4 100644 (file)
--- a/include/deal.II/integrators/laplace.h
+++ b/include/deal.II/integrators/laplace.h
@@ -61,19 +61,25 @@ namespace LocalIntegrators
         {
           const double dx = fe.JxW(k) * factor;
           for (unsigned int i=0; i<n_dofs; ++i)
-            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;
-                  }
-              }
+            {
+              double Mii = 0.0;
+              for (unsigned int d=0; d<n_components; ++d)
+                Mii += dx *
+                       (fe.shape_grad_component(i,k,d) * fe.shape_grad_component(i,k,d));
+
+              M(i,i) += Mii;
+
+              for (unsigned int j=i+1; j<n_dofs; ++j)
+                {
+                  double Mij = 0.0;
+                  for (unsigned int d=0; d<n_components; ++d)
+                    Mij += dx *
+                           (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
+
+                  M(i,j) += Mij;
+                  M(j,i) += Mij;
+                }
+            }
         }
 
     }