From 92fc2a92f1659eeb53e994bab9409b26e7185d61 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 10 Mar 2022 08:02:58 -0500 Subject: [PATCH] Implement DoFRenumbering::support_point_wise(). --- doc/news/changes/major/20220312DavidWells | 5 + include/deal.II/dofs/dof_renumbering.h | 49 ++++ source/dofs/dof_renumbering.cc | 154 +++++++++++ source/dofs/dof_renumbering.inst.in | 4 + tests/dofs/nodal_renumbering_01.cc | 241 ++++++++++++++++++ ...mbering_01.mpirun=2.with_p4est=true.output | 31 +++ ...mbering_01.mpirun=3.with_p4est=true.output | 47 ++++ ...odal_renumbering_01.with_p4est=true.output | 15 ++ 8 files changed, 546 insertions(+) create mode 100644 doc/news/changes/major/20220312DavidWells create mode 100644 tests/dofs/nodal_renumbering_01.cc create mode 100644 tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.output create mode 100644 tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.output create mode 100644 tests/dofs/nodal_renumbering_01.with_p4est=true.output diff --git a/doc/news/changes/major/20220312DavidWells b/doc/news/changes/major/20220312DavidWells new file mode 100644 index 0000000000..786b9d5e36 --- /dev/null +++ b/doc/news/changes/major/20220312DavidWells @@ -0,0 +1,5 @@ +New: Added a new function DoFRenumbering::support_point_wise() +which renumbers DoFs so that DoFs with identical support points +are now numbered adjacently. +
+(David Wells, 2021/03/11) diff --git a/include/deal.II/dofs/dof_renumbering.h b/include/deal.II/dofs/dof_renumbering.h index 82fb88919f..6aa9a9a5a6 100644 --- a/include/deal.II/dofs/dof_renumbering.h +++ b/include/deal.II/dofs/dof_renumbering.h @@ -1214,6 +1214,55 @@ namespace DoFRenumbering * @} */ + /** + * @name Numberings based on properties of the finite element space + * @{ + */ + + /** + * Stably sort DoFs first by component and second by location of their support + * points, so that DoFs with the same support point are listed consecutively + * in the component order. + * + * The primary use of this ordering is that it enables one to interpret a + * vector of FE coefficients as a vector of tensors. For example, suppose `X` + * is a vector containing coordinates (i.e., the sort of vector one would use + * with MappingFEField) and `U` is a vector containing velocities in 2D. Then + * the `k`th support point is mapped to `{X[2*k], X[2*k + 1]}` and the + * velocity there is `{U[2*k], U[2*k + 1]}`. Hence, with this reordering, one + * can read solution data at each support point without additional indexing. + * This is useful for, e.g., passing vectors of FE coefficients to external + * libraries which expect nodal data in this format. + * + * @warning This function only supports finite elements which have the same + * number of DoFs in each component. This is checked with an assertion. + * + * @note This renumbering assumes that the base elements of each vector-valued + * element in @p dof_handler are numbered in the way FESystem currently + * distributes DoFs: i.e., for a given support point, the global dof indices + * for component i should be less than the global dof indices for component i + * + 1. Due to various technical complications (like DoF unification in + * hp-mode) this assumption needs to be satified for this function to work in + * all relevant cases. + */ + template + void + support_point_wise(DoFHandler &dof_handler); + + /** + * Compute the renumbering vector needed by the support_point_wise() function. + * Does not perform the renumbering on the @p DoFHandler dofs but returns the + * renumbering vector. + */ + template + void + compute_support_point_wise( + std::vector &new_dof_indices, + const DoFHandler & dof_handler); + + /** + * @} + */ /** diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc index d4c31794ff..debfb8b4f6 100644 --- a/source/dofs/dof_renumbering.cc +++ b/source/dofs/dof_renumbering.cc @@ -53,6 +53,7 @@ DEAL_II_DISABLE_EXTRA_DIAGNOSTICS #include #include #include + #undef BOOST_BIND_GLOBAL_PLACEHOLDERS DEAL_II_ENABLE_EXTRA_DIAGNOSTICS @@ -2248,6 +2249,159 @@ namespace DoFRenumbering ExcInternalError()); } + + + template + void + support_point_wise(DoFHandler &dof_handler) + { + std::vector renumbering( + dof_handler.n_locally_owned_dofs(), numbers::invalid_dof_index); + compute_support_point_wise(renumbering, dof_handler); + + // If there is only one component then there is nothing to do, so check + // first: + if (Utilities::MPI::max(renumbering.size(), + dof_handler.get_communicator()) > 0) + dof_handler.renumber_dofs(renumbering); + } + + + + template + void + compute_support_point_wise( + std::vector &new_dof_indices, + const DoFHandler & dof_handler) + { + const types::global_dof_index n_dofs = dof_handler.n_locally_owned_dofs(); + Assert(new_dof_indices.size() == n_dofs, + ExcDimensionMismatch(new_dof_indices.size(), n_dofs)); + + // This renumbering occurs in three steps: + // 1. Compute the component-wise renumbering so that all DoFs of component + // i are less than the DoFs of component i + 1. + // 2. Compute a second renumbering component_to_nodal in which the + // renumbering is now, for two components, [u0, v0, u1, v1, ...]: i.e., + // DoFs are first sorted by component and then by support point. + // 3. Compose the two renumberings to obtain the final result. + + // Step 1: + std::vector component_renumbering( + n_dofs, numbers::invalid_dof_index); + compute_component_wise(component_renumbering, + dof_handler.begin_active(), + dof_handler.end(), + std::vector(), + false); + + const std::vector dofs_per_component = + DoFTools::count_dofs_per_fe_component(dof_handler, true); + + // At this point we have no more communication to do - simplify things by + // returning early if possible + if (component_renumbering.size() == 0) + { + new_dof_indices.resize(0); + return; + } + std::fill(new_dof_indices.begin(), + new_dof_indices.end(), + numbers::invalid_dof_index); + // This index set equals what dof_handler.locally_owned_dofs() would be if + // we executed the componentwise renumbering. + IndexSet component_renumbered_dofs(dof_handler.n_dofs()); + component_renumbered_dofs.add_indices(component_renumbering.begin(), + component_renumbering.end()); + for (const auto &dpc : dofs_per_component) + AssertThrow(dofs_per_component[0] == dpc, ExcNotImplemented()); + for (const FiniteElement &fe : + dof_handler.get_fe_collection()) + { + AssertThrow(fe.dofs_per_cell == 0 || fe.has_support_points(), + ExcNotImplemented()); + for (unsigned int i = 0; i < fe.n_base_elements(); ++i) + AssertThrow( + fe.base_element(0).get_unit_support_points() == + fe.base_element(i).get_unit_support_points(), + ExcMessage( + "All base elements should have the same support points.")); + } + + std::vector component_to_nodal( + dof_handler.n_locally_owned_dofs(), numbers::invalid_dof_index); + + // Step 2: + std::vector cell_dofs; + std::vector component_renumbered_cell_dofs; + const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); + // Reuse the original index space for the new renumbering: it is the right + // size and is contiguous on the current processor + auto next_dof_it = locally_owned_dofs.begin(); + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + const FiniteElement &fe = cell->get_fe(); + cell_dofs.resize(fe.dofs_per_cell); + component_renumbered_cell_dofs.resize(fe.dofs_per_cell); + cell->get_dof_indices(cell_dofs); + // Apply the component renumbering while skipping any ghost dofs. This + // algorithm assumes that all locally owned DoFs before the component + // renumbering are still locally owned afterwards (just with a new + // index). + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + { + if (locally_owned_dofs.is_element(cell_dofs[i])) + { + const auto local_index = + locally_owned_dofs.index_within_set(cell_dofs[i]); + component_renumbered_cell_dofs[i] = + component_renumbering[local_index]; + } + else + { + component_renumbered_cell_dofs[i] = + numbers::invalid_dof_index; + } + } + + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + { + if (fe.system_to_component_index(i).first == 0 && + component_renumbered_dofs.is_element( + component_renumbered_cell_dofs[i])) + { + for (unsigned int component = 0; + component < fe.n_components(); + ++component) + { + // Since we are indexing in an odd way here it is much + // simpler to compute the permutation separately and + // combine it at the end instead of doing both at once + const auto local_index = + component_renumbered_dofs.index_within_set( + component_renumbered_cell_dofs[i] + + dofs_per_component[0] * component); + + if (component_to_nodal[local_index] == + numbers::invalid_dof_index) + { + component_to_nodal[local_index] = *next_dof_it; + ++next_dof_it; + } + } + } + } + } + + // Step 3: + for (std::size_t i = 0; i < dof_handler.n_locally_owned_dofs(); ++i) + { + const auto local_index = + component_renumbered_dofs.index_within_set(component_renumbering[i]); + new_dof_indices[i] = component_to_nodal[local_index]; + } + } } // namespace DoFRenumbering diff --git a/source/dofs/dof_renumbering.inst.in b/source/dofs/dof_renumbering.inst.in index 75afeb9820..3f97733f7e 100644 --- a/source/dofs/dof_renumbering.inst.in +++ b/source/dofs/dof_renumbering.inst.in @@ -64,6 +64,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) template void hierarchical(DoFHandler &); + template void + support_point_wise( + DoFHandler &); + \} #endif } diff --git a/tests/dofs/nodal_renumbering_01.cc b/tests/dofs/nodal_renumbering_01.cc new file mode 100644 index 0000000000..4c767bfbd0 --- /dev/null +++ b/tests/dofs/nodal_renumbering_01.cc @@ -0,0 +1,241 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 - 2022 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +#include "../simplex/simplex_grids.h" + +// Test DoFRenumbering::support_point_wise(). The vector position should contain +// (x0, y0, z0, x1, y1, z1, ...) and solution1 should contain (u0, v0, w0, u1, +// v1, w1, ...). Verify that they match when we interpolate. + +using namespace dealii; + +class Test : public Function<2> +{ +public: + Test(const unsigned int n_components) + : Function<2>(n_components) + {} + + double + value(const Point<2> &p, const unsigned int component = 0) const override + { + switch (component) + { + case 0: + return std::sin(p[0]) * std::sin(2.0 * p[1]) + 1.0; + case 1: + return std::cos(p[0]) * std::cos(3.0 * p[1]) + 2.0; + default: + Assert(false, ExcNotImplemented()); + } + return 0.0; + } +}; + +void +test(DoFHandler<2> &dof_handler, const hp::MappingCollection<2> &mappings) +{ + DoFRenumbering::support_point_wise(dof_handler); + const MPI_Comm comm = dof_handler.get_communicator(); + + const IndexSet &local_dofs = dof_handler.locally_owned_dofs(); + deallog << "locally owned dofs =\n"; + local_dofs.print(deallog); + deallog << std::endl; + const IndexSet relevant_dofs = + DoFTools::extract_locally_relevant_dofs(dof_handler); + LinearAlgebra::distributed::Vector position(local_dofs, + relevant_dofs, + comm); + LinearAlgebra::distributed::Vector solution1(local_dofs, + relevant_dofs, + comm); + LinearAlgebra::distributed::Vector solution2(local_dofs, + relevant_dofs, + comm); + + const unsigned int n_components = + dof_handler.get_fe_collection().n_components(); + // It doesn't make sense to treat position as a vector of coordinates unless + // we have enough components + if (n_components == 2) + { + VectorTools::interpolate(mappings, + dof_handler, + Functions::IdentityFunction<2>(), + position); + + Test test(n_components); + if (n_components == 2) + VectorTools::interpolate(mappings, dof_handler, test, solution1); + + // Verify that we get the same results when we interpolate either manually + // or by reading off position data and evaluating the interpolated + // function + for (unsigned int node_n = 0; + node_n < position.locally_owned_size() / n_components; + ++node_n) + { + const auto i0 = node_n * n_components; + const auto i1 = node_n * n_components + 1; + Point<2> p(position.local_element(i0), position.local_element(i1)); + solution2.local_element(i0) = test.value(p, 0); + solution2.local_element(i1) = test.value(p, 1); + } + LinearAlgebra::distributed::Vector difference = solution1; + difference -= solution2; + deallog << "difference norm = " << difference.l2_norm() << std::endl; + } + +#if 0 + DataOut<2> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(position, "X"); + data_out.add_data_vector(solution1, "U1"); + data_out.add_data_vector(solution2, "U2"); + data_out.build_patches(4); + data_out.write_vtu_with_pvtu_record( + "./", "solution", 0, comm, 2, 8); +#endif +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + MPILogInitAll all; + const MPI_Comm comm = MPI_COMM_WORLD; + + // Test with p::s::T, mixed FE, multiple components + { + parallel::shared::Triangulation<2> tria( + comm, + dealii::Triangulation<2>::none, + true, + parallel::shared::Triangulation<2>::partition_zorder); + GridGenerator::cube_and_pyramid(tria); + + hp::FECollection<2> fe(FESystem<2>(FE_Q<2>(1), 2), + FESystem<2>(FE_SimplexP<2>(1), 2)); + hp::MappingCollection<2> mappings(MappingQ<2>(1), + MappingFE<2>(FE_SimplexP<2>(1))); + DoFHandler<2> dof_handler(tria); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + if (cell->reference_cell() == ReferenceCells::Quadrilateral) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); + } + } + dof_handler.distribute_dofs(fe); + + test(dof_handler, mappings); + } + + // Test with a scalar FE + { + parallel::shared::Triangulation<2> tria( + comm, + dealii::Triangulation<2>::none, + true, + parallel::shared::Triangulation<2>::partition_zorder); + GridGenerator::cube_and_pyramid(tria); + + hp::FECollection<2> fe(FE_Q<2>(1), FE_SimplexP<2>(1)); + hp::MappingCollection<2> mappings(MappingQ<2>(1), + MappingFE<2>(FE_SimplexP<2>(1))); + DoFHandler<2> dof_handler(tria); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + if (cell->reference_cell() == ReferenceCells::Quadrilateral) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); + } + } + dof_handler.distribute_dofs(fe); + + test(dof_handler, mappings); + } + + // Test with another scalar FE + { + parallel::shared::Triangulation<2> tria( + comm, + dealii::Triangulation<2>::none, + true, + parallel::shared::Triangulation<2>::partition_zorder); + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FE_Q<2> fe(5); + hp::MappingCollection<2> mappings(MappingQ<2>(1)); + DoFHandler<2> dof_handler(tria); + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + cell->set_active_fe_index(0); + dof_handler.distribute_dofs(fe); + + test(dof_handler, mappings); + } + + // Test with another vector FE + grid refinement + { + parallel::distributed::Triangulation<2> tria(comm); + GridGenerator::hyper_cube(tria); + tria.refine_global(3); + + FESystem<2> fe(FE_Q<2>(5), 2); + DoFHandler<2> dof_handler(tria); + dof_handler.distribute_dofs(fe); + hp::MappingCollection<2> mappings(MappingQ<2>(1)); + for (const auto &cell : dof_handler.active_cell_iterators()) + cell->set_active_fe_index(0); + + test(dof_handler, mappings); + } +} diff --git a/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.output b/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.output new file mode 100644 index 0000000000..7dadd62145 --- /dev/null +++ b/tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.output @@ -0,0 +1,31 @@ + +DEAL:0::locally owned dofs = +{[0,7]} +DEAL:0:: +DEAL:0::difference norm = 0.00000 +DEAL:0::locally owned dofs = +{[0,3]} +DEAL:0:: +DEAL:0::locally owned dofs = +{[0,230]} +DEAL:0:: +DEAL:0::locally owned dofs = +{[0,1721]} +DEAL:0:: +DEAL:0::difference norm = 0.00000 + +DEAL:1::locally owned dofs = +{[8,9]} +DEAL:1:: +DEAL:1::difference norm = 0.00000 +DEAL:1::locally owned dofs = +{4} +DEAL:1:: +DEAL:1::locally owned dofs = +{[231,440]} +DEAL:1:: +DEAL:1::locally owned dofs = +{[1722,3361]} +DEAL:1:: +DEAL:1::difference norm = 0.00000 + diff --git a/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.output b/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.output new file mode 100644 index 0000000000..bd56e24115 --- /dev/null +++ b/tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.output @@ -0,0 +1,47 @@ + +DEAL:0::locally owned dofs = +{} +DEAL:0:: +DEAL:0::difference norm = 0.00000 +DEAL:0::locally owned dofs = +{} +DEAL:0:: +DEAL:0::locally owned dofs = +{[0,120]} +DEAL:0:: +DEAL:0::locally owned dofs = +{[0,1101]} +DEAL:0:: +DEAL:0::difference norm = 0.00000 + +DEAL:1::locally owned dofs = +{[0,7]} +DEAL:1:: +DEAL:1::difference norm = 0.00000 +DEAL:1::locally owned dofs = +{[0,3]} +DEAL:1:: +DEAL:1::locally owned dofs = +{[121,340]} +DEAL:1:: +DEAL:1::locally owned dofs = +{[1102,2361]} +DEAL:1:: +DEAL:1::difference norm = 0.00000 + + +DEAL:2::locally owned dofs = +{[8,9]} +DEAL:2:: +DEAL:2::difference norm = 0.00000 +DEAL:2::locally owned dofs = +{4} +DEAL:2:: +DEAL:2::locally owned dofs = +{[341,440]} +DEAL:2:: +DEAL:2::locally owned dofs = +{[2362,3361]} +DEAL:2:: +DEAL:2::difference norm = 0.00000 + diff --git a/tests/dofs/nodal_renumbering_01.with_p4est=true.output b/tests/dofs/nodal_renumbering_01.with_p4est=true.output new file mode 100644 index 0000000000..3826b2a5bd --- /dev/null +++ b/tests/dofs/nodal_renumbering_01.with_p4est=true.output @@ -0,0 +1,15 @@ + +DEAL:0::locally owned dofs = +{[0,9]} +DEAL:0:: +DEAL:0::difference norm = 0.00000 +DEAL:0::locally owned dofs = +{[0,4]} +DEAL:0:: +DEAL:0::locally owned dofs = +{[0,440]} +DEAL:0:: +DEAL:0::locally owned dofs = +{[0,3361]} +DEAL:0:: +DEAL:0::difference norm = 0.00000 -- 2.39.5