From 7a55b0985464b1263627ab78bf4d4620fe12adf3 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 3 Apr 2024 13:49:21 +0200 Subject: [PATCH] Add DoFTools::extract_level_constant_modes() --- doc/news/changes/minor/20240407Munch | 4 + include/deal.II/dofs/dof_tools.h | 10 + source/dofs/dof_renumbering.cc | 22 +- source/dofs/dof_tools.cc | 308 +++++++++++------- source/dofs/dof_tools.inst.in | 6 + tests/mpi/extract_constant_modes_04.cc | 129 ++++++++ ..._modes_04.mpirun=10.with_p4est=true.output | 8 + ...t_modes_04.mpirun=4.with_p4est=true.output | 8 + 8 files changed, 381 insertions(+), 114 deletions(-) create mode 100644 doc/news/changes/minor/20240407Munch create mode 100644 tests/mpi/extract_constant_modes_04.cc create mode 100644 tests/mpi/extract_constant_modes_04.mpirun=10.with_p4est=true.output create mode 100644 tests/mpi/extract_constant_modes_04.mpirun=4.with_p4est=true.output diff --git a/doc/news/changes/minor/20240407Munch b/doc/news/changes/minor/20240407Munch new file mode 100644 index 0000000000..999d88cca3 --- /dev/null +++ b/doc/news/changes/minor/20240407Munch @@ -0,0 +1,4 @@ +New: The new function DoFTools::extract_level_constant_modes() +allows to extract constant modes on multigrid levels. +
+(Peter Munch, 2024/04/07) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index bb380b6440..7d1d365945 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -1371,6 +1371,16 @@ namespace DoFTools extract_constant_modes(const DoFHandler &dof_handler, const ComponentMask &component_mask, std::vector> &constant_modes); + + /** + * Same as above but for multigrid levels. + */ + template + void + extract_level_constant_modes(const unsigned int level, + const DoFHandler &dof_handler, + const ComponentMask &component_mask, + std::vector> &constant_modes); /** @} */ /** diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc index f53ca729de..d4f060248b 100644 --- a/source/dofs/dof_renumbering.cc +++ b/source/dofs/dof_renumbering.cc @@ -723,12 +723,26 @@ namespace DoFRenumbering renumbering, start, end, component_order_arg, true); (void)result; - Assert(result == 0 || result == dof_handler.n_dofs(level), + // If we don't have a renumbering (i.e., when there is 1 component) then + // return + if (Utilities::MPI::max(renumbering.size(), + dof_handler.get_communicator()) == 0) + return; + + // verify that the last numbered + // degree of freedom is either + // equal to the number of degrees + // of freedom in total (the + // sequential case) or in the + // distributed case at least + // makes sense + Assert((result == dof_handler.locally_owned_mg_dofs(level).n_elements()) || + ((dof_handler.locally_owned_mg_dofs(level).n_elements() < + dof_handler.n_dofs(level)) && + (result <= dof_handler.n_dofs(level))), ExcInternalError()); - if (Utilities::MPI::max(renumbering.size(), - dof_handler.get_communicator()) > 0) - dof_handler.renumber_dofs(level, renumbering); + dof_handler.renumber_dofs(level, renumbering); } diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index 4b0c558b55..3fd1a6641c 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -178,7 +178,7 @@ namespace DoFTools // the output array dofs_by_component lists for each dof the // corresponding vector component. if the DoFHandler is based on a // parallel distributed triangulation then the output array is index by - // dof.locally_owned_dofs().index_within_set(indices[i]) + // dof_handler.locally_owned_dofs().index_within_set(indices[i]) // // if an element is non-primitive then we assign to each degree of // freedom the following component: @@ -190,16 +190,23 @@ namespace DoFTools // in the component_mask that corresponds to this shape function template void - get_component_association(const DoFHandler &dof, - const ComponentMask &component_mask, - std::vector &dofs_by_component) + get_component_association( + const DoFHandler &dof_handler, + const ComponentMask &component_mask, + std::vector &dofs_by_component, + const unsigned int mg_level = numbers::invalid_unsigned_int) { const dealii::hp::FECollection &fe_collection = - dof.get_fe_collection(); + dof_handler.get_fe_collection(); Assert(fe_collection.n_components() < 256, ExcNotImplemented()); - Assert(dofs_by_component.size() == dof.n_locally_owned_dofs(), - ExcDimensionMismatch(dofs_by_component.size(), - dof.n_locally_owned_dofs())); + + const auto &locally_owned_dofs = + (mg_level == numbers::invalid_unsigned_int) ? + dof_handler.locally_owned_dofs() : + dof_handler.locally_owned_mg_dofs(mg_level); + + AssertDimension(dofs_by_component.size(), + locally_owned_dofs.n_elements()); // next set up a table for the degrees of freedom on each of the // cells (regardless of the fact whether it is listed in the @@ -219,18 +226,40 @@ namespace DoFTools // then loop over all cells and do the work std::vector indices; - for (const auto &c : - dof.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) - { - const types::fe_index fe_index = c->active_fe_index(); - const unsigned int dofs_per_cell = c->get_fe().n_dofs_per_cell(); - indices.resize(dofs_per_cell); - c->get_dof_indices(indices); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - if (dof.locally_owned_dofs().is_element(indices[i])) - dofs_by_component[dof.locally_owned_dofs().index_within_set( - indices[i])] = local_component_association[fe_index][i]; - } + + const auto runner = [&](const auto &task) { + if (mg_level == numbers::invalid_unsigned_int) + { + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + indices.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(indices); + + task(cell); + } + } + else + { + for (const auto &cell : + dof_handler.cell_iterators_on_level(mg_level)) + if (cell->is_locally_owned_on_level()) + { + indices.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_mg_dof_indices(indices); + + task(cell); + } + } + }; + + runner([&](const auto &cell) { + const types::fe_index fe_index = cell->active_fe_index(); + for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell(); ++i) + if (locally_owned_dofs.is_element(indices[i])) + dofs_by_component[locally_owned_dofs.index_within_set(indices[i])] = + local_component_association[fe_index][i]; + }); } @@ -1235,104 +1264,133 @@ namespace DoFTools } - - template - void - extract_constant_modes(const DoFHandler &dof_handler, - const ComponentMask &component_mask, - std::vector> &constant_modes) + namespace internal { - // If there are no locally owned DoFs, return with an empty - // constant_modes object: - if (dof_handler.n_locally_owned_dofs() == 0) + template + void + extract_constant_modes(const DoFHandler &dof_handler, + const ComponentMask &component_mask, + const unsigned int mg_level, + std::vector> &constant_modes) + { + const auto &locally_owned_dofs = + (mg_level == numbers::invalid_unsigned_int) ? + dof_handler.locally_owned_dofs() : + dof_handler.locally_owned_mg_dofs(mg_level); + + // If there are no locally owned DoFs, return with an empty + // constant_modes object: + if (locally_owned_dofs.n_elements() == 0) + { + constant_modes = std::vector>(0); + return; + } + + const unsigned int n_components = dof_handler.get_fe(0).n_components(); + Assert(component_mask.represents_n_components(n_components), + ExcDimensionMismatch(n_components, component_mask.size())); + + std::vector dofs_by_component( + locally_owned_dofs.n_elements()); + internal::get_component_association(dof_handler, + component_mask, + dofs_by_component, + mg_level); + unsigned int n_selected_dofs = 0; + for (unsigned int i = 0; i < n_components; ++i) + if (component_mask[i] == true) + n_selected_dofs += + std::count(dofs_by_component.begin(), dofs_by_component.end(), i); + + // Find local numbering within the selected components + std::vector component_numbering( + locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int); + for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements(); + ++i) + if (component_mask[dofs_by_component[i]]) + component_numbering[i] = count++; + + // get the element constant modes and find a translation table between + // index in the constant modes and the components. + // + // TODO: We might be able to extend this also for elements which do not + // have the same constant modes, but that is messy... + const dealii::hp::FECollection &fe_collection = + dof_handler.get_fe_collection(); + std::vector> element_constant_modes; + std::vector>> + constant_mode_to_component_translation(n_components); { - constant_modes = std::vector>(0); - return; - } + unsigned int n_constant_modes = 0; + int first_non_empty_constant_mode = -1; + for (unsigned int f = 0; f < fe_collection.size(); ++f) + { + std::pair, std::vector> data = + fe_collection[f].get_constant_modes(); - const unsigned int n_components = dof_handler.get_fe(0).n_components(); - Assert(component_mask.represents_n_components(n_components), - ExcDimensionMismatch(n_components, component_mask.size())); + // Store the index of the current element if it is the first that + // has non-empty constant modes. + if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0) + { + first_non_empty_constant_mode = f; + // This is the first non-empty constant mode, so we figure out + // the translation between index in the constant modes and the + // components + for (unsigned int i = 0; i < data.second.size(); ++i) + if (component_mask[data.second[i]]) + constant_mode_to_component_translation[data.second[i]] + .emplace_back(n_constant_modes++, i); + } - std::vector dofs_by_component( - dof_handler.n_locally_owned_dofs()); - internal::get_component_association(dof_handler, - component_mask, - dofs_by_component); - unsigned int n_selected_dofs = 0; - for (unsigned int i = 0; i < n_components; ++i) - if (component_mask[i] == true) - n_selected_dofs += - std::count(dofs_by_component.begin(), dofs_by_component.end(), i); - - // Find local numbering within the selected components - const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); - std::vector component_numbering( - locally_owned_dofs.n_elements(), numbers::invalid_unsigned_int); - for (unsigned int i = 0, count = 0; i < locally_owned_dofs.n_elements(); - ++i) - if (component_mask[dofs_by_component[i]]) - component_numbering[i] = count++; - - // get the element constant modes and find a translation table between - // index in the constant modes and the components. - // - // TODO: We might be able to extend this also for elements which do not - // have the same constant modes, but that is messy... - const dealii::hp::FECollection &fe_collection = - dof_handler.get_fe_collection(); - std::vector> element_constant_modes; - std::vector>> - constant_mode_to_component_translation(n_components); - { - unsigned int n_constant_modes = 0; - int first_non_empty_constant_mode = -1; - for (unsigned int f = 0; f < fe_collection.size(); ++f) - { - std::pair, std::vector> data = - fe_collection[f].get_constant_modes(); + // Add the constant modes of this element to the list and assert + // that there are as many constant modes as for the other elements + // (or zero constant modes). + element_constant_modes.push_back(data.first); + Assert(element_constant_modes.back().n_rows() == 0 || + element_constant_modes.back().n_rows() == + element_constant_modes[first_non_empty_constant_mode] + .n_rows(), + ExcInternalError()); + } + AssertIndexRange(first_non_empty_constant_mode, fe_collection.size()); - // Store the index of the current element if it is the first that has - // non-empty constant modes. - if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0) - { - first_non_empty_constant_mode = f; - // This is the first non-empty constant mode, so we figure out the - // translation between index in the constant modes and the - // components - for (unsigned int i = 0; i < data.second.size(); ++i) - if (component_mask[data.second[i]]) - constant_mode_to_component_translation[data.second[i]] - .emplace_back(n_constant_modes++, i); - } + // Now we know the number of constant modes and resize the return vector + // accordingly + constant_modes.clear(); + constant_modes.resize(n_constant_modes, + std::vector(n_selected_dofs, false)); + } - // Add the constant modes of this element to the list and assert that - // there are as many constant modes as for the other elements (or zero - // constant modes). - element_constant_modes.push_back(data.first); - Assert( - element_constant_modes.back().n_rows() == 0 || - element_constant_modes.back().n_rows() == - element_constant_modes[first_non_empty_constant_mode].n_rows(), - ExcInternalError()); - } - AssertIndexRange(first_non_empty_constant_mode, fe_collection.size()); + // Loop over all owned cells and ask the element for the constant modes + std::vector dof_indices; - // Now we know the number of constant modes and resize the return vector - // accordingly - constant_modes.clear(); - constant_modes.resize(n_constant_modes, - std::vector(n_selected_dofs, false)); - } + const auto runner = [&](const auto &task) { + if (mg_level == numbers::invalid_unsigned_int) + { + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + dof_indices.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(dof_indices); - // Loop over all owned cells and ask the element for the constant modes - std::vector dof_indices; - for (const auto &cell : dof_handler.active_cell_iterators() | - IteratorFilters::LocallyOwnedCell()) - { - dof_indices.resize(cell->get_fe().n_dofs_per_cell()); - cell->get_dof_indices(dof_indices); + task(cell); + } + } + else + { + for (const auto &cell : + dof_handler.cell_iterators_on_level(mg_level)) + if (cell->is_locally_owned_on_level()) + { + dof_indices.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_mg_dof_indices(dof_indices); + + task(cell); + } + } + }; + runner([&](const auto &cell) { for (unsigned int i = 0; i < dof_indices.size(); ++i) if (locally_owned_dofs.is_element(dof_indices[i])) { @@ -1347,7 +1405,37 @@ namespace DoFTools element_constant_modes[cell->active_fe_index()]( indices.second, i); } - } + }); + } + } // namespace internal + + + + template + void + extract_constant_modes(const DoFHandler &dof_handler, + const ComponentMask &component_mask, + std::vector> &constant_modes) + { + internal::extract_constant_modes(dof_handler, + component_mask, + numbers::invalid_unsigned_int, + constant_modes); + } + + + + template + void + extract_level_constant_modes(const unsigned int level, + const DoFHandler &dof_handler, + const ComponentMask &component_mask, + std::vector> &constant_modes) + { + internal::extract_constant_modes(dof_handler, + component_mask, + level, + constant_modes); } diff --git a/source/dofs/dof_tools.inst.in b/source/dofs/dof_tools.inst.in index e4613113b1..864e81037a 100644 --- a/source/dofs/dof_tools.inst.in +++ b/source/dofs/dof_tools.inst.in @@ -193,6 +193,12 @@ for (deal_II_dimension : DIMENSIONS) const ComponentMask &selected_components, std::vector> &constant_modes); + template void DoFTools::extract_level_constant_modes( + const unsigned int level, + const DoFHandler &dof_handler, + const ComponentMask &selected_components, + std::vector> &constant_modes); + template void DoFTools::get_active_fe_indices( const DoFHandler &dof_handler, std::vector &active_fe_indices); diff --git a/tests/mpi/extract_constant_modes_04.cc b/tests/mpi/extract_constant_modes_04.cc new file mode 100644 index 0000000000..fa8754db81 --- /dev/null +++ b/tests/mpi/extract_constant_modes_04.cc @@ -0,0 +1,129 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2010 - 2023 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. +// +// ------------------------------------------------------------------------ + + + +// test DoFTools::extract_level_constant_modes (see also +// extract_constant_modes_01.cc). + +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "../tests.h" + + + +template +void +test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + parallel::distributed::Triangulation tr( + MPI_COMM_WORLD, + Triangulation::none, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + + std::vector sub(2); + sub[0] = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + sub[1] = 1; + GridGenerator::subdivided_hyper_rectangle( + static_cast &>(tr), sub, Point<2>(0, 0), Point<2>(1, 1)); + + FESystem fe(FE_Q(1), 2, FE_DGQ(0), 1); + DoFHandler dofh(tr); + dofh.distribute_dofs(fe); + dofh.distribute_mg_dofs(); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + deallog << "Total dofs=" << dofh.n_dofs() << std::endl; + + // extract constant modes and print + if (myid == 0) + for (unsigned int c = 0; c < fe.n_components(); ++c) + { + ComponentMask mask(fe.n_components(), false); + mask.set(c, true); + + std::vector> constant_modes; + DoFTools::extract_level_constant_modes(0, dofh, mask, constant_modes); + + for (unsigned int i = 0; i < constant_modes.size(); ++i) + { + for (unsigned int j = 0; j < constant_modes[i].size(); ++j) + deallog << (constant_modes[i][j] ? '1' : '0') << ' '; + deallog << std::endl; + } + } + + // renumber dofs and do the same again + DoFRenumbering::component_wise(dofh, 0); + if (myid == 0) + for (unsigned int c = 0; c < fe.n_components(); ++c) + { + ComponentMask mask(fe.n_components(), false); + mask.set(c, true); + + std::vector> constant_modes; + DoFTools::extract_level_constant_modes(0, dofh, mask, constant_modes); + + for (unsigned int i = 0; i < constant_modes.size(); ++i) + { + for (unsigned int j = 0; j < constant_modes[i].size(); ++j) + deallog << (constant_modes[i][j] ? '1' : '0') << ' '; + deallog << std::endl; + } + } +} + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); + } + else + { + deallog.push("2d"); + test<2>(); + deallog.pop(); + } +} diff --git a/tests/mpi/extract_constant_modes_04.mpirun=10.with_p4est=true.output b/tests/mpi/extract_constant_modes_04.mpirun=10.with_p4est=true.output new file mode 100644 index 0000000000..ad360700d2 --- /dev/null +++ b/tests/mpi/extract_constant_modes_04.mpirun=10.with_p4est=true.output @@ -0,0 +1,8 @@ + +DEAL:0:2d::Total dofs=54 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 diff --git a/tests/mpi/extract_constant_modes_04.mpirun=4.with_p4est=true.output b/tests/mpi/extract_constant_modes_04.mpirun=4.with_p4est=true.output new file mode 100644 index 0000000000..0378ef06a2 --- /dev/null +++ b/tests/mpi/extract_constant_modes_04.mpirun=4.with_p4est=true.output @@ -0,0 +1,8 @@ + +DEAL:0:2d::Total dofs=24 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 1 1 1 +DEAL:0:2d::1 -- 2.39.5