--- /dev/null
+New: Added a new function DoFRenumbering::support_point_wise()
+which renumbers DoFs so that DoFs with identical support points
+are now numbered adjacently.
+<br>
+(David Wells, 2021/03/11)
* @}
*/
+ /**
+ * @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 <int dim, int spacedim>
+ void
+ support_point_wise(DoFHandler<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ compute_support_point_wise(
+ std::vector<types::global_dof_index> &new_dof_indices,
+ const DoFHandler<dim, spacedim> & dof_handler);
+
+ /**
+ * @}
+ */
/**
#include <boost/graph/properties.hpp>
#include <boost/random.hpp>
#include <boost/random/uniform_int_distribution.hpp>
+
#undef BOOST_BIND_GLOBAL_PLACEHOLDERS
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
ExcInternalError());
}
+
+
+ template <int dim, int spacedim>
+ void
+ support_point_wise(DoFHandler<dim, spacedim> &dof_handler)
+ {
+ std::vector<types::global_dof_index> 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 <int dim, int spacedim>
+ void
+ compute_support_point_wise(
+ std::vector<types::global_dof_index> &new_dof_indices,
+ const DoFHandler<dim, spacedim> & 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<types::global_dof_index> component_renumbering(
+ n_dofs, numbers::invalid_dof_index);
+ compute_component_wise<dim, spacedim>(component_renumbering,
+ dof_handler.begin_active(),
+ dof_handler.end(),
+ std::vector<unsigned int>(),
+ false);
+
+ const std::vector<types::global_dof_index> 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<dim, spacedim> &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<types::global_dof_index> component_to_nodal(
+ dof_handler.n_locally_owned_dofs(), numbers::invalid_dof_index);
+
+ // Step 2:
+ std::vector<types::global_dof_index> cell_dofs;
+ std::vector<types::global_dof_index> 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<dim, spacedim> &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
template void
hierarchical(DoFHandler<deal_II_dimension, deal_II_space_dimension> &);
+ template void
+ support_point_wise(
+ DoFHandler<deal_II_dimension, deal_II_space_dimension> &);
+
\}
#endif
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/function_lib.h>
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools_interpolate.h>
+
+#include <boost/algorithm/apply_permutation.hpp>
+
+#include <iostream>
+#include <map>
+
+#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<double> position(local_dofs,
+ relevant_dofs,
+ comm);
+ LinearAlgebra::distributed::Vector<double> solution1(local_dofs,
+ relevant_dofs,
+ comm);
+ LinearAlgebra::distributed::Vector<double> 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<double> 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);
+ }
+}
--- /dev/null
+
+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
+
--- /dev/null
+
+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
+
--- /dev/null
+
+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
projects / dealii.git / commitdiff