]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement DoFRenumbering::support_point_wise().
authorDavid Wells <drwells@email.unc.edu>
Thu, 10 Mar 2022 13:02:58 +0000 (08:02 -0500)
committerDavid Wells <drwells@email.unc.edu>
Fri, 25 Mar 2022 22:36:10 +0000 (18:36 -0400)
doc/news/changes/major/20220312DavidWells [new file with mode: 0644]
include/deal.II/dofs/dof_renumbering.h
source/dofs/dof_renumbering.cc
source/dofs/dof_renumbering.inst.in
tests/dofs/nodal_renumbering_01.cc [new file with mode: 0644]
tests/dofs/nodal_renumbering_01.mpirun=2.with_p4est=true.output [new file with mode: 0644]
tests/dofs/nodal_renumbering_01.mpirun=3.with_p4est=true.output [new file with mode: 0644]
tests/dofs/nodal_renumbering_01.with_p4est=true.output [new file with mode: 0644]

diff --git a/doc/news/changes/major/20220312DavidWells b/doc/news/changes/major/20220312DavidWells
new file mode 100644 (file)
index 0000000..786b9d5
--- /dev/null
@@ -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.
+<br>
+(David Wells, 2021/03/11)
index 82fb88919f6762ddc835436f4143eca2e08964be..6aa9a9a5a6b7d28ee85bcc74dc243468c8f1d40e 100644 (file)
@@ -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 <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);
+
+  /**
+   * @}
+   */
 
 
   /**
index d4c31794ff9c7e98788ec76695161bcf8180f26c..debfb8b4f67a876d5eee6f8b12ed5abb9d6eac6d 100644 (file)
@@ -53,6 +53,7 @@ DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
 #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
 
@@ -2248,6 +2249,159 @@ namespace DoFRenumbering
            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
 
 
index 75afeb9820d643754eb4b1c41adbb49c0366cd3c..3f97733f7e5c5ed27990722cee71b55b7c72baed 100644 (file)
@@ -64,6 +64,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
       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
   }
diff --git a/tests/dofs/nodal_renumbering_01.cc b/tests/dofs/nodal_renumbering_01.cc
new file mode 100644 (file)
index 0000000..4c767bf
--- /dev/null
@@ -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 <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);
+  }
+}
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 (file)
index 0000000..7dadd62
--- /dev/null
@@ -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 (file)
index 0000000..bd56e24
--- /dev/null
@@ -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 (file)
index 0000000..3826b2a
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.