From da351d22df7e85205cd5448f895d3cb97c62c8c7 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 17 Aug 2021 09:19:51 +0200 Subject: [PATCH] Enable MGTwoLevelTransfer::reinit_geometric_transfer() for refinment levels --- .../multigrid/mg_transfer_global_coarsening.h | 11 +- .../mg_transfer_global_coarsening.templates.h | 285 +++++++++++++----- .../multigrid_a_01.cc | 4 +- .../multigrid_a_02.cc | 158 ++++++++++ ..._a_02.mpirun=1.with_trilinos=true_1.output | 10 + .../multigrid_p_02.cc | 3 +- .../multigrid_util.h | 19 +- 7 files changed, 406 insertions(+), 84 deletions(-) create mode 100644 tests/multigrid-global-coarsening/multigrid_a_02.cc create mode 100644 tests/multigrid-global-coarsening/multigrid_a_02.mpirun=1.with_trilinos=true_1.output diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h index 1edbf0d7bf..f2ca2ebcf4 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h @@ -193,10 +193,13 @@ public: * can be only performed on active levels. */ void - reinit_geometric_transfer(const DoFHandler & dof_handler_fine, - const DoFHandler & dof_handler_coarse, - const AffineConstraints &constraint_fine, - const AffineConstraints &constraint_coarse); + reinit_geometric_transfer( + const DoFHandler & dof_handler_fine, + const DoFHandler & dof_handler_coarse, + const AffineConstraints &constraint_fine, + const AffineConstraints &constraint_coarse, + const unsigned int mg_level_fine = numbers::invalid_unsigned_int, + const unsigned int mg_level_coarse = numbers::invalid_unsigned_int); /** * Set up polynomial coarsening between the given DoFHandler objects ( diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index 77da761097..3ec986903f 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -40,6 +40,7 @@ #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -681,9 +682,15 @@ namespace internal ghost_indices.end()), ghost_indices.end()); - this->is_extended_locally_owned = dof_handler_dst.locally_owned_dofs(); + this->is_extended_locally_owned = + mg_level_fine == numbers::invalid_unsigned_int ? + dof_handler_dst.locally_owned_dofs() : + dof_handler_dst.locally_owned_mg_dofs(mg_level_fine); - this->is_extended_ghosts = IndexSet(dof_handler_dst.n_dofs()); + this->is_extended_ghosts = + IndexSet(mg_level_fine == numbers::invalid_unsigned_int ? + dof_handler_dst.n_dofs() : + dof_handler_dst.n_dofs(mg_level_fine)); this->is_extended_ghosts.add_indices(ghost_indices.begin(), ghost_indices.end()); this->is_extended_ghosts.subtract_set(this->is_extended_locally_owned); @@ -852,10 +859,16 @@ namespace internal { public: GlobalCoarseningFineDoFHandlerView(const DoFHandler &dof_handler_dst, - const DoFHandler &dof_handler_src) - : FineDoFHandlerView< - dim>(dof_handler_dst, dof_handler_src, numbers::invalid_unsigned_int /*global coarsening only possible on active levels*/) + const DoFHandler &dof_handler_src, + const unsigned int mg_level_fine, + const unsigned int mg_level_coarse) + : FineDoFHandlerView(dof_handler_dst, dof_handler_src, mg_level_fine) { + Assert((mg_level_fine == numbers::invalid_unsigned_int && + mg_level_coarse == numbers::invalid_unsigned_int) || + (mg_level_coarse + 1 == mg_level_fine), + ExcNotImplemented()); + // get reference to triangulations const auto &tria_dst = dof_handler_dst.get_triangulation(); const auto &tria_src = dof_handler_src.get_triangulation(); @@ -865,28 +878,46 @@ namespace internal IndexSet is_dst_remote(this->cell_id_translator.size()); IndexSet is_src_locally_owned(this->cell_id_translator.size()); - for (const auto &cell : tria_dst.active_cell_iterators()) - if (cell->is_locally_owned()) - is_dst_locally_owned.add_index( - this->cell_id_translator.translate(cell)); + const auto fine_operation = [&](const auto &cell) { + is_dst_locally_owned.add_index( + this->cell_id_translator.translate(cell)); + }; + const auto coarse_operation = [&](const auto &cell) { + is_src_locally_owned.add_index( + this->cell_id_translator.translate(cell)); - for (const auto &cell : tria_src.active_cell_iterators()) - if (cell->is_locally_owned()) - { - is_src_locally_owned.add_index( - this->cell_id_translator.translate(cell)); - is_dst_remote.add_index(this->cell_id_translator.translate(cell)); + // in the case of global coarsening identity transfer is possible + if (mg_level_coarse == numbers::invalid_unsigned_int) + is_dst_remote.add_index(this->cell_id_translator.translate(cell)); - if (cell->level() + 1u == tria_dst.n_global_levels()) - continue; + if (cell->level() + 1u == tria_dst.n_global_levels()) + return; - for (unsigned int i = 0; - i < GeometryInfo::max_children_per_cell; - ++i) - is_dst_remote.add_index( - this->cell_id_translator.translate(cell, i)); + for (unsigned int i = 0; i < GeometryInfo::max_children_per_cell; + ++i) + is_dst_remote.add_index(this->cell_id_translator.translate(cell, i)); + }; + + const auto loop = [&](const auto &tria, + const auto level, + const auto &op) { + if (level == numbers::invalid_unsigned_int) + { + for (const auto &cell : tria.active_cell_iterators()) + if (!cell->is_artificial() && cell->is_locally_owned()) + op(cell); } + else + { + for (const auto &cell : tria.cell_iterators_on_level(level)) + if (cell->level_subdomain_id() == tria.locally_owned_subdomain()) + op(cell); + } + }; + + loop(tria_dst, mg_level_fine, fine_operation); + loop(tria_src, mg_level_coarse, coarse_operation); this->reinit(is_dst_locally_owned, is_dst_remote, @@ -960,12 +991,16 @@ namespace internal precompute_constraints( const DoFHandler & dof_handler_coarse, const dealii::AffineConstraints &constraint_coarse, + const unsigned int mg_level_coarse, MGTwoLevelTransfer> &transfer) { std::vector dependencies; - const auto &locally_owned_dofs = dof_handler_coarse.locally_owned_dofs(); + const auto &locally_owned_dofs = + mg_level_coarse == numbers::invalid_unsigned_int ? + dof_handler_coarse.locally_owned_dofs() : + dof_handler_coarse.locally_owned_mg_dofs(mg_level_coarse); for (const auto i : locally_owned_dofs) { @@ -1134,7 +1169,7 @@ namespace internal else { for (const auto &cell : - dof_handler_fine.mg_cell_iterators_on_level(mg_level_fine)) + dof_handler_fine.cell_iterators_on_level(mg_level_fine)) { if (!cell->is_locally_owned_on_level()) continue; @@ -1175,14 +1210,38 @@ namespace internal const DoFHandler & dof_handler_coarse, const dealii::AffineConstraints &constraint_fine, const dealii::AffineConstraints &constraint_coarse, + const unsigned int mg_level_fine, + const unsigned int mg_level_coarse, MGTwoLevelTransfer> &transfer) { + Assert((mg_level_fine == numbers::invalid_unsigned_int && + mg_level_coarse == numbers::invalid_unsigned_int) || + (mg_level_coarse + 1 == mg_level_fine), + ExcNotImplemented()); + + if (mg_level_fine != numbers::invalid_unsigned_int) + AssertIndexRange(mg_level_fine, + MGTools::max_level_for_coarse_mesh( + dof_handler_fine.get_triangulation()) + + 1); + + if (mg_level_coarse != numbers::invalid_unsigned_int) + AssertIndexRange(mg_level_coarse, + MGTools::max_level_for_coarse_mesh( + dof_handler_coarse.get_triangulation()) + + 1); + const GlobalCoarseningFineDoFHandlerView view(dof_handler_fine, - dof_handler_coarse); + dof_handler_coarse, + mg_level_fine, + mg_level_coarse); // copy constrain matrix; TODO: why only for the coarse level? - precompute_constraints(dof_handler_coarse, constraint_coarse, transfer); + precompute_constraints(dof_handler_coarse, + constraint_coarse, + mg_level_coarse, + transfer); // gather ranges for active FE indices on both fine and coarse dofhandlers std::array min_active_fe_indices = { @@ -1190,23 +1249,53 @@ namespace internal std::numeric_limits::max()}}; std::array max_active_fe_indices = {{0, 0}}; - for (const auto &cell : dof_handler_fine.active_cell_iterators()) - if (cell->is_locally_owned()) - { - min_active_fe_indices[0] = - std::min(min_active_fe_indices[0], cell->active_fe_index()); - max_active_fe_indices[0] = - std::max(max_active_fe_indices[0], cell->active_fe_index()); - } + if (mg_level_fine == numbers::invalid_unsigned_int) + { + for (const auto &cell : dof_handler_fine.active_cell_iterators()) + if (cell->is_locally_owned()) + { + min_active_fe_indices[0] = + std::min(min_active_fe_indices[0], cell->active_fe_index()); + max_active_fe_indices[0] = + std::max(max_active_fe_indices[0], cell->active_fe_index()); + } + } + else + { + for (const auto &cell : + dof_handler_fine.mg_cell_iterators_on_level(mg_level_fine)) + if (cell->is_locally_owned_on_level()) + { + min_active_fe_indices[0] = + std::min(min_active_fe_indices[0], cell->active_fe_index()); + max_active_fe_indices[0] = + std::max(max_active_fe_indices[0], cell->active_fe_index()); + } + } - for (const auto &cell : dof_handler_coarse.active_cell_iterators()) - if (cell->is_locally_owned()) - { - min_active_fe_indices[1] = - std::min(min_active_fe_indices[1], cell->active_fe_index()); - max_active_fe_indices[1] = - std::max(max_active_fe_indices[1], cell->active_fe_index()); - } + if (mg_level_fine == numbers::invalid_unsigned_int) + { + for (const auto &cell : dof_handler_coarse.active_cell_iterators()) + if (cell->is_locally_owned()) + { + min_active_fe_indices[1] = + std::min(min_active_fe_indices[1], cell->active_fe_index()); + max_active_fe_indices[1] = + std::max(max_active_fe_indices[1], cell->active_fe_index()); + } + } + else + { + for (const auto &cell : + dof_handler_coarse.mg_cell_iterators_on_level(mg_level_coarse)) + if (cell->is_locally_owned_on_level()) + { + min_active_fe_indices[1] = + std::min(min_active_fe_indices[1], cell->active_fe_index()); + max_active_fe_indices[1] = + std::max(max_active_fe_indices[1], cell->active_fe_index()); + } + } const auto comm = dof_handler_fine.get_communicator(); @@ -1249,12 +1338,19 @@ namespace internal // ... coarse mesh (needed since user vector might be const) { IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs(dof_handler_coarse, - locally_relevant_dofs); + + if (mg_level_coarse == numbers::invalid_unsigned_int) + DoFTools::extract_locally_relevant_dofs(dof_handler_coarse, + locally_relevant_dofs); + else + DoFTools::extract_locally_relevant_level_dofs( + dof_handler_coarse, mg_level_coarse, locally_relevant_dofs); transfer.partitioner_coarse = std::make_shared( - dof_handler_coarse.locally_owned_dofs(), + mg_level_coarse == numbers::invalid_unsigned_int ? + dof_handler_coarse.locally_owned_dofs() : + dof_handler_coarse.locally_owned_mg_dofs(mg_level_coarse), locally_relevant_dofs, dof_handler_coarse.get_communicator()); transfer.vec_coarse.reinit(transfer.partitioner_coarse); @@ -1267,31 +1363,57 @@ namespace internal // helper function: to process the fine level cells; function @p fu_non_refined is // performed on cells that are not refined and @fu_refined is performed on // children of cells that are refined - const auto process_cells = - [&view, &dof_handler_coarse](const auto &fu_non_refined, + const auto process_cells = [&](const auto &fu_non_refined, const auto &fu_refined) { - // loop over all active locally-owned cells - for (const auto &cell_coarse : - dof_handler_coarse.active_cell_iterators()) - { - if (!cell_coarse->is_locally_owned()) - continue; + // loop over all active locally-owned cells + if (mg_level_coarse == numbers::invalid_unsigned_int) + { + for (const auto &cell_coarse : + dof_handler_coarse.active_cell_iterators()) + { + if (cell_coarse->is_artificial() == true || + cell_coarse->is_locally_owned() == false) + continue; - // get a reference to the equivalent cell on the fine - // triangulation - const auto cell_coarse_on_fine_mesh = view.get_cell(cell_coarse); - - // check if cell has children - if (cell_coarse_on_fine_mesh.has_children()) - // ... cell has children -> process children - for (unsigned int c = 0; - c < GeometryInfo::max_children_per_cell; - c++) - fu_refined(cell_coarse, view.get_cell(cell_coarse, c), c); - else // ... cell has no children -> process cell - fu_non_refined(cell_coarse, cell_coarse_on_fine_mesh); - } - }; + // get a reference to the equivalent cell on the fine + // triangulation + const auto cell_coarse_on_fine_mesh = + view.get_cell(cell_coarse); + + // check if cell has children + if (cell_coarse_on_fine_mesh.has_children()) + // ... cell has children -> process children + for (unsigned int c = 0; + c < GeometryInfo::max_children_per_cell; + c++) + fu_refined(cell_coarse, view.get_cell(cell_coarse, c), c); + else // ... cell has no children -> process cell + fu_non_refined(cell_coarse, cell_coarse_on_fine_mesh); + } + } + else + { + for (const auto &cell_coarse : + dof_handler_coarse.mg_cell_iterators_on_level(mg_level_coarse)) + { + if (cell_coarse->is_locally_owned_on_level() == false) + continue; + + // get a reference to the equivalent cell on the fine + // triangulation + // const auto cell_coarse_on_fine_mesh = + // view.get_cell(cell_coarse); + + // check if cell has children + if (cell_coarse->has_children()) + // ... cell has children -> process children + for (unsigned int c = 0; + c < GeometryInfo::max_children_per_cell; + c++) + fu_refined(cell_coarse, view.get_cell(cell_coarse, c), c); + } + } + }; // set up two mg-schemes // (0) no refinement -> identity @@ -1395,7 +1517,11 @@ namespace internal [&](const auto &cell_coarse, const auto &cell_fine) { // parent { - cell_coarse->get_dof_indices(local_dof_indices); + if (mg_level_coarse == numbers::invalid_unsigned_int) + cell_coarse->get_dof_indices(local_dof_indices); + else + cell_coarse->get_mg_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < transfer.schemes[0].dofs_per_cell_coarse; i++) @@ -1427,7 +1553,11 @@ namespace internal // parent (only once at the beginning) if (c == 0) { - cell_coarse->get_dof_indices(local_dof_indices); + if (mg_level_coarse == numbers::invalid_unsigned_int) + cell_coarse->get_dof_indices(local_dof_indices); + else + cell_coarse->get_mg_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < transfer.schemes[1].dofs_per_cell_coarse; i++) @@ -1569,7 +1699,7 @@ namespace internal // compute weights globally LinearAlgebra::distributed::Vector weight_vector; compute_weights(dof_handler_fine, - numbers::invalid_unsigned_int /*active level*/, + mg_level_fine, constraint_fine, transfer, weight_vector); @@ -1742,7 +1872,10 @@ namespace internal dof_handler_fine.get_fe(fe_index_pair.first.second).degree; } - precompute_constraints(dof_handler_coarse, constraint_coarse, transfer); + precompute_constraints(dof_handler_coarse, + constraint_coarse, + mg_level_coarse, + transfer); const auto comm = dof_handler_coarse.get_communicator(); { @@ -2767,13 +2900,17 @@ MGTwoLevelTransfer>:: reinit_geometric_transfer(const DoFHandler & dof_handler_fine, const DoFHandler & dof_handler_coarse, const AffineConstraints &constraint_fine, - const AffineConstraints &constraint_coarse) + const AffineConstraints &constraint_coarse, + const unsigned int mg_level_fine, + const unsigned int mg_level_coarse) { internal::MGTwoLevelTransferImplementation::reinit_geometric_transfer( dof_handler_fine, dof_handler_coarse, constraint_fine, constraint_coarse, + mg_level_fine, + mg_level_coarse, *this); } @@ -2886,6 +3023,8 @@ MGTwoLevelTransfer>::reinit( dof_handler_coarse, constraint_fine, constraint_coarse, + mg_level_fine, + mg_level_coarse, *this); } diff --git a/tests/multigrid-global-coarsening/multigrid_a_01.cc b/tests/multigrid-global-coarsening/multigrid_a_01.cc index 62c04637a3..302aafb0a9 100644 --- a/tests/multigrid-global-coarsening/multigrid_a_01.cc +++ b/tests/multigrid-global-coarsening/multigrid_a_01.cc @@ -15,8 +15,8 @@ /** - * Test p-multigrid for a uniformly refined mesh both for simplex and - * hypercube mesh. + * Test global-coarsening multigrid for a uniformly refined mesh both for + * simplex and hypercube mesh. */ #include "multigrid_util.h" diff --git a/tests/multigrid-global-coarsening/multigrid_a_02.cc b/tests/multigrid-global-coarsening/multigrid_a_02.cc new file mode 100644 index 0000000000..b8ccd60ae2 --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_a_02.cc @@ -0,0 +1,158 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +/** + * Test global-coarsening multigrid so that it works also on refinement levels. + */ + +#include "multigrid_util.h" + +template +void +test(const unsigned int n_refinements, const unsigned int fe_degree_fine) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + parallel::distributed::Triangulation tria( + MPI_COMM_WORLD, + Triangulation::MeshSmoothing::none, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + GridGenerator::subdivided_hyper_cube(tria, 2); + tria.refine_global(n_refinements); + + const unsigned int min_level = 0; + const unsigned int max_level = n_refinements; + + MGLevelObject> constraints(min_level, max_level); + MGLevelObject> transfers(min_level, + max_level); + MGLevelObject> operators(min_level, max_level); + + std::unique_ptr> fe; + std::unique_ptr> quad; + std::unique_ptr> mapping; + + fe = std::make_unique>(fe_degree_fine); + quad = std::make_unique>(fe_degree_fine + 1); + mapping = std::make_unique>(FE_Q(1)); + + // set up dofhandler + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(*fe); + dof_handler.distribute_mg_dofs(); + + std::set dirichlet_boundary; + dirichlet_boundary.insert(0); + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, + dirichlet_boundary); + + // set up levels + for (auto l = min_level; l <= max_level; ++l) + { + auto &constraint = constraints[l]; + auto &op = operators[l]; + + // set up constraints + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, + l, + relevant_dofs); + constraint.reinit(relevant_dofs); + constraint.add_lines(mg_constrained_dofs.get_boundary_indices(l)); + constraint.close(); + + constraint.print(std::cout); + + // set up operator + op.reinit(*mapping, dof_handler, *quad, constraint, l); + } + + // set up transfer operator + for (unsigned int l = min_level; l < max_level; ++l) + transfers[l + 1].reinit_geometric_transfer( + dof_handler, dof_handler, constraints[l + 1], constraints[l], l + 1, l); + + MGTransferGlobalCoarsening transfer( + transfers, + [&](const auto l, auto &vec) { operators[l].initialize_dof_vector(vec); }); + + GMGParameters mg_data; // TODO + + VectorType dst, src; + operators[max_level].initialize_dof_vector(dst); + operators[max_level].initialize_dof_vector(src); + + operators[max_level].rhs(src); + + ReductionControl solver_control( + mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false); + + mg_solve(solver_control, + dst, + src, + mg_data, + dof_handler, + operators[max_level], + operators, + transfer); + + constraints[max_level].distribute(dst); + + deallog << dim << " " << fe_degree_fine << " " << n_refinements << " " + << solver_control.last_step() << std::endl; + + return; + + static unsigned int counter = 0; + + MGLevelObject results(min_level, max_level); + + transfer.interpolate_to_mg(dof_handler, results, dst); + + for (unsigned int l = min_level; l <= max_level; ++l) + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector( + results[l], + "solution", + DataOut_DoFData::DataVectorType::type_dof_data); + data_out.build_patches(*mapping, 2); + + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter) + "." + std::to_string(l) + + ".vtk"); + data_out.write_vtk(output); + } + + counter++; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + deallog.precision(8); + + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 2; degree <= 4; ++degree) + test<2>(n_refinements, degree); +} diff --git a/tests/multigrid-global-coarsening/multigrid_a_02.mpirun=1.with_trilinos=true_1.output b/tests/multigrid-global-coarsening/multigrid_a_02.mpirun=1.with_trilinos=true_1.output new file mode 100644 index 0000000000..e460d32df9 --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_a_02.mpirun=1.with_trilinos=true_1.output @@ -0,0 +1,10 @@ + +DEAL:0::2 2 2 3 +DEAL:0::2 3 2 3 +DEAL:0::2 4 2 3 +DEAL:0::2 2 3 3 +DEAL:0::2 3 3 3 +DEAL:0::2 4 3 3 +DEAL:0::2 2 4 3 +DEAL:0::2 3 4 3 +DEAL:0::2 4 4 3 diff --git a/tests/multigrid-global-coarsening/multigrid_p_02.cc b/tests/multigrid-global-coarsening/multigrid_p_02.cc index 00aa4430ed..cc30f93f23 100644 --- a/tests/multigrid-global-coarsening/multigrid_p_02.cc +++ b/tests/multigrid-global-coarsening/multigrid_p_02.cc @@ -15,8 +15,7 @@ /** - * Test p-multigrid for a uniformly and locally refined meshes both for - * simplex and hypercube mesh. + * Test p-multigrid so that it also works on refinement level 0. */ #include "multigrid_util.h" diff --git a/tests/multigrid-global-coarsening/multigrid_util.h b/tests/multigrid-global-coarsening/multigrid_util.h index e70980ae21..db312c29b2 100644 --- a/tests/multigrid-global-coarsening/multigrid_util.h +++ b/tests/multigrid-global-coarsening/multigrid_util.h @@ -96,7 +96,11 @@ public: virtual types::global_dof_index m() const { - return matrix_free.get_dof_handler().n_dofs(); + if (this->matrix_free.get_mg_level() != numbers::invalid_unsigned_int) + return this->matrix_free.get_dof_handler().n_dofs( + this->matrix_free.get_mg_level()); + else + return this->matrix_free.get_dof_handler().n_dofs(); } Number @@ -149,10 +153,19 @@ public: const auto &dof_handler = this->matrix_free.get_dof_handler(); TrilinosWrappers::SparsityPattern dsp( - dof_handler.locally_owned_dofs(), + this->matrix_free.get_mg_level() != numbers::invalid_unsigned_int ? + dof_handler.locally_owned_mg_dofs( + this->matrix_free.get_mg_level()) : + dof_handler.locally_owned_dofs(), matrix_free.get_task_info().communicator); - DoFTools::make_sparsity_pattern(dof_handler, dsp, this->constraints); + if (this->matrix_free.get_mg_level() != numbers::invalid_unsigned_int) + MGTools::make_sparsity_pattern(dof_handler, + dsp, + this->matrix_free.get_mg_level(), + this->constraints); + else + DoFTools::make_sparsity_pattern(dof_handler, dsp, this->constraints); dsp.compress(); system_matrix.reinit(dsp); -- 2.39.5