From: Denis Davydov <davydden@gmail.com>
Date: Fri, 29 May 2015 20:24:02 +0000 (+0200)
Subject: accumulate sums over n_components
X-Git-Tag: v8.3.0-rc1~139^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F978%2Fhead;p=dealii.git

accumulate sums over n_components
---

diff --git a/include/deal.II/integrators/l2.h b/include/deal.II/integrators/l2.h
index 41672153ee..82071cc40a 100644
--- 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 8c74227f87..a802e13a16 100644
--- 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;
+                }
+            }
         }
 
     }