]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Renumber lexycographic + face patch dof numbering 18413/head
authorMichał Wichrowski <mtwichrowski@gmail.com>
Fri, 4 Jul 2025 13:35:25 +0000 (15:35 +0200)
committerMichał Wichrowski <mtwichrowski@gmail.com>
Fri, 4 Jul 2025 13:35:25 +0000 (15:35 +0200)
12 files changed:
doc/news/changes/minor/20250502Wichrowski [new file with mode: 0644]
doc/news/changes/minor/20250701Wichrowski [new file with mode: 0644]
include/deal.II/dofs/dof_renumbering.h
include/deal.II/fe/fe_tools.h
include/deal.II/fe/fe_tools.templates.h
source/dofs/dof_renumbering.cc
source/dofs/dof_renumbering.inst.in
source/fe/fe_tools.inst.in
tests/dofs/dof_renumbering_11.cc [new file with mode: 0644]
tests/dofs/dof_renumbering_11.output [new file with mode: 0644]
tests/fe/cell_face_lex.cc [new file with mode: 0644]
tests/fe/cell_face_lex.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20250502Wichrowski b/doc/news/changes/minor/20250502Wichrowski
new file mode 100644 (file)
index 0000000..796bc13
--- /dev/null
@@ -0,0 +1,4 @@
+New: A new function FETools::cell_to_face_lexicographic() to generate a
+mapping from cell-local DoFs to lexicographic ordering of DoFs on two adjacent cells
+has been added.
+(Michał Wichrowski, 2025/05/03)
diff --git a/doc/news/changes/minor/20250701Wichrowski b/doc/news/changes/minor/20250701Wichrowski
new file mode 100644 (file)
index 0000000..0f11fec
--- /dev/null
@@ -0,0 +1,2 @@
+New: The ordering strategy DoFRenumbering::lexicographic() has been added.
+(Michał Wichrowski, 2025/07/01)
index a2c5a9ca0bef824b52b01536aa78d529c7200480..39a9a36dfaa1b66427b53a49f3ba3d0547a3e124 100644 (file)
@@ -1273,6 +1273,38 @@ namespace DoFRenumbering
     std::vector<types::global_dof_index> &new_dof_indices,
     const DoFHandler<dim, spacedim>      &dof_handler);
 
+  /**
+   * Reorder the degrees of freedom in lexicographic order of support points.
+   * The tolerance is used to determine whether two support points are equal.
+   *
+   * This is particularly useful for finite element spaces that are defined on
+   * tensor product grids, such as the ones generated by
+   * GridGenerator::hyper_rectangle() or GridGenerator::hyper_cube() or
+   * their subdivided variants.
+   * With this ordering certain higher-dimension finite element matrices can be
+   * expressed as Kronecker products of lower-dimensional matrices.
+   *
+   * @warning parallel::distributed::Triangulation objects are not supported.
+   * The finite element must have support points.
+   */
+  template <int dim>
+  void
+  lexicographic(DoFHandler<dim> &dof_handler, const double tolerance = 1e-12);
+
+  /**
+   * Compute the renumbering vector needed by the lexicographic() function.
+   * Does not perform the renumbering on the @p DoFHandler dofs but returns the
+   * renumbering vector.
+   * The tolerance is treated as an absolute tolerance for comparing support
+   * point coordinates. To use a relative tolerance, multiply a relative value
+   * by dof_handler.begin(0)->diameter() or a similar characteristic length.
+   */
+  template <int dim>
+  void
+  compute_lexicographic(std::vector<types::global_dof_index> &new_dof_indices,
+                        const DoFHandler<dim>                &handler,
+                        const double tolerance = 1e-12);
+
   /**
    * @}
    */
index 69ed37fe1c4b5a6bbaaa458efea42f5c3aae2a39..c34cec1ab62f5b0ee626d19ec8dcdaf1d5d4aca7 100644 (file)
@@ -931,6 +931,34 @@ namespace FETools
   std::vector<unsigned int>
   lexicographic_to_hierarchic_numbering(unsigned int degree);
 
+  /**
+   * Given the polynomial degree, direction, and numbering options, this
+   * function returns a pair of vectors mapping cell DoFs to face patch DoFs
+   * in lexicographic order.
+   *
+   * This function works for scalar finite elements, specifically those derived
+   * from FE_Q and FE_DGQ. It computes the mapping from the DoFs of
+   * two adjacent cells sharing a face perpendicular to @p direction in
+   * reference space, to the DoFs on the face patch. The result is a pair of
+   * vectors: the first for the DoFs on the first cell, the second for the
+   * DoFs on the second cell. The @p cell_hierarchical_numbering flag
+   * determines whether hierarchical or lexicographic numbering is assumed for
+   * the cell DoFs. The @p is_continuous flag indicates whether the finite
+   * element space is continuous (i.e., FE_Q).
+   *
+   * @param degree Polynomial degree of the finite element.
+   * @param direction Axis perpendicular to the considered face.
+   * @param cell_hierarchical_numbering If true, use hierarchical numbering for cell DoFs.
+   * @param is_continuous If true, assumes the finite element space is continuous (FE_Q).
+   * @return A pair of vectors mapping cell DoFs to face patch DoFs in lexicographic order.
+   */
+  template <int dim>
+  std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
+  cell_to_face_patch(const unsigned int &degree,
+                     const unsigned int &direction,
+                     const bool         &cell_hierarchical_numbering,
+                     const bool         &is_continuous);
+
   /**
    * A namespace that contains functions that help setting up internal
    * data structures when implementing FiniteElement which are build
index 420be5f268997d908ca096d42f2b617f500b7cca..65167f6c5329294dab45303680c51a8a9ddc1699 100644 (file)
@@ -3243,6 +3243,104 @@ namespace FETools
     return Utilities::invert_permutation(
       hierarchic_to_lexicographic_numbering<dim>(degree));
   }
+
+
+
+  template <int dim>
+  std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
+  cell_to_face_patch(const unsigned int &degree,
+                     const unsigned int &direction,
+                     const bool         &cell_hierarchical_numbering,
+                     const bool         &is_continuous)
+  {
+    AssertThrow(dim > 1, ExcInvalidFEDimension(1, dim));
+    AssertThrow(direction < dim, ExcInvalidFEDimension(direction, dim));
+
+    // n_1d: DoFs per 1D edge (parallel to face)
+    const unsigned int n_1d = degree + 1;
+    // n_ort: DoFs in the direction orthogonal to the face in the combined face
+    // space For discontinuous: 2*n (n for cell0 face, n for cell1 face) For
+    // continuous: 2*n - 1 (shared DoFs)
+    const unsigned int n_ort_1d = 2 * n_1d - (is_continuous ? 1 : 0);
+
+    const unsigned int dofs_per_cell = Utilities::fixed_power<dim>(n_1d);
+    const unsigned int dofs_per_face_patch =
+      Utilities::fixed_power<dim - 1>(n_1d) * n_ort_1d;
+
+    std::array<unsigned int, dim> face_patch_dofs_in_direction;
+    for (unsigned int d = 0; d < dim; ++d)
+      face_patch_dofs_in_direction[d] = n_1d;
+    face_patch_dofs_in_direction[direction] = n_ort_1d;
+
+    // Initialize maps for cell0 and cell1 DoFs to face DoFs
+    std::array<std::vector<unsigned int>, 2> results;
+    results[0].resize(dofs_per_cell, numbers::invalid_unsigned_int);
+    results[1].resize(dofs_per_cell, numbers::invalid_unsigned_int);
+
+    // Compute with lexicographic first
+    for (unsigned int i = 0; i < dofs_per_face_patch; ++i)
+      {
+        // Decompose i into multiindex
+        std::array<unsigned int, dim> face_multiindex;
+        unsigned int                  temp = i;
+        for (unsigned int d = 0; d < dim; ++d)
+          {
+            face_multiindex[d] = temp % face_patch_dofs_in_direction[d];
+            temp /= face_patch_dofs_in_direction[d];
+          }
+
+        std::array<unsigned int, dim> cell_multiindex = face_multiindex;
+
+        // compute corresponding cell
+        unsigned int cell          = face_multiindex[direction] / n_1d;
+        cell_multiindex[direction] = face_multiindex[direction] % n_1d;
+        if (is_continuous)
+          cell_multiindex[direction] += cell;
+
+        {
+          unsigned int cell_index = 0;
+          unsigned int stride     = 1;
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              cell_index += cell_multiindex[d] * stride;
+              stride *= n_1d;
+            }
+          results[cell][cell_index] = i;
+        }
+
+        if (face_multiindex[direction] == n_1d - 1 && is_continuous)
+          {
+            // fill the other cell
+            unsigned other_cell        = 1;
+            cell_multiindex[direction] = 0;
+            unsigned int cell_index    = 0;
+            unsigned int stride        = 1;
+            for (unsigned int d = 0; d < dim; ++d)
+              {
+                cell_index += cell_multiindex[d] * stride;
+                stride *= n_1d;
+              }
+            results[other_cell][cell_index] = i;
+          }
+      }
+    // Now we have the DoFs in lexicographic order
+    if (!cell_hierarchical_numbering)
+      return std::make_pair(results[0], results[1]);
+
+    // Renumbering required for hierarchical numbering
+    std::vector<unsigned int> lex_to_hie =
+      lexicographic_to_hierarchic_numbering<dim>(degree);
+
+    for (auto &result : results)
+      {
+        std::vector<unsigned int> renumbered_result = result;
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          renumbered_result[lex_to_hie[i]] = result[i];
+        result = renumbered_result;
+      }
+    return std::make_pair(results[0], results[1]);
+  }
+
 } // namespace FETools
 
 
index a0b0f3240e90b1f1e2f2611b9755d710ace1a117..8600dbf6d5e2ea7365c1e688291b2fa4b587d239 100644 (file)
@@ -2442,6 +2442,72 @@ namespace DoFRenumbering
 
 
 
+  template <int dim>
+  void
+  lexicographic(DoFHandler<dim> &dof_handler, const double tolerance)
+  {
+    std::vector<types::global_dof_index> renumbering;
+    compute_lexicographic(renumbering, dof_handler, tolerance);
+    dof_handler.renumber_dofs(renumbering);
+  }
+
+
+  template <int dim>
+  void
+  compute_lexicographic(std::vector<types::global_dof_index> &new_dof_indices,
+                        const DoFHandler<dim>                &dof_handler,
+                        const double                          tolerance)
+  {
+    Assert(dynamic_cast<const parallel::distributed::Triangulation<dim> *>(
+             &dof_handler.get_triangulation()) == nullptr,
+           ExcMessage(
+             "Lexicographic renumbering is not implemented for distributed "
+             "triangulations."));
+
+    std::map<types::global_dof_index, Point<dim>> dof_location_map;
+    DoFTools::map_dofs_to_support_points(MappingQ1<dim>(),
+                                         dof_handler,
+                                         dof_location_map);
+    std::vector<std::pair<types::global_dof_index, Point<dim>>>
+      dof_location_vector;
+
+    dof_location_vector.reserve(dof_location_map.size());
+    for (const auto &s : dof_location_map)
+      dof_location_vector.push_back(s);
+
+    std::sort(
+      dof_location_vector.begin(),
+      dof_location_vector.end(),
+      [&](const std::pair<types::global_dof_index, Point<dim>> &p1,
+          const std::pair<types::global_dof_index, Point<dim>> &p2) -> bool {
+        for (int i = dim - 1; i >= 0; --i)
+          {
+            const double diff = p1.second(i) - p2.second(i);
+            // Check if p1 is significantly smaller than p2 in this dimension
+            if (diff < -tolerance)
+              return true;
+            // Check if p1 is significantly larger than p2 in this dimension
+            if (diff > tolerance)
+              return false;
+            // Otherwise, coordinates are considered equal within tolerance for
+            // this dimension, continue to the next dimension.
+          }
+        // All coordinates are within tolerance, points are considered
+        // equivalent. Use the original DoF index as a tie-breaker to preserve
+        // relative order.
+        return p1.first < p2.first;
+      });
+
+    new_dof_indices.resize(dof_location_vector.size(),
+                           numbers::invalid_dof_index);
+
+    for (types::global_dof_index dof = 0; dof < dof_location_vector.size();
+         ++dof)
+      new_dof_indices[dof_location_vector[dof].first] = dof;
+  }
+
+
+
   template <int dim,
             int spacedim,
             typename Number,
index 7d141aa34345daf442066a506917d37a344874f1..f3a0f0bd1adb478bc877e9c3f99717ae3799415a 100644 (file)
@@ -258,6 +258,9 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         const unsigned int,
         const bool,
         const std::vector<types::global_dof_index> &);
+
+      template void
+      lexicographic(DoFHandler<deal_II_dimension> &, const double);
     \}
 #endif
   }
index 7c9e101c04ee93aeaed06e244c6e1386179917ee..7dfec11c67e99ea81108142bc2aa0a93c2e3eb13 100644 (file)
@@ -262,6 +262,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
       template std::vector<unsigned int>
       lexicographic_to_hierarchic_numbering<deal_II_dimension>(unsigned int);
 
+      template std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
+      cell_to_face_patch<deal_II_dimension>(const unsigned int &,
+                                            const unsigned int &,
+                                            const bool &,
+                                            const bool &);
+
 #endif
     \}
   }
diff --git a/tests/dofs/dof_renumbering_11.cc b/tests/dofs/dof_renumbering_11.cc
new file mode 100644 (file)
index 0000000..bbf3c9e
--- /dev/null
@@ -0,0 +1,93 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2019 - 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 lexicographic renumbering on a patch
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_accessor.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_system.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test(const unsigned int fe_degree, unsigned int n_components)
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+  tria.refine_global(1);
+
+  deallog << "Testing dim=" << dim << ", degree=" << fe_degree
+          << ", components=" << n_components << std::endl;
+
+  FESystem<dim>   fe(FE_Q<dim>(fe_degree), n_components);
+  DoFHandler<dim> dof_handler(tria);
+
+  const unsigned int n_dof_1d = (2 * fe_degree + 1);
+
+  dof_handler.distribute_dofs(fe);
+  DoFRenumbering::lexicographic(dof_handler);
+  DoFRenumbering::component_wise(dof_handler);
+
+  std::vector<types::global_dof_index> local_dof_indices(fe.dofs_per_cell);
+
+  std::map<types::global_dof_index, Point<dim>> support_points;
+  MappingQ1<dim>                                mapping;
+  DoFTools::map_dofs_to_support_points(mapping, dof_handler, support_points);
+  deallog << "Support points:" << std::endl;
+  for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
+    {
+      deallog << "  " << i << " (" << support_points[i] << ")";
+      if ((i + 1) % n_dof_1d == 0)
+        deallog << std::endl;
+    }
+
+  deallog << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  deallog.precision(2);
+
+  test<2>(1, 1);
+  test<2>(2, 1);
+
+  test<2>(1, 2);
+  test<2>(2, 2);
+
+  test<3>(2, 1);
+
+  return 0;
+}
diff --git a/tests/dofs/dof_renumbering_11.output b/tests/dofs/dof_renumbering_11.output
new file mode 100644 (file)
index 0000000..c3a2bc6
--- /dev/null
@@ -0,0 +1,65 @@
+
+DEAL::Testing dim=2, degree=1, components=1
+DEAL::Support points:
+DEAL::  0 (0.0 0.0)  1 (0.50 0.0)  2 (1.0 0.0)
+DEAL::  3 (0.0 0.50)  4 (0.50 0.50)  5 (1.0 0.50)
+DEAL::  6 (0.0 1.0)  7 (0.50 1.0)  8 (1.0 1.0)
+DEAL::
+DEAL::Testing dim=2, degree=2, components=1
+DEAL::Support points:
+DEAL::  0 (0.0 0.0)  1 (0.25 0.0)  2 (0.50 0.0)  3 (0.75 0.0)  4 (1.0 0.0)
+DEAL::  5 (0.0 0.25)  6 (0.25 0.25)  7 (0.50 0.25)  8 (0.75 0.25)  9 (1.0 0.25)
+DEAL::  10 (0.0 0.50)  11 (0.25 0.50)  12 (0.50 0.50)  13 (0.75 0.50)  14 (1.0 0.50)
+DEAL::  15 (0.0 0.75)  16 (0.25 0.75)  17 (0.50 0.75)  18 (0.75 0.75)  19 (1.0 0.75)
+DEAL::  20 (0.0 1.0)  21 (0.25 1.0)  22 (0.50 1.0)  23 (0.75 1.0)  24 (1.0 1.0)
+DEAL::
+DEAL::Testing dim=2, degree=1, components=2
+DEAL::Support points:
+DEAL::  0 (0.0 0.0)  1 (0.50 0.0)  2 (1.0 0.0)
+DEAL::  3 (0.0 0.50)  4 (0.50 0.50)  5 (1.0 0.50)
+DEAL::  6 (0.0 1.0)  7 (0.50 1.0)  8 (1.0 1.0)
+DEAL::  9 (0.0 0.0)  10 (0.50 0.0)  11 (1.0 0.0)
+DEAL::  12 (0.0 0.50)  13 (0.50 0.50)  14 (1.0 0.50)
+DEAL::  15 (0.0 1.0)  16 (0.50 1.0)  17 (1.0 1.0)
+DEAL::
+DEAL::Testing dim=2, degree=2, components=2
+DEAL::Support points:
+DEAL::  0 (0.0 0.0)  1 (0.25 0.0)  2 (0.50 0.0)  3 (0.75 0.0)  4 (1.0 0.0)
+DEAL::  5 (0.0 0.25)  6 (0.25 0.25)  7 (0.50 0.25)  8 (0.75 0.25)  9 (1.0 0.25)
+DEAL::  10 (0.0 0.50)  11 (0.25 0.50)  12 (0.50 0.50)  13 (0.75 0.50)  14 (1.0 0.50)
+DEAL::  15 (0.0 0.75)  16 (0.25 0.75)  17 (0.50 0.75)  18 (0.75 0.75)  19 (1.0 0.75)
+DEAL::  20 (0.0 1.0)  21 (0.25 1.0)  22 (0.50 1.0)  23 (0.75 1.0)  24 (1.0 1.0)
+DEAL::  25 (0.0 0.0)  26 (0.25 0.0)  27 (0.50 0.0)  28 (0.75 0.0)  29 (1.0 0.0)
+DEAL::  30 (0.0 0.25)  31 (0.25 0.25)  32 (0.50 0.25)  33 (0.75 0.25)  34 (1.0 0.25)
+DEAL::  35 (0.0 0.50)  36 (0.25 0.50)  37 (0.50 0.50)  38 (0.75 0.50)  39 (1.0 0.50)
+DEAL::  40 (0.0 0.75)  41 (0.25 0.75)  42 (0.50 0.75)  43 (0.75 0.75)  44 (1.0 0.75)
+DEAL::  45 (0.0 1.0)  46 (0.25 1.0)  47 (0.50 1.0)  48 (0.75 1.0)  49 (1.0 1.0)
+DEAL::
+DEAL::Testing dim=3, degree=2, components=1
+DEAL::Support points:
+DEAL::  0 (0.0 0.0 0.0)  1 (0.25 0.0 0.0)  2 (0.50 0.0 0.0)  3 (0.75 0.0 0.0)  4 (1.0 0.0 0.0)
+DEAL::  5 (0.0 0.25 0.0)  6 (0.25 0.25 0.0)  7 (0.50 0.25 0.0)  8 (0.75 0.25 0.0)  9 (1.0 0.25 0.0)
+DEAL::  10 (0.0 0.50 0.0)  11 (0.25 0.50 0.0)  12 (0.50 0.50 0.0)  13 (0.75 0.50 0.0)  14 (1.0 0.50 0.0)
+DEAL::  15 (0.0 0.75 0.0)  16 (0.25 0.75 0.0)  17 (0.50 0.75 0.0)  18 (0.75 0.75 0.0)  19 (1.0 0.75 0.0)
+DEAL::  20 (0.0 1.0 0.0)  21 (0.25 1.0 0.0)  22 (0.50 1.0 0.0)  23 (0.75 1.0 0.0)  24 (1.0 1.0 0.0)
+DEAL::  25 (0.0 0.0 0.25)  26 (0.25 0.0 0.25)  27 (0.50 0.0 0.25)  28 (0.75 0.0 0.25)  29 (1.0 0.0 0.25)
+DEAL::  30 (0.0 0.25 0.25)  31 (0.25 0.25 0.25)  32 (0.50 0.25 0.25)  33 (0.75 0.25 0.25)  34 (1.0 0.25 0.25)
+DEAL::  35 (0.0 0.50 0.25)  36 (0.25 0.50 0.25)  37 (0.50 0.50 0.25)  38 (0.75 0.50 0.25)  39 (1.0 0.50 0.25)
+DEAL::  40 (0.0 0.75 0.25)  41 (0.25 0.75 0.25)  42 (0.50 0.75 0.25)  43 (0.75 0.75 0.25)  44 (1.0 0.75 0.25)
+DEAL::  45 (0.0 1.0 0.25)  46 (0.25 1.0 0.25)  47 (0.50 1.0 0.25)  48 (0.75 1.0 0.25)  49 (1.0 1.0 0.25)
+DEAL::  50 (0.0 0.0 0.50)  51 (0.25 0.0 0.50)  52 (0.50 0.0 0.50)  53 (0.75 0.0 0.50)  54 (1.0 0.0 0.50)
+DEAL::  55 (0.0 0.25 0.50)  56 (0.25 0.25 0.50)  57 (0.50 0.25 0.50)  58 (0.75 0.25 0.50)  59 (1.0 0.25 0.50)
+DEAL::  60 (0.0 0.50 0.50)  61 (0.25 0.50 0.50)  62 (0.50 0.50 0.50)  63 (0.75 0.50 0.50)  64 (1.0 0.50 0.50)
+DEAL::  65 (0.0 0.75 0.50)  66 (0.25 0.75 0.50)  67 (0.50 0.75 0.50)  68 (0.75 0.75 0.50)  69 (1.0 0.75 0.50)
+DEAL::  70 (0.0 1.0 0.50)  71 (0.25 1.0 0.50)  72 (0.50 1.0 0.50)  73 (0.75 1.0 0.50)  74 (1.0 1.0 0.50)
+DEAL::  75 (0.0 0.0 0.75)  76 (0.25 0.0 0.75)  77 (0.50 0.0 0.75)  78 (0.75 0.0 0.75)  79 (1.0 0.0 0.75)
+DEAL::  80 (0.0 0.25 0.75)  81 (0.25 0.25 0.75)  82 (0.50 0.25 0.75)  83 (0.75 0.25 0.75)  84 (1.0 0.25 0.75)
+DEAL::  85 (0.0 0.50 0.75)  86 (0.25 0.50 0.75)  87 (0.50 0.50 0.75)  88 (0.75 0.50 0.75)  89 (1.0 0.50 0.75)
+DEAL::  90 (0.0 0.75 0.75)  91 (0.25 0.75 0.75)  92 (0.50 0.75 0.75)  93 (0.75 0.75 0.75)  94 (1.0 0.75 0.75)
+DEAL::  95 (0.0 1.0 0.75)  96 (0.25 1.0 0.75)  97 (0.50 1.0 0.75)  98 (0.75 1.0 0.75)  99 (1.0 1.0 0.75)
+DEAL::  100 (0.0 0.0 1.0)  101 (0.25 0.0 1.0)  102 (0.50 0.0 1.0)  103 (0.75 0.0 1.0)  104 (1.0 0.0 1.0)
+DEAL::  105 (0.0 0.25 1.0)  106 (0.25 0.25 1.0)  107 (0.50 0.25 1.0)  108 (0.75 0.25 1.0)  109 (1.0 0.25 1.0)
+DEAL::  110 (0.0 0.50 1.0)  111 (0.25 0.50 1.0)  112 (0.50 0.50 1.0)  113 (0.75 0.50 1.0)  114 (1.0 0.50 1.0)
+DEAL::  115 (0.0 0.75 1.0)  116 (0.25 0.75 1.0)  117 (0.50 0.75 1.0)  118 (0.75 0.75 1.0)  119 (1.0 0.75 1.0)
+DEAL::  120 (0.0 1.0 1.0)  121 (0.25 1.0 1.0)  122 (0.50 1.0 1.0)  123 (0.75 1.0 1.0)  124 (1.0 1.0 1.0)
+DEAL::
diff --git a/tests/fe/cell_face_lex.cc b/tests/fe/cell_face_lex.cc
new file mode 100644 (file)
index 0000000..aa39456
--- /dev/null
@@ -0,0 +1,91 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2009 - 2021 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.
+//
+// ------------------------------------------------------------------------
+
+
+
+// since early 2009, the FEValues objects try to be more efficient by only
+// recomputing things like gradients of shape functions if the cell on which
+// we are is not a translation of the previous one. in this series of tests we
+// make sure that this actually works the way it's supposed to be; in
+// particular, if we create a mesh of two identical cells but one has a curved
+// boundary, then they are the same if we use a Q1 mapping, but not a Q2
+// mapping. so we test that the mass matrix we get from each of these cells is
+// actually different in the latter case, but the same in the former
+//
+// Check cell to face lexicographic ordering
+
+
+#include <deal.II/fe/fe_tools.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test(unsigned int degree,
+     unsigned int direction,
+     bool         cell_hierarchical_numbering,
+     bool         is_continuous)
+{
+  std::pair<std::vector<unsigned int>, std::vector<unsigned int>> result =
+    FETools::cell_to_face_patch<dim>(degree,
+                                     direction,
+                                     cell_hierarchical_numbering,
+                                     is_continuous);
+
+  deallog << "Results for degree = " << degree << ", dim = " << dim
+          << ", direction = " << direction
+          << ", is_continuous = " << (is_continuous ? "true" : "false")
+          << ", cell_hierarchical_numbering = "
+          << (cell_hierarchical_numbering ? "true" : "false") << std::endl;
+
+  deallog << "Cell 0: ";
+  for (unsigned int val : result.first)
+    deallog << val << " ";
+  deallog << std::endl;
+
+  deallog << "Cell 1: ";
+  for (unsigned int val : result.second)
+    deallog << val << " ";
+  deallog << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(3);
+
+  test<2>(2, 0, false, false);
+  test<2>(2, 0, false, true);
+  test<2>(2, 0, true, true);
+
+  test<2>(4, 0, false, true);
+
+  test<2>(2, 1, false, false);
+  test<2>(2, 1, false, true);
+
+
+  test<3>(2, 0, false, false);
+  test<3>(2, 0, false, true);
+  test<3>(2, 0, true, true);
+
+  test<3>(2, 1, false, false);
+  test<3>(2, 1, false, true);
+
+  test<3>(2, 2, false, false);
+  test<3>(2, 2, false, true);
+}
diff --git a/tests/fe/cell_face_lex.output b/tests/fe/cell_face_lex.output
new file mode 100644 (file)
index 0000000..d8264bf
--- /dev/null
@@ -0,0 +1,40 @@
+
+DEAL::Results for degree = 2, dim = 2, direction = 0, is_continuous = false, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 6 7 8 12 13 14 
+DEAL::Cell 1: 3 4 5 9 10 11 15 16 17 
+DEAL::Results for degree = 2, dim = 2, direction = 0, is_continuous = true, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 5 6 7 10 11 12 
+DEAL::Cell 1: 2 3 4 7 8 9 12 13 14 
+DEAL::Results for degree = 2, dim = 2, direction = 0, is_continuous = true, cell_hierarchical_numbering = true
+DEAL::Cell 0: 0 2 10 12 5 7 1 11 6 
+DEAL::Cell 1: 2 4 12 14 7 9 3 13 8 
+DEAL::Results for degree = 4, dim = 2, direction = 0, is_continuous = true, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 3 4 9 10 11 12 13 18 19 20 21 22 27 28 29 30 31 36 37 38 39 40 
+DEAL::Cell 1: 4 5 6 7 8 13 14 15 16 17 22 23 24 25 26 31 32 33 34 35 40 41 42 43 44 
+DEAL::Results for degree = 2, dim = 2, direction = 1, is_continuous = false, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 3 4 5 6 7 8 
+DEAL::Cell 1: 9 10 11 12 13 14 15 16 17 
+DEAL::Results for degree = 2, dim = 2, direction = 1, is_continuous = true, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 3 4 5 6 7 8 
+DEAL::Cell 1: 6 7 8 9 10 11 12 13 14 
+DEAL::Results for degree = 2, dim = 3, direction = 0, is_continuous = false, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 6 7 8 12 13 14 18 19 20 24 25 26 30 31 32 36 37 38 42 43 44 48 49 50 
+DEAL::Cell 1: 3 4 5 9 10 11 15 16 17 21 22 23 27 28 29 33 34 35 39 40 41 45 46 47 51 52 53 
+DEAL::Results for degree = 2, dim = 3, direction = 0, is_continuous = true, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 5 6 7 10 11 12 15 16 17 20 21 22 25 26 27 30 31 32 35 36 37 40 41 42 
+DEAL::Cell 1: 2 3 4 7 8 9 12 13 14 17 18 19 22 23 24 27 28 29 32 33 34 37 38 39 42 43 44 
+DEAL::Results for degree = 2, dim = 3, direction = 0, is_continuous = true, cell_hierarchical_numbering = true
+DEAL::Cell 0: 0 2 10 12 30 32 40 42 5 7 1 11 35 37 31 41 15 17 25 27 20 22 16 26 6 36 21 
+DEAL::Cell 1: 2 4 12 14 32 34 42 44 7 9 3 13 37 39 33 43 17 19 27 29 22 24 18 28 8 38 23 
+DEAL::Results for degree = 2, dim = 3, direction = 1, is_continuous = false, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 3 4 5 6 7 8 18 19 20 21 22 23 24 25 26 36 37 38 39 40 41 42 43 44 
+DEAL::Cell 1: 9 10 11 12 13 14 15 16 17 27 28 29 30 31 32 33 34 35 45 46 47 48 49 50 51 52 53 
+DEAL::Results for degree = 2, dim = 3, direction = 1, is_continuous = true, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 3 4 5 6 7 8 15 16 17 18 19 20 21 22 23 30 31 32 33 34 35 36 37 38 
+DEAL::Cell 1: 6 7 8 9 10 11 12 13 14 21 22 23 24 25 26 27 28 29 36 37 38 39 40 41 42 43 44 
+DEAL::Results for degree = 2, dim = 3, direction = 2, is_continuous = false, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
+DEAL::Cell 1: 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 
+DEAL::Results for degree = 2, dim = 3, direction = 2, is_continuous = true, cell_hierarchical_numbering = false
+DEAL::Cell 0: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
+DEAL::Cell 1: 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

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.