From 2fd63ac5e23af30addd5ce6cfb250f26862e0837 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 1 May 2020 07:28:02 +0200 Subject: [PATCH] Remove quartic complexity from create_coupling_mass_matrix --- source/non_matching/coupling.cc | 44 +++++++++++++++++---------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/source/non_matching/coupling.cc b/source/non_matching/coupling.cc index c71f3a57b1..93bb108b8c 100644 --- 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, -- 2.39.5