]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add basic periodicity for tetrahedra with a mesh generator and a test. 18081/head
authorKyle Schwiebert <kjschwie@gmail.com>
Thu, 10 Apr 2025 00:15:11 +0000 (20:15 -0400)
committerKyle Schwiebert <kjschwie@gmail.com>
Thu, 10 Apr 2025 00:15:11 +0000 (20:15 -0400)
12 files changed:
doc/news/changes/minor/20250404KyleSchwiebert [new file with mode: 0644]
include/deal.II/fe/fe_simplex_p.h
include/deal.II/grid/grid_generator.h
source/fe/fe_simplex_p.cc
source/grid/grid_generator.cc
source/grid/grid_generator.inst.in
source/grid/grid_generator_from_name.cc
source/grid/grid_tools_dof_handlers.cc
tests/grid/grid_generator_from_name_and_argument_02.cc
tests/grid/grid_generator_from_name_and_argument_02.output
tests/simplex/periodic_tet.cc [new file with mode: 0644]
tests/simplex/periodic_tet.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20250404KyleSchwiebert b/doc/news/changes/minor/20250404KyleSchwiebert
new file mode 100644 (file)
index 0000000..f25ad59
--- /dev/null
@@ -0,0 +1,9 @@
+New: FE_Simplex_Poly now implements FiniteElement::face_to_cell_index(),
+enabling periodicity for derived classes, with restrictions similar to the
+implementation in FE_Q_Base. Additionally,
+GridGenerator::hyper_rectangle_with_simplicies() now accepts an additional
+boolean argument which, in 3D, breaks each cubical cell into six tetrahedra
+instead of five.  This mesh meets the restrictions of
+FE_Simplex_Poly::face_to_cell_index().
+<br>
+(Kyle Schwiebert, 2025/04/04)
index f493c16649210d15792c408f839a9262d94c6374..7ce1999318e1f4c180fdc8e22ab0debcc47262ef 100644 (file)
@@ -19,6 +19,7 @@
 
 #include <deal.II/base/mutex.h>
 #include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/base/types.h>
 
 #include <deal.II/fe/fe_poly.h>
 
@@ -65,6 +66,22 @@ public:
     const RefinementCase<dim> &refinement_case =
       RefinementCase<dim>::isotropic_refinement) const override;
 
+  /**
+   * @see FiniteElement::face_to_cell_index()
+   *
+   * @note This function has some limitations in 3D. If @p dim is 3, this
+   * function works if either:
+   * 1. @p combined_orientation is the default orientation.
+   * 2. @p combined_orientation is `(false, false, false)` and the triangular
+   * face has no more than one degree of freedom in its interior.
+   * Otherwise this function throws an error.
+   */
+  virtual unsigned int
+  face_to_cell_index(const unsigned int                 face_dof_index,
+                     const unsigned int                 face,
+                     const types::geometric_orientation combined_orientation =
+                       numbers::default_geometric_orientation) const override;
+
   /**
    * @copydoc dealii::FiniteElement::get_restriction_matrix()
    *
index 73a4415d6e01b5395a18a15f81722dd1a4012caa..c081235abc16e724992de926b1d0b67b231f3f12 100644 (file)
@@ -2651,9 +2651,14 @@ namespace GridGenerator
    * The number of vertices in coordinate
    * direction @p i is given by <tt>repetitions[i]+1</tt>.
    *
-   * @note This function takes the mesh produced by subdivided_hyper_rectangle()
-   * and further subdivides each cell into 2 triangles (for @p dim 2) or
-   * 5 tetrahedra (for @p dim 3), respectively.
+   * This function takes the mesh produced by subdivided_hyper_rectangle() and
+   * further subdivides each cell. For @p dim 2, it subdivides each cell into 2
+   * triangles.  For @p dim 3, it subdivides each cell into 5 or 6 tetrahedra
+   * based on the value of @p periodic. If @p periodic is true, then we split
+   * each cell into 6 cells so that each face of the rectangular prism has
+   * the same stencil, enabling periodicity. If @p periodic is false, we
+   * instead subdivide each hexahedral cell into 5 tetrahedra. If @p dim is not
+   * 3, then @p periodic has no effect.
    *
    * @note Currently, this function only works for `dim==spacedim`.
    *
@@ -2667,7 +2672,8 @@ namespace GridGenerator
     const std::vector<unsigned int> &repetitions,
     const Point<dim>                &p1,
     const Point<dim>                &p2,
-    const bool                       colorize = false);
+    const bool                       colorize = false,
+    const bool                       periodic = false);
 
   /**
    * Initialize the given triangulation with a hypercube (square in 2d and
@@ -2679,9 +2685,14 @@ namespace GridGenerator
    * the limits are given as arguments. They default to zero and unity, then
    * producing the unit hypercube.
    *
-   * @note This function takes the mesh produced by subdivided_hyper_cube()
-   * and further subdivides each cell into 2 triangles (for @p dim 2) or
-   * 5 tetrahedra (for @p dim 3), respectively.
+   * This function takes the mesh produced by subdivided_hyper_cube() and
+   * further subdivides each cell. For @p dim 2, it subdivides each cell into 2
+   * triangles.  For @p dim 3, it subdivides each cell into 5 or 6 tetrahedra
+   * based on the value of @p periodic. If @p periodic is true, then we split
+   * each cell into 6 cells so that each face of the rectangular prism has
+   * the same stencil, enabling periodicity. If @p periodic is false, we
+   * instead subdivide each hexahedral cell into 5 tetrahedra. If @p dim is not
+   * 3, then @p periodic has no effect.
    *
    * Also see
    * @ref simplex "Simplex support".
@@ -2692,7 +2703,8 @@ namespace GridGenerator
                                        const unsigned int repetitions,
                                        const double       p1       = 0.0,
                                        const double       p2       = 1.0,
-                                       const bool         colorize = false);
+                                       const bool         colorize = false,
+                                       const bool         periodic = false);
 
   /** @} */
 
index e958760b5c7d4b927a83a3e963883df5b70d2860..5c2a65e1fcb4ccc8daabdf5942aeb099785c191f 100644 (file)
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/exceptions.h>
 #include <deal.II/base/polynomials_barycentric.h>
 #include <deal.II/base/qprojector.h>
+#include <deal.II/base/types.h>
 
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_nothing.h>
@@ -24,6 +26,7 @@
 #include <deal.II/fe/fe_tools.h>
 #include <deal.II/fe/mapping.h>
 
+#include "deal.II/grid/tria_orientation.h"
 #include <deal.II/grid/grid_generator.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -436,6 +439,140 @@ FE_SimplexPoly<dim, spacedim>::get_prolongation_matrix(
 
 
 
+template <int dim, int spacedim>
+unsigned int
+FE_SimplexPoly<dim, spacedim>::face_to_cell_index(
+  const unsigned int                 face_dof_index,
+  const unsigned int                 face,
+  const types::geometric_orientation combined_orientation) const
+{
+  // Function largely lifted from FE_Q_Base::face_to_cell_index()
+  AssertIndexRange(face_dof_index, this->n_dofs_per_face(face));
+  AssertIndexRange(face, this->reference_cell().n_faces());
+
+  // TODO: once the default orientation is switched to 0 then we can remove this
+  // special case for 1D.
+  if (dim == 1)
+    Assert(combined_orientation == numbers::default_geometric_orientation,
+           ExcMessage("In 1D, all faces must have the default orientation."));
+  else
+    AssertIndexRange(combined_orientation,
+                     this->reference_cell().n_face_orientations(face));
+
+  // we need to distinguish between DoFs on vertices, lines and in 3d quads.
+  // do so in a sequence of if-else statements
+  if (face_dof_index < this->get_first_face_line_index(face))
+    // DoF is on a vertex
+    {
+      // get the number of the vertex on the face that corresponds to this DoF,
+      // along with the number of the DoF on this vertex
+      const unsigned int face_vertex =
+        face_dof_index / this->n_dofs_per_vertex();
+      const unsigned int dof_index_on_vertex =
+        face_dof_index % this->n_dofs_per_vertex();
+
+      // then get the number of this vertex on the cell and translate
+      // this to a DoF number on the cell
+      return this->reference_cell().face_to_cell_vertices(
+               face, face_vertex, combined_orientation) *
+               this->n_dofs_per_vertex() +
+             dof_index_on_vertex;
+    }
+  else if (face_dof_index < this->get_first_face_quad_index(face))
+    // DoF is on a line
+    {
+      // do the same kind of translation as before. we need to only consider
+      // DoFs on the lines, i.e., ignoring those on the vertices
+      const unsigned int index =
+        face_dof_index - this->get_first_face_line_index(face);
+
+      const unsigned int face_line         = index / this->n_dofs_per_line();
+      const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
+
+      // we now also need to adjust the line index for the case of
+      // face orientation, face flips, etc
+      unsigned int adjusted_dof_index_on_line = 0;
+      switch (dim)
+        {
+          case 1:
+            DEAL_II_ASSERT_UNREACHABLE();
+            break;
+
+          case 2:
+            if (combined_orientation == numbers::default_geometric_orientation)
+              adjusted_dof_index_on_line = dof_index_on_line;
+            else
+              adjusted_dof_index_on_line =
+                this->n_dofs_per_line() - dof_index_on_line - 1;
+            break;
+
+          case 3:
+            // in 3d, things are difficult. someone will have to think
+            // about how this code here should look like, by drawing a bunch
+            // of pictures of how all the faces can look like with the various
+            // flips and rotations.
+            //
+            // that said, we can implement a couple other situations easily:
+            // if the face orientation is
+            // numbers::default_combined_face_orientation then things are
+            // simple. Likewise if the face orientation is
+            // internal::combined_face_orientation(false,false,false) then we
+            // just know we need to reverse the DoF order on each line.
+            //
+            // this at least allows for tetrahedral meshes where periodic face
+            // pairs are each in one of the two above configurations. Doing so
+            // may turn out to be always possible or it may not, but at least
+            // this is less restrictive.
+            Assert((this->n_dofs_per_line() <= 1) ||
+                     combined_orientation ==
+                       numbers::default_geometric_orientation ||
+                     combined_orientation ==
+                       internal::combined_face_orientation(false, false, false),
+                   ExcNotImplemented());
+
+            if (combined_orientation == numbers::default_geometric_orientation)
+              {
+                adjusted_dof_index_on_line = dof_index_on_line;
+              }
+            else
+              {
+                adjusted_dof_index_on_line =
+                  this->n_dofs_per_line() - dof_index_on_line - 1;
+              }
+            break;
+
+          default:
+            DEAL_II_ASSERT_UNREACHABLE();
+        }
+
+      return (this->get_first_line_index() +
+              this->reference_cell().face_to_cell_lines(face,
+                                                        face_line,
+                                                        combined_orientation) *
+                this->n_dofs_per_line() +
+              adjusted_dof_index_on_line);
+    }
+  else
+    // DoF is on a quad
+    {
+      Assert(dim >= 3, ExcInternalError());
+
+      // ignore vertex and line dofs
+      const unsigned int index =
+        face_dof_index - this->get_first_face_quad_index(face);
+
+      // the same is true here as above for the 3d case -- someone will
+      // just have to draw a bunch of pictures. in the meantime,
+      // we can implement the degree <= 3 case in which it is simple
+      Assert((this->n_dofs_per_quad(face) <= 1) ||
+               combined_orientation == numbers::default_geometric_orientation,
+             ExcNotImplemented());
+      return (this->get_first_quad_index(face) + index);
+    }
+}
+
+
+
 template <int dim, int spacedim>
 const FullMatrix<double> &
 FE_SimplexPoly<dim, spacedim>::get_restriction_matrix(
index 5acf678de89558ec9e4d43b9455ee4196c8572e1..40d87b361496bcaceb0702f05c43dedee56d5597 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/grid/filtered_iterator.h>
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_topology.h>
 #include <deal.II/grid/intergrid_map.h>
 #include <deal.II/grid/manifold_lib.h>
 #include <deal.II/grid/tria.h>
@@ -8940,179 +8941,348 @@ namespace GridGenerator
 
 
 
-  template <int dim, int spacedim>
-  void
-  subdivided_hyper_rectangle_with_simplices(
-    Triangulation<dim, spacedim>    &tria,
-    const std::vector<unsigned int> &repetitions,
-    const Point<dim>                &p1,
-    const Point<dim>                &p2,
-    const bool                       colorize)
+  // Hide the implementation for two cases of
+  // subdivided_hyper_rectangle_with_simplices in an anonymous namespace.
+  namespace
   {
-    AssertDimension(dim, spacedim);
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_rectangle_with_simplices_no_periodic(
+      Triangulation<dim, spacedim>    &tria,
+      const std::vector<unsigned int> &repetitions,
+      const Point<dim>                &p1,
+      const Point<dim>                &p2,
+      const bool                       colorize)
+    {
+      AssertDimension(dim, spacedim);
 
-    std::vector<Point<spacedim>> vertices;
-    std::vector<CellData<dim>>   cells;
+      std::vector<Point<spacedim>> vertices;
+      std::vector<CellData<dim>>   cells;
 
-    if (dim == 2)
-      {
-        // determine cell sizes
-        const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
-                            (p2[1] - p1[1]) / repetitions[1]);
+      if (dim == 2)
+        {
+          // determine cell sizes
+          const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+                              (p2[1] - p1[1]) / repetitions[1]);
+
+          // create vertices
+          for (unsigned int j = 0; j <= repetitions[1]; ++j)
+            for (unsigned int i = 0; i <= repetitions[0]; ++i)
+              vertices.push_back(
+                Point<spacedim>(p1[0] + dx[0] * i, p1[1] + dx[1] * j));
+
+          // create cells
+          for (unsigned int j = 0; j < repetitions[1]; ++j)
+            for (unsigned int i = 0; i < repetitions[0]; ++i)
+              {
+                // create reference QUAD cell
+                std::array<unsigned int, 4> quad{{
+                  (j + 0) * (repetitions[0] + 1) + i + 0, //
+                  (j + 0) * (repetitions[0] + 1) + i + 1, //
+                  (j + 1) * (repetitions[0] + 1) + i + 0, //
+                  (j + 1) * (repetitions[0] + 1) + i + 1  //
+                }};                                       //
+
+                // TRI cell 0
+                {
+                  CellData<dim> tri;
+                  tri.vertices = {quad[0], quad[1], quad[2]};
+                  cells.push_back(tri);
+                }
+
+                // TRI cell 1
+                {
+                  CellData<dim> tri;
+                  tri.vertices = {quad[3], quad[2], quad[1]};
+                  cells.push_back(tri);
+                }
+              }
+        }
+      else if (dim == 3)
+        {
+          // determine cell sizes
+          const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+                              (p2[1] - p1[1]) / repetitions[1],
+                              (p2[2] - p1[2]) / repetitions[2]);
+
+          // create vertices
+          for (unsigned int k = 0; k <= repetitions[2]; ++k)
+            for (unsigned int j = 0; j <= repetitions[1]; ++j)
+              for (unsigned int i = 0; i <= repetitions[0]; ++i)
+                vertices.push_back(Point<spacedim>(p1[0] + dx[0] * i,
+                                                   p1[1] + dx[1] * j,
+                                                   p1[2] + dx[2] * k));
+
+          // create cells
+          for (unsigned int k = 0; k < repetitions[2]; ++k)
+            for (unsigned int j = 0; j < repetitions[1]; ++j)
+              for (unsigned int i = 0; i < repetitions[0]; ++i)
+                {
+                  // create reference HEX cell
+                  std::array<unsigned int, 8> quad{
+                    {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 0) * (repetitions[0] + 1) + i + 0,
+                     (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 0) * (repetitions[0] + 1) + i + 1,
+                     (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 1) * (repetitions[0] + 1) + i + 0,
+                     (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 1) * (repetitions[0] + 1) + i + 1,
+                     (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 0) * (repetitions[0] + 1) + i + 0,
+                     (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 0) * (repetitions[0] + 1) + i + 1,
+                     (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 1) * (repetitions[0] + 1) + i + 0,
+                     (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                       (j + 1) * (repetitions[0] + 1) + i + 1}};
+
+                  // TET cell 0
+                  {
+                    CellData<dim> cell;
+                    if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+                      cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}};
+                    else
+                      cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}};
+
+                    cells.push_back(cell);
+                  }
 
-        // create vertices
+                  // TET cell 1
+                  {
+                    CellData<dim> cell;
+                    if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+                      cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}};
+                    else
+                      cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}};
+                    cells.push_back(cell);
+                  }
+
+                  // TET cell 2
+                  {
+                    CellData<dim> cell;
+                    if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+                      cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}};
+                    else
+                      cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}};
+                    cells.push_back(cell);
+                  }
+
+                  // TET cell 3
+                  {
+                    CellData<dim> cell;
+                    if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+                      cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}};
+                    else
+                      cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}};
+                    cells.push_back(cell);
+                  }
+
+                  // TET cell 4
+                  {
+                    CellData<dim> cell;
+                    if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+                      cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}};
+                    else
+                      cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}};
+                    cells.push_back(cell);
+                  }
+                }
+        }
+      else
+        {
+          AssertThrow(false, ExcNotImplemented());
+        }
+
+      // actually create triangulation
+      tria.create_triangulation(vertices, cells, SubCellData());
+
+      if (colorize)
+        {
+          // to colorize, run through all
+          // faces of all cells and set
+          // boundary indicator to the
+          // correct value if it was 0.
+
+          // use a large epsilon to
+          // compare numbers to avoid
+          // roundoff problems.
+          double epsilon = std::numeric_limits<double>::max();
+          for (unsigned int i = 0; i < dim; ++i)
+            epsilon =
+              std::min(epsilon,
+                       0.01 * (std::abs(p2[i] - p1[i]) / repetitions[i]));
+          Assert(epsilon > 0,
+                 ExcMessage(
+                   "The distance between corner points must be positive."));
+
+          // actual code is external since
+          // 1-D is different from 2/3d.
+          colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
+        }
+    }
+
+
+
+    // This function is only needed in 3D.
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_rectangle_with_simplices_periodic(
+      Triangulation<dim, spacedim>    &tria,
+      const std::vector<unsigned int> &repetitions,
+      const Point<dim>                &p1,
+      const Point<dim>                &p2,
+      const bool                       colorize)
+    {
+      // This function is only needed in 3D (and hypothetically in higher
+      // dimension), so library internals should ensure it is never called
+      // unless dim == 3.
+      Assert(dim == 3, ExcInternalError());
+      AssertDimension(dim, spacedim);
+
+      std::vector<Point<spacedim>> vertices;
+      std::vector<CellData<dim>>   cells;
+
+      // determine cell sizes
+      const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+                          (p2[1] - p1[1]) / repetitions[1],
+                          (p2[2] - p1[2]) / repetitions[2]);
+
+      // create vertices
+      for (unsigned int k = 0; k <= repetitions[2]; ++k)
         for (unsigned int j = 0; j <= repetitions[1]; ++j)
           for (unsigned int i = 0; i <= repetitions[0]; ++i)
-            vertices.push_back(
-              Point<spacedim>(p1[0] + dx[0] * i, p1[1] + dx[1] * j));
+            vertices.push_back(Point<spacedim>(p1[0] + dx[0] * i,
+                                               p1[1] + dx[1] * j,
+                                               p1[2] + dx[2] * k));
 
-        // create cells
+      // create cells
+      for (unsigned int k = 0; k < repetitions[2]; ++k)
         for (unsigned int j = 0; j < repetitions[1]; ++j)
           for (unsigned int i = 0; i < repetitions[0]; ++i)
             {
-              // create reference QUAD cell
-              std::array<unsigned int, 4> quad{{
-                (j + 0) * (repetitions[0] + 1) + i + 0, //
-                (j + 0) * (repetitions[0] + 1) + i + 1, //
-                (j + 1) * (repetitions[0] + 1) + i + 0, //
-                (j + 1) * (repetitions[0] + 1) + i + 1  //
-              }};                                       //
-
-              // TRI cell 0
+              // create reference HEX cell
+              std::array<unsigned int, 8> quad{
+                {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 0) * (repetitions[0] + 1) + i + 0,
+                 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 0) * (repetitions[0] + 1) + i + 1,
+                 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 1) * (repetitions[0] + 1) + i + 0,
+                 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 1) * (repetitions[0] + 1) + i + 1,
+                 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 0) * (repetitions[0] + 1) + i + 0,
+                 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 0) * (repetitions[0] + 1) + i + 1,
+                 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 1) * (repetitions[0] + 1) + i + 0,
+                 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+                   (j + 1) * (repetitions[0] + 1) + i + 1}};
+
+              // TET cell 0
               {
-                CellData<dim> tri;
-                tri.vertices = {quad[0], quad[1], quad[2]};
-                cells.push_back(tri);
+                CellData<dim> cell;
+                cell.vertices = {{quad[0], quad[1], quad[3], quad[7]}};
+                cells.push_back(cell);
               }
 
-              // TRI cell 1
+              // TET cell 1
               {
-                CellData<dim> tri;
-                tri.vertices = {quad[3], quad[2], quad[1]};
-                cells.push_back(tri);
+                CellData<dim> cell;
+                cell.vertices = {{quad[0], quad[1], quad[7], quad[5]}};
+                cells.push_back(cell);
               }
-            }
-      }
-    else if (dim == 3)
-      {
-        // determine cell sizes
-        const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
-                            (p2[1] - p1[1]) / repetitions[1],
-                            (p2[2] - p1[2]) / repetitions[2]);
 
-        // create vertices
-        for (unsigned int k = 0; k <= repetitions[2]; ++k)
-          for (unsigned int j = 0; j <= repetitions[1]; ++j)
-            for (unsigned int i = 0; i <= repetitions[0]; ++i)
-              vertices.push_back(Point<spacedim>(p1[0] + dx[0] * i,
-                                                 p1[1] + dx[1] * j,
-                                                 p1[2] + dx[2] * k));
+              // TET cell 2
+              {
+                CellData<dim> cell;
+                cell.vertices = {{quad[0], quad[7], quad[3], quad[2]}};
+                cells.push_back(cell);
+              }
 
-        // create cells
-        for (unsigned int k = 0; k < repetitions[2]; ++k)
-          for (unsigned int j = 0; j < repetitions[1]; ++j)
-            for (unsigned int i = 0; i < repetitions[0]; ++i)
+              // TET cell 3
               {
-                // create reference HEX cell
-                std::array<unsigned int, 8> quad{
-                  {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 0) * (repetitions[0] + 1) + i + 0,
-                   (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 0) * (repetitions[0] + 1) + i + 1,
-                   (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 1) * (repetitions[0] + 1) + i + 0,
-                   (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 1) * (repetitions[0] + 1) + i + 1,
-                   (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 0) * (repetitions[0] + 1) + i + 0,
-                   (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 0) * (repetitions[0] + 1) + i + 1,
-                   (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 1) * (repetitions[0] + 1) + i + 0,
-                   (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
-                     (j + 1) * (repetitions[0] + 1) + i + 1}};
-
-                // TET cell 0
-                {
-                  CellData<dim> cell;
-                  if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
-                    cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}};
-                  else
-                    cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}};
+                CellData<dim> cell;
+                cell.vertices = {{quad[2], quad[6], quad[0], quad[7]}};
+                cells.push_back(cell);
+              }
 
-                  cells.push_back(cell);
-                }
+              // TET cell 4
+              {
+                CellData<dim> cell;
+                cell.vertices = {{quad[4], quad[7], quad[5], quad[0]}};
+                cells.push_back(cell);
+              }
 
-                // TET cell 1
-                {
-                  CellData<dim> cell;
-                  if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
-                    cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}};
-                  else
-                    cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}};
-                  cells.push_back(cell);
-                }
+              // TET cell 5
+              {
+                CellData<dim> cell;
+                cell.vertices = {{quad[4], quad[6], quad[7], quad[0]}};
+                cells.push_back(cell);
+              }
+            }
 
-                // TET cell 2
-                {
-                  CellData<dim> cell;
-                  if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
-                    cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}};
-                  else
-                    cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}};
-                  cells.push_back(cell);
-                }
+      // actually create triangulation
+      tria.create_triangulation(vertices, cells, SubCellData());
 
-                // TET cell 3
-                {
-                  CellData<dim> cell;
-                  if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
-                    cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}};
-                  else
-                    cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}};
-                  cells.push_back(cell);
-                }
+      if (colorize)
+        {
+          // to colorize, run through all
+          // faces of all cells and set
+          // boundary indicator to the
+          // correct value if it was 0.
+
+          // use a large epsilon to
+          // compare numbers to avoid
+          // roundoff problems.
+          double epsilon = std::numeric_limits<double>::max();
+          for (unsigned int i = 0; i < dim; ++i)
+            epsilon =
+              std::min(epsilon,
+                       0.01 * (std::abs(p2[i] - p1[i]) / repetitions[i]));
+          Assert(epsilon > 0,
+                 ExcMessage(
+                   "The distance between corner points must be positive."));
 
-                // TET cell 4
-                {
-                  CellData<dim> cell;
-                  if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
-                    cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}};
-                  else
-                    cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}};
-                  cells.push_back(cell);
-                }
-              }
+          // actual code is external since
+          // 1-D is different from 2/3d.
+          colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
+        }
+    }
+  } // namespace
+
+
+
+  template <int dim, int spacedim>
+  void
+  subdivided_hyper_rectangle_with_simplices(
+    Triangulation<dim, spacedim>    &tria,
+    const std::vector<unsigned int> &repetitions,
+    const Point<dim>                &p1,
+    const Point<dim>                &p2,
+    const bool                       colorize,
+    const bool                       periodic)
+  {
+    // We only need to call the "periodic" variant if it was requested and we
+    // are in 3D.
+    if (dim != 3)
+      {
+        subdivided_hyper_rectangle_with_simplices_no_periodic(
+          tria, repetitions, p1, p2, colorize);
+        return;
       }
-    else
+    else if (!periodic)
       {
-        AssertThrow(false, ExcNotImplemented());
+        subdivided_hyper_rectangle_with_simplices_no_periodic(
+          tria, repetitions, p1, p2, colorize);
+        return;
       }
-
-    // actually create triangulation
-    tria.create_triangulation(vertices, cells, SubCellData());
-
-    if (colorize)
+    else
       {
-        // to colorize, run through all
-        // faces of all cells and set
-        // boundary indicator to the
-        // correct value if it was 0.
-
-        // use a large epsilon to
-        // compare numbers to avoid
-        // roundoff problems.
-        double epsilon = std::numeric_limits<double>::max();
-        for (unsigned int i = 0; i < dim; ++i)
-          epsilon = std::min(epsilon,
-                             0.01 * (std::abs(p2[i] - p1[i]) / repetitions[i]));
-        Assert(epsilon > 0,
-               ExcMessage(
-                 "The distance between corner points must be positive."));
-
-        // actual code is external since
-        // 1-D is different from 2/3d.
-        colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
+        subdivided_hyper_rectangle_with_simplices_periodic(
+          tria, repetitions, p1, p2, colorize);
+        return;
       }
   }
 
@@ -9124,12 +9294,17 @@ namespace GridGenerator
                                        const unsigned int repetitions,
                                        const double       p1,
                                        const double       p2,
-                                       const bool         colorize)
+                                       const bool         colorize,
+                                       const bool         periodic)
   {
     if (dim == 2)
       {
-        subdivided_hyper_rectangle_with_simplices(
-          tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize);
+        subdivided_hyper_rectangle_with_simplices(tria,
+                                                  {{repetitions, repetitions}},
+                                                  {p1, p1},
+                                                  {p2, p2},
+                                                  colorize,
+                                                  periodic);
       }
     else if (dim == 3)
       {
@@ -9138,7 +9313,8 @@ namespace GridGenerator
           {{repetitions, repetitions, repetitions}},
           {p1, p1, p1},
           {p2, p2, p2},
-          colorize);
+          colorize,
+          periodic);
       }
     else
       {
index 189ce57a84bb7ef32ba9adb8114c3cb789f730e3..90d15de511b35062eee86db71604c537355daa5a 100644 (file)
@@ -302,12 +302,14 @@ for (deal_II_dimension : DIMENSIONS)
       const std::vector<unsigned int> &repetitions,
       const Point<deal_II_dimension>  &p1,
       const Point<deal_II_dimension>  &p2,
-      const bool                       colorize);
+      const bool                       colorize,
+      const bool                       periodic);
 
     template void GridGenerator::subdivided_hyper_cube_with_simplices(
       Triangulation<deal_II_dimension> & tria,
       const unsigned int repetitions,
       const double       p1,
       const double       p2,
-      const bool         colorize);
+      const bool         colorize,
+      const bool         periodic);
   }
index 361915d5f26086c6e130946c77389ea96483d5a1..bc0b89a38fc887635168890f59b1b6c7f63572c2 100644 (file)
@@ -232,7 +232,7 @@ namespace GridGenerator
                          bool>(concentric_hyper_shells, arguments, tria);
 
       else if (name == "subdivided_hyper_cube_with_simplices")
-        parse_and_create<dim, dim, unsigned int, double, double, bool>(
+        parse_and_create<dim, dim, unsigned int, double, double, bool, bool>(
           subdivided_hyper_cube_with_simplices, arguments, tria);
 
       else if (name == "subdivided_hyper_rectangle_with_simplices")
@@ -241,6 +241,7 @@ namespace GridGenerator
                          const std::vector<unsigned int> &,
                          const Point<dim> &,
                          const Point<dim> &,
+                         bool,
                          bool>(subdivided_hyper_rectangle_with_simplices,
                                arguments,
                                tria);
index 009e8d32743bb59d4df3e8b7c40962e473962b97..3f7da7e914eb96831044688a1c4cea56318f7bd4 100644 (file)
@@ -15,6 +15,7 @@
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
+#include <deal.II/base/types.h>
 
 #include <deal.II/distributed/fully_distributed_tria.h>
 #include <deal.II/distributed/shared_tria.h>
@@ -2432,19 +2433,17 @@ namespace GridTools
   {
     Assert(matrix.m() == matrix.n(),
            ExcMessage("The supplied matrix must be a square matrix"));
-
-    static const int dim = FaceIterator::AccessorType::dimension;
+    Assert(face1->reference_cell() == face2->reference_cell(),
+           ExcMessage(
+             "The faces to be matched must have the same reference cell."));
 
     // Do a full matching of the face vertices:
-
-    std::array<unsigned int, GeometryInfo<dim>::vertices_per_face>
-      face1_vertices, face2_vertices;
-
-    face1_vertices.fill(numbers::invalid_unsigned_int);
-    face2_vertices.fill(numbers::invalid_unsigned_int);
-
     AssertDimension(face1->n_vertices(), face2->n_vertices());
 
+    std::vector<unsigned int> face1_vertices(face1->n_vertices(),
+                                             numbers::invalid_unsigned_int),
+      face2_vertices(face2->n_vertices(), numbers::invalid_unsigned_int);
+
     std::set<unsigned int> face2_vertices_set;
     for (unsigned int i = 0; i < face1->n_vertices(); ++i)
       face2_vertices_set.insert(i);
index 2883b09d6e0fb6b135a4b838a319e3bcd292f93b..40d9b1929743f1b9da7c68779fe0da46a7139cc3 100644 (file)
@@ -54,13 +54,15 @@ main()
   test<2, 2>("hyper_ball_balanced", "0,0 : 1");
   test<3, 3>("hyper_ball_balanced", "0,0,0 : 1");
 
-  test<2, 2>("subdivided_hyper_cube_with_simplices", "2 : 0.0 : 1.0 : false");
-  test<3, 3>("subdivided_hyper_cube_with_simplices", "2 : 0.0 : 1.0 : false");
+  test<2, 2>("subdivided_hyper_cube_with_simplices",
+             "2 : 0.0 : 1.0 : false : false");
+  test<3, 3>("subdivided_hyper_cube_with_simplices",
+             "2 : 0.0 : 1.0 : false : false");
 
   test<2, 2>("subdivided_hyper_rectangle_with_simplices",
-             "2, 2 : 0.0, 0.0 : 1.0, 2.0 : false");
+             "2, 2 : 0.0, 0.0 : 1.0, 2.0 : false : false");
   test<3, 3>("subdivided_hyper_rectangle_with_simplices",
-             "2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false");
+             "2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false : false");
 
   test<2, 2>("subdivided_hyper_L", "5, 5 : 0, 0 : 1, 1 : 2, 3");
   test<3, 3>("subdivided_hyper_L", "5, 5, 5 : 0, 0, 0 : 1, 1, 1 : 2, 2, 3");
index e5a39f8d22d620f85b0f1d5ee3662186acfd5eb0..31b233c3f52ec1e61fe1390dc5d04b41f5e0148d 100644 (file)
@@ -191,7 +191,7 @@ $ELM
 31 5 0 0 8 18 16 38 37 27 25 52 53 
 32 5 0 0 8 43 36 48 50 37 38 52 53 
 $ENDELM
-DEAL::Generating Triangulation<2, 2> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false)
+DEAL::Generating Triangulation<2, 2> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false : false)
 $NOD
 9
 1  0.00000 0.00000 0
@@ -215,7 +215,7 @@ $ELM
 7 2 0 0 3 5 6 8 
 8 2 0 0 3 9 8 6 
 $ENDELM
-DEAL::Generating Triangulation<3, 3> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false)
+DEAL::Generating Triangulation<3, 3> : subdivided_hyper_cube_with_simplices(2 : 0.0 : 1.0 : false : false)
 $NOD
 27
 1  0.00000 0.00000 0.00000
@@ -289,7 +289,7 @@ $ELM
 39 4 0 0 4 18 24 27 26 
 40 4 0 0 4 14 18 26 24 
 $ENDELM
-DEAL::Generating Triangulation<2, 2> : subdivided_hyper_rectangle_with_simplices(2, 2 : 0.0, 0.0 : 1.0, 2.0 : false)
+DEAL::Generating Triangulation<2, 2> : subdivided_hyper_rectangle_with_simplices(2, 2 : 0.0, 0.0 : 1.0, 2.0 : false : false)
 $NOD
 9
 1  0.00000 0.00000 0
@@ -313,7 +313,7 @@ $ELM
 7 2 0 0 3 5 6 8 
 8 2 0 0 3 9 8 6 
 $ENDELM
-DEAL::Generating Triangulation<3, 3> : subdivided_hyper_rectangle_with_simplices(2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false)
+DEAL::Generating Triangulation<3, 3> : subdivided_hyper_rectangle_with_simplices(2, 2, 3 : 0.0, 0.0, 1.0 : 1.0, 2.0, 3.0 : false : false)
 $NOD
 36
 1  0.00000 0.00000 1.00000
diff --git a/tests/simplex/periodic_tet.cc b/tests/simplex/periodic_tet.cc
new file mode 100644 (file)
index 0000000..d3745af
--- /dev/null
@@ -0,0 +1,300 @@
+/* ------------------------------------------------------------------------
+ *
+ * SPDX-License-Identifier: LGPL-2.1-or-later
+ * Copyright (C) 2020 - 2024 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 that the function
+// GridGenerator::subdivided_hyper_rectable_with_simplices() can give a mesh
+// which is suitible for periodicity with tetrahedra.
+
+#include <deal.II/base/types.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_orientation.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include "../tests.h"
+
+void
+check()
+{
+  Triangulation<3>          tria;
+  Point<3>                  corner_p1     = {-0.5, 1.0, 0.0};
+  Point<3>                  corner_p2     = {0.5, 2.0, 1.2};
+  std::vector<unsigned int> sub_divisions = {2, 2, 2};
+  GridGenerator::subdivided_hyper_rectangle_with_simplices(
+    tria, sub_divisions, corner_p1, corner_p2, true, true);
+
+  // Check that we have made a grid with the expected number of cells and only
+  // one level of refinement.
+  if (tria.n_cells() != 48)
+    {
+      deallog << "FAIL! Wrong number of cells." << std::endl;
+    }
+
+  if (tria.n_levels() != 1)
+    {
+      deallog << "FAIL! Mesh has more than one refinement level." << std::endl;
+    }
+
+  // Check that each cell has a positive measure.
+  {
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        double measure = cell->measure();
+        if (measure <= 0.0)
+          {
+            deallog << "FAIL! Found a cell with nonpositive measure: "
+                    << cell->active_cell_index() << std::endl;
+          }
+      }
+  }
+
+  // Collect periodic face pairs and ensure that they make sense and we have
+  // the right numbers of them.
+  std::array<std::vector<GridTools::PeriodicFacePair<
+               typename Triangulation<3>::cell_iterator>>,
+             3>
+    pairs;
+  GridTools::collect_periodic_faces(tria, 0, 1, 0, pairs[0]);
+  GridTools::collect_periodic_faces(tria, 2, 3, 1, pairs[1]);
+  GridTools::collect_periodic_faces(tria, 4, 5, 2, pairs[2]);
+  tria.add_periodicity(pairs[0]);
+  tria.add_periodicity(pairs[1]);
+  tria.add_periodicity(pairs[2]);
+
+  // To enable periodicity the easy way, we need each boundary pair to be in
+  // default orientation. Check that this holds.
+  {
+    IteratorFilters::AtBoundary        filter;
+    const types::geometric_orientation default_orientation =
+      numbers::default_geometric_orientation;
+    const types::geometric_orientation default_mirror_orientation =
+      internal::combined_face_orientation(false, false, false);
+    for (const auto &pair : pairs)
+      {
+        for (const auto &match : pair)
+          {
+            if (match.orientation != default_orientation &&
+                match.orientation != default_mirror_orientation)
+              {
+                deallog << "FAIL! Found a face pair which is not in default "
+                        << "orientation." << std::endl;
+                deallog << "Orientation is "
+                        << static_cast<unsigned int>(match.orientation)
+                        << " vs default "
+                        << static_cast<unsigned int>(default_orientation)
+                        << std::endl;
+                deallog << "Cells " << match.cell[0]->active_cell_index()
+                        << " and " << match.cell[1]->active_cell_index()
+                        << std::endl;
+                deallog << "Faces " << match.face_idx[0] << " and "
+                        << match.face_idx[1] << std::endl;
+                return;
+              }
+          }
+      }
+  }
+
+  {
+    // Check that we have gotten all the periodicity.
+    bool         mismatched_faces   = false;
+    unsigned int changing_direction = 0;
+    for (const auto &pair : pairs)
+      {
+        if (pair.size() != 8)
+          {
+            deallog << "FAIL! There should be 8 matches on each face."
+                    << std::endl;
+            return;
+          }
+        for (int face = 0; face < static_cast<int>(pair.size()); ++face)
+          {
+            Point<3> p1(
+              pair[face].cell[0]->face(pair[face].face_idx[0])->center());
+            Point<3> p2(
+              pair[face].cell[1]->face(pair[face].face_idx[1])->center());
+            const std::array<bool, 3> results = {
+              {std::abs(p1[0] - p2[0]) > 1e-10,
+               std::abs(p1[1] - p2[1]) > 1e-10,
+               std::abs(p1[2] - p2[2]) > 1e-10}};
+            bool fail = false;
+
+            for (unsigned int i = 0; i < 3; ++i)
+              {
+                if (i == changing_direction)
+                  {
+                    continue;
+                  }
+                else
+                  {
+                    fail = fail || results[i];
+                  }
+              }
+
+            if (fail)
+              {
+                deallog << "FAIL! Found non matching faces." << std::endl;
+                return;
+              }
+          }
+        ++changing_direction;
+      }
+  }
+
+  // write 2 outputs (total mesh and only surface mesh)
+  const auto grid_out = [](const auto &tria, const bool surface_mesh_only) {
+    GridOutFlags::Vtk flags;
+
+    if (surface_mesh_only)
+      {
+        flags.output_cells         = false;
+        flags.output_faces         = true;
+        flags.output_edges         = false;
+        flags.output_only_relevant = false;
+      }
+
+    GridOut grid_out;
+    grid_out.set_flags(flags);
+
+    grid_out.write_vtk(tria, deallog.get_file_stream());
+  };
+
+  grid_out(tria, false); // total mesh
+  grid_out(tria, true);  // only surface mesh
+
+  // If we got this far, we have a mesh which supports periodicity.
+  // Let's see if we can proscribe periodicity with actual DoFs now.
+  // Check for C0 tetrahedra up to degree 3. (Highest implemented at the time of
+  // writing.)
+  for (const unsigned int degree : {1, 2, 3})
+    {
+      DoFHandler<3>  dof_handler(tria);
+      FE_SimplexP<3> fe(degree);
+      dof_handler.distribute_dofs(fe);
+      std::array<AffineConstraints<double>, 3> constraints;
+
+      // Periodic in the x direction:
+      DoFTools::make_periodicity_constraints(
+        dof_handler, 0, 1, 0, constraints[0]);
+      constraints[0].close();
+
+      // Periodic in the y direction:
+      DoFTools::make_periodicity_constraints(
+        dof_handler, 2, 3, 1, constraints[1]);
+      constraints[1].close();
+
+      // Periodic in the z direction:
+      DoFTools::make_periodicity_constraints(
+        dof_handler, 4, 5, 2, constraints[2]);
+      constraints[2].close();
+
+      // Some information for checking constraints.
+      MappingFE<3> mapping(FE_SimplexP<3>(1));
+      auto map = DoFTools::map_dofs_to_support_points(mapping, dof_handler);
+      unsigned int changing_direction = 0;
+
+      for (const auto &cons : constraints)
+        {
+          // Check the total number of constraints.
+          if (cons.n_constraints() != (2 * degree + 1) * (2 * degree + 1))
+            {
+              deallog << "FAIL! Wrong number of constraints!" << std::endl;
+            }
+
+          // Loop through the constraints and verify that the correct components
+          // of the DoF's support points are equal.
+          for (const auto &constraint : cons.get_lines())
+            {
+              if (constraint.entries.size() == 1)
+                {
+                  if (std::abs(constraint.entries[0].second - 1.0) > 1e-10)
+                    {
+                      deallog << "FAIL! All constraint entries should be one."
+                              << std::endl;
+                      continue;
+                    }
+
+                  const types::global_dof_index dof_1 = constraint.index;
+                  const types::global_dof_index dof_2 =
+                    constraint.entries[0].first;
+                  const std::array<bool, 3> results = {
+                    {std::abs(map[dof_1][0] - map[dof_2][0]) > 1e-10,
+                     std::abs(map[dof_1][1] - map[dof_2][1]) > 1e-10,
+                     std::abs(map[dof_1][2] - map[dof_2][2]) > 1e-10}};
+                  bool fail = false;
+
+                  for (unsigned int i = 0; i < 3; ++i)
+                    {
+                      if (i == changing_direction)
+                        {
+                          continue;
+                        }
+                      else
+                        {
+                          fail = fail || results[i];
+                        }
+                    }
+
+                  if (fail)
+                    {
+                      deallog << "FAIL! Found non matching constraint."
+                              << std::endl;
+                      deallog << "(" << map[dof_1][0] << "," << map[dof_1][1]
+                              << "," << map[dof_1][2] << ")" << std::endl;
+                      deallog << "(" << map[dof_2][0] << "," << map[dof_2][1]
+                              << "," << map[dof_2][2] << ")" << std::endl;
+                    }
+                }
+              else
+                {
+                  deallog
+                    << "FAIL! "
+                    << "Found constraint with the wrong number of entries."
+                    << std::endl;
+                  return;
+                }
+            }
+          ++changing_direction;
+
+          // If degree is 3, put down a little more data for the output.
+          if (degree == 3)
+            {
+              cons.print(deallog.get_file_stream());
+            }
+        }
+    }
+
+  deallog << "OK!" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  check();
+}
diff --git a/tests/simplex/periodic_tet.output b/tests/simplex/periodic_tet.output
new file mode 100644 (file)
index 0000000..ca3a39b
--- /dev/null
@@ -0,0 +1,385 @@
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 27 double
+-0.500000 1.00000 0.00000
+0.00000 1.00000 0.00000
+0.500000 1.00000 0.00000
+-0.500000 1.50000 0.00000
+0.00000 1.50000 0.00000
+0.500000 1.50000 0.00000
+-0.500000 2.00000 0.00000
+0.00000 2.00000 0.00000
+0.500000 2.00000 0.00000
+-0.500000 1.00000 0.600000
+0.00000 1.00000 0.600000
+0.500000 1.00000 0.600000
+-0.500000 1.50000 0.600000
+0.00000 1.50000 0.600000
+0.500000 1.50000 0.600000
+-0.500000 2.00000 0.600000
+0.00000 2.00000 0.600000
+0.500000 2.00000 0.600000
+-0.500000 1.00000 1.20000
+0.00000 1.00000 1.20000
+0.500000 1.00000 1.20000
+-0.500000 1.50000 1.20000
+0.00000 1.50000 1.20000
+0.500000 1.50000 1.20000
+-0.500000 2.00000 1.20000
+0.00000 2.00000 1.20000
+0.500000 2.00000 1.20000
+
+CELLS 88 400
+4 0 1 4 13
+4 0 1 13 10
+4 0 13 4 3
+4 3 12 0 13
+4 9 13 10 0
+4 9 12 13 0
+4 1 2 5 14
+4 1 2 14 11
+4 1 14 5 4
+4 4 13 1 14
+4 10 14 11 1
+4 10 13 14 1
+4 3 4 7 16
+4 3 4 16 13
+4 3 16 7 6
+4 6 15 3 16
+4 12 16 13 3
+4 12 15 16 3
+4 4 5 8 17
+4 4 5 17 14
+4 4 17 8 7
+4 7 16 4 17
+4 13 17 14 4
+4 13 16 17 4
+4 9 10 13 22
+4 9 10 22 19
+4 9 22 13 12
+4 12 21 9 22
+4 18 22 19 9
+4 18 21 22 9
+4 10 11 14 23
+4 10 11 23 20
+4 10 23 14 13
+4 13 22 10 23
+4 19 23 20 10
+4 19 22 23 10
+4 12 13 16 25
+4 12 13 25 22
+4 12 25 16 15
+4 15 24 12 25
+4 21 25 22 12
+4 21 24 25 12
+4 13 14 17 26
+4 13 14 26 23
+4 13 26 17 16
+4 16 25 13 26
+4 22 26 23 13
+4 22 25 26 13
+3 1 0 10
+3 0 1 4
+3 2 1 11
+3 1 2 5
+3 5 2 14
+3 14 2 11
+3 0 4 3
+3 3 4 7
+3 1 5 4
+3 4 5 8
+3 8 5 17
+3 17 5 14
+3 15 6 16
+3 3 7 6
+3 7 16 6
+3 16 7 17
+3 4 8 7
+3 8 17 7
+3 9 10 0
+3 10 9 19
+3 10 11 1
+3 11 10 20
+3 23 11 20
+3 14 11 23
+3 26 14 23
+3 24 15 25
+3 16 25 15
+3 25 16 26
+3 17 14 26
+3 17 26 16
+3 18 21 22
+3 18 22 19
+3 18 19 9
+3 19 23 20
+3 19 20 10
+3 21 24 25
+3 21 25 22
+3 19 22 23
+3 22 26 23
+3 22 25 26
+
+CELL_TYPES 88
+10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 
+5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
+
+
+CELL_DATA 88
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+6 6 6 6 6 6 7 7 7 7 7 7 6 6 6 6 6 6 7 7 7 7 7 7 6 6 6 6 6 6 7 7 7 7 7 7 6 6 6 6 6 6 7 7 7 7 7 7 
+2 4 2 4 1 1 4 4 4 4 1 1 3 4 3 3 4 3 2 2 2 2 1 1 1 3 3 3 1 3 5 5 2 5 2 5 5 5 5 5 
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 27 double
+-0.500000 1.00000 0.00000
+0.00000 1.00000 0.00000
+0.500000 1.00000 0.00000
+-0.500000 1.50000 0.00000
+0.00000 1.50000 0.00000
+0.500000 1.50000 0.00000
+-0.500000 2.00000 0.00000
+0.00000 2.00000 0.00000
+0.500000 2.00000 0.00000
+-0.500000 1.00000 0.600000
+0.00000 1.00000 0.600000
+0.500000 1.00000 0.600000
+-0.500000 1.50000 0.600000
+0.00000 1.50000 0.600000
+0.500000 1.50000 0.600000
+-0.500000 2.00000 0.600000
+0.00000 2.00000 0.600000
+0.500000 2.00000 0.600000
+-0.500000 1.00000 1.20000
+0.00000 1.00000 1.20000
+0.500000 1.00000 1.20000
+-0.500000 1.50000 1.20000
+0.00000 1.50000 1.20000
+0.500000 1.50000 1.20000
+-0.500000 2.00000 1.20000
+0.00000 2.00000 1.20000
+0.500000 2.00000 1.20000
+
+CELLS 48 192
+3 1 0 10
+3 0 1 4
+3 2 1 11
+3 1 2 5
+3 5 2 14
+3 14 2 11
+3 3 12 0
+3 0 4 3
+3 3 4 7
+3 1 5 4
+3 4 5 8
+3 8 5 17
+3 17 5 14
+3 6 15 3
+3 15 6 16
+3 3 7 6
+3 7 16 6
+3 16 7 17
+3 4 8 7
+3 8 17 7
+3 12 9 0
+3 9 10 0
+3 10 9 19
+3 10 11 1
+3 11 10 20
+3 23 11 20
+3 15 12 3
+3 12 21 9
+3 14 11 23
+3 26 14 23
+3 15 24 12
+3 24 15 25
+3 16 25 15
+3 25 16 26
+3 17 14 26
+3 17 26 16
+3 21 18 9
+3 18 21 22
+3 18 22 19
+3 18 19 9
+3 19 23 20
+3 19 20 10
+3 24 21 12
+3 21 24 25
+3 21 25 22
+3 19 22 23
+3 22 26 23
+3 22 25 26
+
+CELL_TYPES 48
+5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
+
+
+CELL_DATA 48
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+2 4 2 4 1 1 0 4 4 4 4 1 1 0 3 4 3 3 4 3 0 2 2 2 2 1 0 0 1 1 0 3 3 3 1 3 0 5 5 2 5 2 0 5 5 5 5 5 
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
+    64 0:  1.00000
+    65 30:  1.00000
+    66 40:  1.00000
+    69 31:  1.00000
+    70 32:  1.00000
+    75 44:  1.00000
+    76 43:  1.00000
+    77 41:  1.00000
+    78 42:  1.00000
+    82 47:  1.00000
+    83 50:  1.00000
+    86 56:  1.00000
+    87 55:  1.00000
+    88 61:  1.00000
+    89 60:  1.00000
+    92 63:  1.00000
+    160 132:  1.00000
+    161 142:  1.00000
+    162 133:  1.00000
+    163 134:  1.00000
+    168 146:  1.00000
+    169 145:  1.00000
+    170 143:  1.00000
+    171 144:  1.00000
+    175 149:  1.00000
+    176 157:  1.00000
+    177 156:  1.00000
+    179 159:  1.00000
+    244 220:  1.00000
+    247 224:  1.00000
+    248 223:  1.00000
+    249 221:  1.00000
+    250 222:  1.00000
+    253 227:  1.00000
+    254 230:  1.00000
+    257 236:  1.00000
+    258 235:  1.00000
+    259 241:  1.00000
+    260 240:  1.00000
+    263 243:  1.00000
+    316 298:  1.00000
+    319 302:  1.00000
+    320 301:  1.00000
+    321 299:  1.00000
+    322 300:  1.00000
+    325 305:  1.00000
+    326 313:  1.00000
+    327 312:  1.00000
+    329 315:  1.00000
+    112 1:  1.00000
+    113 20:  1.00000
+    122 23:  1.00000
+    123 24:  1.00000
+    132 0:  1.00000
+    135 22:  1.00000
+    136 21:  1.00000
+    137 5:  1.00000
+    138 4:  1.00000
+    141 27:  1.00000
+    142 50:  1.00000
+    143 56:  1.00000
+    144 55:  1.00000
+    147 54:  1.00000
+    148 53:  1.00000
+    150 59:  1.00000
+    160 64:  1.00000
+    161 83:  1.00000
+    170 86:  1.00000
+    171 87:  1.00000
+    180 85:  1.00000
+    181 84:  1.00000
+    182 68:  1.00000
+    183 67:  1.00000
+    186 90:  1.00000
+    187 107:  1.00000
+    188 106:  1.00000
+    189 110:  1.00000
+    280 206:  1.00000
+    285 209:  1.00000
+    286 210:  1.00000
+    294 208:  1.00000
+    295 207:  1.00000
+    297 213:  1.00000
+    298 230:  1.00000
+    299 236:  1.00000
+    300 235:  1.00000
+    303 234:  1.00000
+    304 233:  1.00000
+    306 239:  1.00000
+    316 254:  1.00000
+    321 257:  1.00000
+    322 258:  1.00000
+    330 256:  1.00000
+    331 255:  1.00000
+    333 261:  1.00000
+    334 275:  1.00000
+    335 274:  1.00000
+    336 278:  1.00000
+    196 2:  1.00000
+    206 1:  1.00000
+    211 7:  1.00000
+    212 6:  1.00000
+    220 30:  1.00000
+    225 36:  1.00000
+    226 35:  1.00000
+    230 0:  1.00000
+    231 9:  1.00000
+    232 8:  1.00000
+    233 5:  1.00000
+    234 4:  1.00000
+    237 16:  1.00000
+    240 31:  1.00000
+    241 32:  1.00000
+    242 38:  1.00000
+    244 65:  1.00000
+    254 64:  1.00000
+    259 70:  1.00000
+    260 69:  1.00000
+    268 96:  1.00000
+    269 95:  1.00000
+    272 72:  1.00000
+    273 71:  1.00000
+    274 68:  1.00000
+    275 67:  1.00000
+    276 79:  1.00000
+    279 98:  1.00000
+    280 112:  1.00000
+    290 115:  1.00000
+    291 114:  1.00000
+    298 132:  1.00000
+    303 138:  1.00000
+    304 137:  1.00000
+    308 117:  1.00000
+    309 116:  1.00000
+    310 124:  1.00000
+    312 133:  1.00000
+    313 134:  1.00000
+    314 140:  1.00000
+    316 160:  1.00000
+    326 163:  1.00000
+    327 162:  1.00000
+    334 183:  1.00000
+    335 182:  1.00000
+    338 165:  1.00000
+    339 164:  1.00000
+    340 172:  1.00000
+    342 185:  1.00000
+DEAL::OK!

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.