From 6e250cd147deb3dc026879cc94e1cf93bee7c16f Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sun, 5 Jan 2025 23:53:25 -0500 Subject: [PATCH] support DoF coupling in FESystem --- source/fe/fe_system.cc | 55 +++++ ...le_q_iso_q1.cc => assemble_q_iso_q1_01.cc} | 0 ..._q1.output => assemble_q_iso_q1_01.output} | 0 tests/fe/assemble_q_iso_q1_02.cc | 230 ++++++++++++++++++ tests/fe/assemble_q_iso_q1_02.output | 4 + tests/fe/deformed_projection.h | 11 +- tests/fe/local_sparsity_01.cc | 67 +++++ tests/fe/local_sparsity_01.output | 36 +++ tests/fe/local_sparsity_q_iso_q1.cc | 2 +- 9 files changed, 395 insertions(+), 10 deletions(-) rename tests/fe/{assemble_q_iso_q1.cc => assemble_q_iso_q1_01.cc} (100%) rename tests/fe/{assemble_q_iso_q1.output => assemble_q_iso_q1_01.output} (100%) create mode 100644 tests/fe/assemble_q_iso_q1_02.cc create mode 100644 tests/fe/assemble_q_iso_q1_02.output create mode 100644 tests/fe/local_sparsity_01.cc create mode 100644 tests/fe/local_sparsity_01.output diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 80872c440d..4de69bdb6d 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -2009,6 +2009,61 @@ FESystem::initialize( Assert(index == this->n_dofs_per_line(), ExcInternalError()); }); + + // Compute local_dof_sparsity_pattern if any of our base elements contains a + // non-empty one (empty denotes the default of all DoFs coupling within a + // cell). Note the we currently only handle coupling within a base element and + // not between two different base elements. Handling the latter could be + // doable if the underlying element happens to be identical, but we currently + // have no functionality to compute the coupling between different elements + // with a pattern (for example FE_Q_iso_Q1 with different degrees). + { + // Does any of our base elements not couple all DoFs? + const bool have_nonempty = [&]() -> bool { + for (unsigned int b = 0; b < this->n_base_elements(); ++b) + { + if (!this->base_element(b).get_local_dof_sparsity_pattern().empty() && + (this->element_multiplicity(b) > 0)) + return true; + } + return false; + }(); + + if (have_nonempty) + { + this->local_dof_sparsity_pattern.reinit(this->n_dofs_per_cell(), + this->n_dofs_per_cell()); + + // by default, everything couples: + this->local_dof_sparsity_pattern.fill(true); + + // Find shape functions within the same base element. If we do, grab the + // coupling from that base element pattern: + for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) + for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) + { + const auto vi = this->system_to_base_index(i); + const auto vj = this->system_to_base_index(j); + + const auto base_index_i = vi.first.first; + const auto base_index_j = vj.first.first; + if (base_index_i == base_index_j) + { + const auto shape_index_i = vi.second; + const auto shape_index_j = vj.second; + + const auto &pattern = this->base_element(base_index_i) + .get_local_dof_sparsity_pattern(); + + if (!pattern.empty()) + this->local_dof_sparsity_pattern(i, j) = + pattern(shape_index_i, shape_index_j); + } + } + } + } + + // wait for all of this to finish init_tasks.join_all(); } diff --git a/tests/fe/assemble_q_iso_q1.cc b/tests/fe/assemble_q_iso_q1_01.cc similarity index 100% rename from tests/fe/assemble_q_iso_q1.cc rename to tests/fe/assemble_q_iso_q1_01.cc diff --git a/tests/fe/assemble_q_iso_q1.output b/tests/fe/assemble_q_iso_q1_01.output similarity index 100% rename from tests/fe/assemble_q_iso_q1.output rename to tests/fe/assemble_q_iso_q1_01.output diff --git a/tests/fe/assemble_q_iso_q1_02.cc b/tests/fe/assemble_q_iso_q1_02.cc new file mode 100644 index 0000000000..fbf5a332f3 --- /dev/null +++ b/tests/fe/assemble_q_iso_q1_02.cc @@ -0,0 +1,230 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2013 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +// Assemble a 1d,2d,3d vector Poisson problem with FE_Q_iso_Q1 elements to +// check that the sparsity pattern is computed correctly. + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + + + +template +class Step4 +{ +public: + Step4(); + void + run(); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_preconditioner(); + + Triangulation triangulation; + + FESystem fe_precondition; + DoFHandler dof_handler_precondition; + + AffineConstraints constraints; + + SparsityPattern prec_pattern; + SparseMatrix preconditioner_matrix; + + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; +}; + + + +template +double +RightHandSide::value(const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i = 0; i < dim; ++i) + return_value += 4 * std::pow(p[i], 4); + + return return_value; +} + + +template +Step4::Step4() + : fe_precondition(FE_Q_iso_Q1(2), dim) + , dof_handler_precondition(triangulation) +{} + + +template +void +Step4::make_grid() +{ + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(1); +} + + + +template +void +Step4::setup_system() +{ + dof_handler_precondition.distribute_dofs(fe_precondition); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler_precondition, + 0, + Functions::ZeroFunction(dim), + constraints); + constraints.close(); + + { + DynamicSparsityPattern c_sparsity(dof_handler_precondition.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler_precondition, + c_sparsity, + constraints, + false); + prec_pattern.copy_from(c_sparsity); + preconditioner_matrix.reinit(prec_pattern); + } + + system_rhs.reinit(dof_handler_precondition.n_dofs()); +} + + + +template +void +Step4::assemble_preconditioner() +{ + QIterated quadrature_formula(QGauss<1>(2), 3); + + const RightHandSide right_hand_side; + + FEValues fe_values(fe_precondition, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe_precondition.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler_precondition.begin_active(), + endc = dof_handler_precondition.end(); + + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += + (fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point)); + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + local_dof_indices, + preconditioner_matrix); + } + preconditioner_matrix.compress(VectorOperation::add); +} + + + +template +void +Step4::run() +{ + deallog.push("dim " + std::to_string(dim)); + + make_grid(); + + setup_system(); + assemble_preconditioner(); + deallog << "nnz: " << preconditioner_matrix.n_nonzero_elements() << std::endl; + + deallog.pop(); +} + + +int +main(int argc, char **argv) +{ + initlog(true); + + { + Step4<1> test; + test.run(); + } + { + Step4<2> test; + test.run(); + } + { + Step4<3> test; + test.run(); + } +} diff --git a/tests/fe/assemble_q_iso_q1_02.output b/tests/fe/assemble_q_iso_q1_02.output new file mode 100644 index 0000000000..aa054de986 --- /dev/null +++ b/tests/fe/assemble_q_iso_q1_02.output @@ -0,0 +1,4 @@ + +DEAL:dim 1::nnz: 11 +DEAL:dim 2::nnz: 228 +DEAL:dim 3::nnz: 3381 diff --git a/tests/fe/deformed_projection.h b/tests/fe/deformed_projection.h index 9f4cb3061e..6d07432acd 100644 --- a/tests/fe/deformed_projection.h +++ b/tests/fe/deformed_projection.h @@ -365,17 +365,10 @@ create_mass_matrix(const Mapping &mapping, } } } - // transfer everything into the - // global object. lock the - // matrix meanwhile + // transfer everything into the global object for (unsigned int i = 0; i < dofs_per_cell; ++i) for (unsigned int j = 0; j < dofs_per_cell; ++j) - if ((n_components == 1) || (cell_matrix(i, j) != 0.0)) - /* - (fe.system_to_component_index(i).first == - fe.system_to_component_index(j).first)) - */ - matrix.add(dof_indices[i], dof_indices[j], cell_matrix(i, j)); + matrix.add(dof_indices[i], dof_indices[j], cell_matrix(i, j)); for (unsigned int i = 0; i < dofs_per_cell; ++i) rhs_vector(dof_indices[i]) += cell_vector(i); diff --git a/tests/fe/local_sparsity_01.cc b/tests/fe/local_sparsity_01.cc new file mode 100644 index 0000000000..7e9396baf0 --- /dev/null +++ b/tests/fe/local_sparsity_01.cc @@ -0,0 +1,67 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2024 - 2025 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +// Check FiniteElement::get_local_dof_sparsity_pattern() for FESystem. + +#include +#include + +#include "deal.II/fe/fe_system.h" +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +test(const FiniteElement &fe, const unsigned int n_subdivisions = 1) +{ + const auto &pattern = fe.get_local_dof_sparsity_pattern(); + deallog << fe.get_name() << ":\n"; + + if (!pattern.empty()) + { + for (unsigned int i = 0; i < pattern.size(0); ++i) + { + for (unsigned int j = 0; j < pattern.size(1); ++j) + deallog << (pattern(i, j) == true ? "X" : "."); + deallog << "\n"; + } + deallog << std::endl; + } +} + +int +main() +{ + initlog(); + + // The simplest case is a single iso Q1: + test<1>(FE_Q_iso_Q1<1>(2), 1); + // The 2 iso Q1 elements couple using their pattern: + test<1>(FESystem<1, 1>(FE_DGQ<1>(1), 1, FE_Q_iso_Q1<1>(2), 2), 1); + // The coupling between the first two to the third copy is a full coupling + // currently, because we don't detect this yet: + test<1>(FESystem<1, 1>(FE_Q_iso_Q1<1>(2), 2, FE_Q_iso_Q1<1>(2), 1), 1); + // Different iso_Q1 degrees always couple fully (off diagonal blocks): + test<1>(FESystem<1, 1>(FE_Q_iso_Q1<1>(2), 1, FE_Q_iso_Q1<1>(3), 1), 1); +} diff --git a/tests/fe/local_sparsity_01.output b/tests/fe/local_sparsity_01.output new file mode 100644 index 0000000000..b5a8236e7c --- /dev/null +++ b/tests/fe/local_sparsity_01.output @@ -0,0 +1,36 @@ + +DEAL::FE_Q_iso_Q1<1>(2): +X.X +.XX +XXX + +DEAL::FESystem<1>[FE_DGQ<1>(1)-FE_Q_iso_Q1<1>(2)^2]: +XX..XXXX +XX..XXXX +..XXXXXX +..XXXXXX +XXXXXXXX +XXXXXXXX +XXXXXXXX +XXXXXXXX + +DEAL::FESystem<1>[FE_Q_iso_Q1<1>(2)^2-FE_Q_iso_Q1<1>(2)]: +XXX..XXXX +XXX..XXXX +XXXXX.XXX +..XXXXXXX +..XXXXXXX +XX.XXXXXX +XXXXXXXXX +XXXXXXXXX +XXXXXXXXX + +DEAL::FESystem<1>[FE_Q_iso_Q1<1>(2)-FE_Q_iso_Q1<1>(3)]: +XX.XXXX +XXX.XX. +.XXXXXX +X.XXX.X +XXXXXXX +XXX.XXX +X.XXXXX + diff --git a/tests/fe/local_sparsity_q_iso_q1.cc b/tests/fe/local_sparsity_q_iso_q1.cc index 5698e80f5d..c01f55e59c 100644 --- a/tests/fe/local_sparsity_q_iso_q1.cc +++ b/tests/fe/local_sparsity_q_iso_q1.cc @@ -12,7 +12,7 @@ // // ------------------------------------------------------------------------ -// Assemble the sparsity pattern for a Q1 finite element space and a an iso Q +// Assemble the sparsity pattern for a Q1 finite element space and an iso Q // element space with the same number of subdivisions, and check we get the same // sparsity patterns. -- 2.39.5