Add Simplex::GridGenerator

diff --git a/include/deal.II/simplex/grid_generator.h b/include/deal.II/simplex/grid_generator.h
new file mode 100644 (file)
index 0000000..f35f98e
--- /dev/null
+++ b/include/deal.II/simplex/grid_generator.h
@@ -0,0 +1,248 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_simplex_grid_generator_h
+#define dealii_simplex_grid_generator_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+  /**
+   * This namespace provides a collection of functions to generate simplex
+   * triangulations for some basic geometries.
+   */
+  namespace GridGenerator
+  {
+    /**
+     * Create a coordinate-parallel brick from the two diagonally opposite
+     * corner points @p p1 and @p p2. The number of vertices in coordinate
+     * direction @p i is given by <tt>repetitions[i]+1</tt>.
+     *
+     * @note This function connects internally 4/8 vertices to quadrilateral/
+     *   hexahedral cells and subdivides these into 2/5 triangular/
+     *   tetrahedral cells.
+     *
+     * @note Currently, this function only works for `dim==spacedim`.
+     */
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_rectangle(Triangulation<dim, spacedim> &   tria,
+                               const std::vector<unsigned int> &repetitions,
+                               const Point<dim> &               p1,
+                               const Point<dim> &               p2,
+                               const bool colorize = false)
+    {
+      AssertDimension(dim, spacedim);
+
+      AssertThrow(colorize == false, ExcNotImplemented());
+
+      std::vector<Point<spacedim>> vertices;
+      std::vector<CellData<dim>>   cells;
+
+      if (dim == 2)
+        {
+          // determine cell sizes
+          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
+          Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+                        (p2[1] - p1[1]) / repetitions[1],
+                        (p2[2] - p1[2]) / repetitions[1]);
+
+          // 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);
+                  }
+
+                  // 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(colorize == false, ExcNotImplemented());
+        }
+
+      // actually create triangulation
+      tria.create_triangulation(vertices, cells, SubCellData());
+    }
+
+    /**
+     * Initialize the given triangulation with a hypercube (square in 2D and
+     * cube in 3D) consisting of @p repetitions cells in each direction.
+     * The hypercube volume is the tensor product interval
+     * $[left,right]^{\text{dim}}$ in the present number of dimensions, where
+     * the limits are given as arguments. They default to zero and unity, then
+     * producing the unit hypercube.
+     *
+     * @note This function connects internally 4/8 vertices to quadrilateral/
+     *   hexahedral cells and subdivides these into 2/5 triangular/
+     *   tetrahedral cells.
+     */
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_cube(Triangulation<dim, spacedim> &tria,
+                          const unsigned int            repetitions,
+                          const double                  p1       = 0.0,
+                          const double                  p2       = 1.0,
+                          const bool                    colorize = false)
+    {
+      if (dim == 2)
+        {
+          subdivided_hyper_rectangle(
+            tria, {repetitions, repetitions}, {p1, p1}, {p2, p2}, colorize);
+        }
+      else if (dim == 3)
+        {
+          subdivided_hyper_rectangle(tria,
+                                     {repetitions, repetitions, repetitions},
+                                     {p1, p1, p1},
+                                     {p2, p2, p2},
+                                     colorize);
+        }
+      else
+        {
+          AssertThrow(false, ExcNotImplemented())
+        }
+    }
+  } // namespace GridGenerator
+} // namespace Simplex
+
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif