From: Timo Heister Date: Thu, 6 Mar 2025 02:00:47 +0000 (-0500) Subject: fix local sparsity for Q_iso_Q1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F18198%2Fhead;p=dealii.git fix local sparsity for Q_iso_Q1 Add local sparsity pattern logic to FE_Q_iso_Q1 also when custom support points are given. --- diff --git a/source/fe/fe_q_iso_q1.cc b/source/fe/fe_q_iso_q1.cc index 31bd9ee88b..ef626fd23c 100644 --- a/source/fe/fe_q_iso_q1.cc +++ b/source/fe/fe_q_iso_q1.cc @@ -30,6 +30,49 @@ DEAL_II_NAMESPACE_OPEN +namespace internal +{ + template + Table<2, bool> + make_local_sparsity_pattern(unsigned int n_points_1d) + { + const unsigned int n_per_cell = Utilities::pow(n_points_1d, dim); + + const std::vector R = + FETools::lexicographic_to_hierarchic_numbering(n_points_1d - 1); + + Table<2, bool> table; + table.reinit(n_per_cell, n_per_cell); + table.fill(false); + + const int N1d = n_points_1d; + for (unsigned int i = 0; i < n_per_cell; ++i) + for (unsigned int j = 0; j < n_per_cell; ++j) + { + // compute l1 distance: + int distance = 0; + + int xi = i; + int xj = j; + for (unsigned int d = 0; d < dim; ++d) + { + int current_distance = std::abs((xi % N1d) - (xj % N1d)); + xi /= N1d; + xj /= N1d; + distance = std::max(distance, current_distance); + if (distance > 1) + break; + } + + if (distance <= 1) + table(R[i], R[j]) = true; + } + + return table; + } +} // namespace internal + + template FE_Q_iso_Q1::FE_Q_iso_Q1(const unsigned int subdivisions) @@ -52,42 +95,8 @@ FE_Q_iso_Q1::FE_Q_iso_Q1(const unsigned int subdivisions) const QIterated<1> points(trapez, subdivisions); this->initialize(points.get_points()); - - { - const std::vector R = - FETools::lexicographic_to_hierarchic_numbering(subdivisions); - - this->local_dof_sparsity_pattern.reinit(this->n_dofs_per_cell(), - this->n_dofs_per_cell()); - - { - this->local_dof_sparsity_pattern.fill(false); - - const unsigned int N = Utilities::pow(points.size(), dim); - const int N1d = points.size(); - for (unsigned int i = 0; i < N; ++i) - for (unsigned int j = 0; j < N; ++j) - { - // compute l1 distance: - int distance = 0; - - int xi = i; - int xj = j; - for (unsigned int d = 0; d < dim; ++d) - { - int current_distance = std::abs((xi % N1d) - (xj % N1d)); - xi /= N1d; - xj /= N1d; - distance = std::max(distance, current_distance); - if (distance > 1) - break; - } - - if (distance <= 1) - this->local_dof_sparsity_pattern(R[i], R[j]) = true; - } - } - } + this->local_dof_sparsity_pattern = + internal::make_local_sparsity_pattern(subdivisions + 1); } @@ -110,6 +119,8 @@ FE_Q_iso_Q1::FE_Q_iso_Q1( "subelements")); this->initialize(support_points); + this->local_dof_sparsity_pattern = + internal::make_local_sparsity_pattern(support_points.size()); }