Remove quartic complexity from create_coupling_mass_matrix

diff --git a/source/non_matching/coupling.cc b/source/non_matching/coupling.cc
index c71f3a57b1e0cd59c84886a26c2a96f62cd6af99..93bb108b8c17e5e16e0a1a589b66fb8e9595ae19 100644 (file)
--- a/source/non_matching/coupling.cc
+++ b/source/non_matching/coupling.cc
@@ -851,27 +851,29 @@ namespace NonMatching
         {
           kernel.set_center(fev0.quadrature_point(q0));
           kernel.value_list(fev1.get_quadrature_points(), kernel_values);
-          for (unsigned int q1 = 0; q1 < quadrature1.size(); ++q1)
-            if (kernel_values[q1] != 0.0)
-              {
-                for (unsigned int i = 0; i < fe0.dofs_per_cell; ++i)
-                  {
-                    const auto comp_i = fe0.system_to_component_index(i).first;
-                    if (gtl0[comp_i] != numbers::invalid_unsigned_int)
-                      for (unsigned int j = 0; j < fe1.dofs_per_cell; ++j)
-                        {
-                          const auto comp_j =
-                            fe1.system_to_component_index(j).first;
-                          if (gtl1[comp_j] == gtl0[comp_i])
-                            {
-                              cell_matrix(i, j) += fev0.shape_value(i, q0) *
-                                                   fev1.shape_value(j, q1) *
-                                                   kernel_values[q1] *
-                                                   fev0.JxW(q0) * fev1.JxW(q1);
-                            }
-                        }
-                  }
-              }
+          for (unsigned int j = 0; j < fe1.dofs_per_cell; ++j)
+            {
+              const auto comp_j = fe1.system_to_component_index(j).first;
+
+              // First compute the part of the integral that does not
+              // depend on i
+              typename Matrix::value_type sum_q1 = {};
+              for (unsigned int q1 = 0; q1 < quadrature1.size(); ++q1)
+                sum_q1 +=
+                  fev1.shape_value(j, q1) * kernel_values[q1] * fev1.JxW(q1);
+              sum_q1 *= fev0.JxW(q0);
+
+              // Now compute the main integral with the sum over q1 already
+              // completed - this gives a cubic complexity as usual rather
+              // than a quartic one with naive loops
+              for (unsigned int i = 0; i < fe0.dofs_per_cell; ++i)
+                {
+                  const auto comp_i = fe0.system_to_component_index(i).first;
+                  if (gtl0[comp_i] != numbers::invalid_unsigned_int &&
+                      gtl1[comp_j] == gtl0[comp_i])
+                    cell_matrix(i, j) += fev0.shape_value(i, q0) * sum_q1;
+                }
+            }
         }
 
       constraints0.distribute_local_to_global(cell_matrix,