]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make GridGenerator::hyper_ball() work in dim=2, spacedim=3.
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 23 Jan 2025 22:48:37 +0000 (15:48 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 27 Jan 2025 17:56:17 +0000 (10:56 -0700)
include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc
source/grid/grid_generator.inst.in

index 90f9c9fa98a92bcab0a4d6055949f7becc53f2d6..e574a82fbd1cec318c799d4b5e65a09cd6f4d955 100644 (file)
@@ -730,12 +730,12 @@ namespace GridGenerator
   /**
    * Initialize the given triangulation with several
    * @ref GlossCoarseMesh "coarse mesh cells"
-   * that cover a hyperball, i.e. a circle in 2d or a
-   * ball in 3d, around @p center with given @p radius. The function is
+   * that cover a hyperball, i.e. a circular disk if `dim==2` or a
+   * ball if `dim==3`, around @p center with given @p radius. The function is
    * used in step-6.
    *
-   * In order to avoid degenerate cells at the boundaries, the circle is
-   * triangulated by five cells, whereas in 3d the ball is subdivided by
+   * In order to avoid degenerate cells at the boundaries, the circular disk is
+   * triangulated by five cells, whereas in 3d the ball is subdivided into
    * seven cells. Specifically, these
    * cells are one cell in the center plus one "cap" cell on each of the faces
    * of this center cell. This ensures that under repeated refinement, none
@@ -747,7 +747,10 @@ namespace GridGenerator
    * after one refinement is optimized.
    *
    * This function is declared to exist for triangulations of all space
-   * dimensions, but throws an error if called in 1d.
+   * dimensions, but throws an error if called if `dim==1`. If `spacedim>dim`
+   * (which is only possible if `dim==2, spacedim==3`) then the function
+   * creates a circular disk for which all vertices have a $z$ value equal
+   * to the $z$ coordinate of the center point provided by the second argument.
    *
    * By default, the manifold_id is set to 0 on the boundary faces, 1 on the
    * boundary cells, and numbers::flat_manifold_id on the central cell and on
@@ -806,11 +809,11 @@ namespace GridGenerator
    * @note This function is available through the python interface as
    * `triangulation.generate_hyper_ball(point, radius = 1.)`.
    */
-  template <int dim>
+  template <int dim, int spacedim>
   void
-  hyper_ball(Triangulation<dim> &tria,
-             const Point<dim>   &center = Point<dim>(),
-             const double        radius = 1.,
+  hyper_ball(Triangulation<dim, spacedim> &tria,
+             const Point<spacedim>        &center                   = {},
+             const double                  radius                   = 1.,
              const bool attach_spherical_manifold_on_boundary_cells = false);
 
   /**
index f15ecb478d5a79847a08b403ad1c10644ce2bbdb..e1b25bc712e0442c214198fbd65b01587ab6592d 100644 (file)
@@ -3929,15 +3929,6 @@ namespace GridGenerator
 
 
 
-  template <>
-  void
-  hyper_ball(Triangulation<1> &, const Point<1> &, const double, const bool)
-  {
-    DEAL_II_NOT_IMPLEMENTED();
-  }
-
-
-
   template <>
   void
   hyper_ball_balanced(Triangulation<1> &, const Point<1> &, const double)
@@ -4293,48 +4284,125 @@ namespace GridGenerator
 
 
 
-  // Implementation for 2d only
-  template <>
+  template <int dim, int spacedim>
   void
-  hyper_ball(Triangulation<2> &tria,
-             const Point<2>   &p,
-             const double      radius,
-             const bool        internal_manifolds)
+  hyper_ball(Triangulation<dim, spacedim> &tria,
+             const Point<spacedim>        &p,
+             const double                  radius,
+             const bool                    internal_manifolds)
   {
-    // equilibrate cell sizes at
-    // transition from the inner part
-    // to the radial cells
-    const double   a           = 1. / (1 + std::sqrt(2.0));
-    const Point<2> vertices[8] = {
-      p + Point<2>(-1, -1) * (radius / std::sqrt(2.0)),
-      p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)),
-      p + Point<2>(-1, -1) * (radius / std::sqrt(2.0) * a),
-      p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a),
-      p + Point<2>(-1, +1) * (radius / std::sqrt(2.0) * a),
-      p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a),
-      p + Point<2>(-1, +1) * (radius / std::sqrt(2.0)),
-      p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))};
+    if constexpr (dim == 2)
+      {
+        const auto embed_point = [](const double x,
+                                    const double y) -> Point<spacedim> {
+          if constexpr (spacedim == 2)
+            return {x, y};
+          else if constexpr (spacedim == 3)
+            return {x, y, 0};
+          else
+            DEAL_II_NOT_IMPLEMENTED();
+        };
 
-    std::vector<CellData<2>> cells(5, CellData<2>());
 
-    for (unsigned int i = 0; i < 5; ++i)
-      {
-        for (unsigned int j = 0; j < 4; ++j)
-          cells[i].vertices[j] = circle_cell_vertices[i][j];
-        cells[i].material_id = 0;
-        cells[i].manifold_id = i == 2 ? numbers::flat_manifold_id : 1;
-      }
+        // Equilibrate cell sizes at transition from the inner part
+        // to the radial cells
+        const double          a           = 1. / (1 + std::sqrt(2.0));
+        const Point<spacedim> vertices[8] = {
+          p + embed_point(-1, -1) * (radius / std::sqrt(2.0)),
+          p + embed_point(+1, -1) * (radius / std::sqrt(2.0)),
+          p + embed_point(-1, -1) * (radius / std::sqrt(2.0) * a),
+          p + embed_point(+1, -1) * (radius / std::sqrt(2.0) * a),
+          p + embed_point(-1, +1) * (radius / std::sqrt(2.0) * a),
+          p + embed_point(+1, +1) * (radius / std::sqrt(2.0) * a),
+          p + embed_point(-1, +1) * (radius / std::sqrt(2.0)),
+          p + embed_point(+1, +1) * (radius / std::sqrt(2.0))};
 
-    tria.create_triangulation(std::vector<Point<2>>(std::begin(vertices),
-                                                    std::end(vertices)),
-                              cells,
-                              SubCellData()); // no boundary information
-    tria.set_all_manifold_ids_on_boundary(0);
-    tria.set_manifold(0, SphericalManifold<2>(p));
-    if (internal_manifolds)
-      tria.set_manifold(1, SphericalManifold<2>(p));
+        std::vector<CellData<2>> cells(5, CellData<2>());
+
+        for (unsigned int i = 0; i < 5; ++i)
+          {
+            for (unsigned int j = 0; j < 4; ++j)
+              cells[i].vertices[j] = circle_cell_vertices[i][j];
+            cells[i].material_id = 0;
+            cells[i].manifold_id = i == 2 ? numbers::flat_manifold_id : 1;
+          }
+
+        tria.create_triangulation(std::vector<Point<spacedim>>(
+                                    std::begin(vertices), std::end(vertices)),
+                                  cells,
+                                  SubCellData()); // no boundary information
+        tria.set_all_manifold_ids_on_boundary(0);
+        tria.set_manifold(0, SphericalManifold<dim, spacedim>(p));
+        if (internal_manifolds)
+          tria.set_manifold(1, SphericalManifold<dim, spacedim>(p));
+        else
+          tria.set_manifold(1, FlatManifold<dim, spacedim>());
+      }
+    else if constexpr (dim == 3)
+      {
+        const double a =
+          1. / (1 + std::sqrt(3.0)); // equilibrate cell sizes at transition
+        // from the inner part to the radial
+        // cells
+        const unsigned int n_vertices           = 16;
+        const Point<3>     vertices[n_vertices] = {
+          // first the vertices of the inner
+          // cell
+          p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0) * a),
+          p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0) * a),
+          p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0) * a),
+          p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0) * a),
+          p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0) * a),
+          p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0) * a),
+          p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0) * a),
+          p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0) * a),
+          // now the eight vertices at
+          // the outer sphere
+          p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0)),
+          p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0)),
+          p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0)),
+          p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0)),
+          p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0)),
+          p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0)),
+          p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0)),
+          p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0)),
+        };
+
+        // one needs to draw the seven cubes to
+        // understand what's going on here
+        const unsigned int n_cells                   = 7;
+        const int          cell_vertices[n_cells][8] = {
+          {0, 1, 4, 5, 3, 2, 7, 6},      // center
+          {8, 9, 12, 13, 0, 1, 4, 5},    // bottom
+          {9, 13, 1, 5, 10, 14, 2, 6},   // right
+          {11, 10, 3, 2, 15, 14, 7, 6},  // top
+          {8, 0, 12, 4, 11, 3, 15, 7},   // left
+          {8, 9, 0, 1, 11, 10, 3, 2},    // front
+          {12, 4, 13, 5, 15, 7, 14, 6}}; // back
+
+        std::vector<CellData<3>> cells(n_cells, CellData<3>());
+
+        for (unsigned int i = 0; i < n_cells; ++i)
+          {
+            for (const unsigned int j : GeometryInfo<3>::vertex_indices())
+              cells[i].vertices[j] = cell_vertices[i][j];
+            cells[i].material_id = 0;
+            cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1;
+          }
+
+        tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
+                                                        std::end(vertices)),
+                                  cells,
+                                  SubCellData()); // no boundary information
+        tria.set_all_manifold_ids_on_boundary(0);
+        tria.set_manifold(0, SphericalManifold<3>(p));
+        if (internal_manifolds)
+          tria.set_manifold(1, SphericalManifold<3>(p));
+        else
+          tria.set_manifold(1, FlatManifold<3>());
+      }
     else
-      tria.set_manifold(1, FlatManifold<2>());
+      DEAL_II_NOT_IMPLEMENTED();
   }
 
 
@@ -5079,78 +5147,6 @@ namespace GridGenerator
 
 
 
-  // Implementation for 3d only
-  template <>
-  void
-  hyper_ball(Triangulation<3> &tria,
-             const Point<3>   &p,
-             const double      radius,
-             const bool        internal_manifold)
-  {
-    const double a =
-      1. / (1 + std::sqrt(3.0)); // equilibrate cell sizes at transition
-    // from the inner part to the radial
-    // cells
-    const unsigned int n_vertices           = 16;
-    const Point<3>     vertices[n_vertices] = {
-      // first the vertices of the inner
-      // cell
-      p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0) * a),
-      p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0) * a),
-      p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0) * a),
-      p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0) * a),
-      p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0) * a),
-      p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0) * a),
-      p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0) * a),
-      p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0) * a),
-      // now the eight vertices at
-      // the outer sphere
-      p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0)),
-      p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0)),
-      p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0)),
-      p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0)),
-      p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0)),
-      p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0)),
-      p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0)),
-      p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0)),
-    };
-
-    // one needs to draw the seven cubes to
-    // understand what's going on here
-    const unsigned int n_cells                   = 7;
-    const int          cell_vertices[n_cells][8] = {
-      {0, 1, 4, 5, 3, 2, 7, 6},      // center
-      {8, 9, 12, 13, 0, 1, 4, 5},    // bottom
-      {9, 13, 1, 5, 10, 14, 2, 6},   // right
-      {11, 10, 3, 2, 15, 14, 7, 6},  // top
-      {8, 0, 12, 4, 11, 3, 15, 7},   // left
-      {8, 9, 0, 1, 11, 10, 3, 2},    // front
-      {12, 4, 13, 5, 15, 7, 14, 6}}; // back
-
-    std::vector<CellData<3>> cells(n_cells, CellData<3>());
-
-    for (unsigned int i = 0; i < n_cells; ++i)
-      {
-        for (const unsigned int j : GeometryInfo<3>::vertex_indices())
-          cells[i].vertices[j] = cell_vertices[i][j];
-        cells[i].material_id = 0;
-        cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1;
-      }
-
-    tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
-                                                    std::end(vertices)),
-                              cells,
-                              SubCellData()); // no boundary information
-    tria.set_all_manifold_ids_on_boundary(0);
-    tria.set_manifold(0, SphericalManifold<3>(p));
-    if (internal_manifold)
-      tria.set_manifold(1, SphericalManifold<3>(p));
-    else
-      tria.set_manifold(1, FlatManifold<3>());
-  }
-
-
-
   void
   non_standard_orientation_mesh(Triangulation<2>  &tria,
                                 const unsigned int n_rotate_middle_square)
index 75a7d06b57795da0ba738bdf7ca3e6cf16efd0f2..189ce57a84bb7ef32ba9adb8114c3cb789f730e3 100644 (file)
@@ -40,7 +40,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         const double,
         const bool);
 
-
       template void
       subdivided_hyper_rectangle<deal_II_dimension, deal_II_space_dimension>(
         Triangulation<deal_II_dimension, deal_II_space_dimension> &,
@@ -72,6 +71,13 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         const Point<deal_II_dimension> &,
         const std::vector<int> &);
 
+      template void
+      hyper_ball<deal_II_dimension, deal_II_space_dimension>(
+        Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+        const Point<deal_II_space_dimension> &,
+        const double,
+        const bool);
+
       template void
       cheese<deal_II_dimension, deal_II_space_dimension>(
         Triangulation<deal_II_dimension, deal_II_space_dimension> &,

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.