Add alternative GridGenerator::hyper_ball with better balance of cells

diff --git a/doc/doxygen/images/hyper_ball_balanced_2d.png b/doc/doxygen/images/hyper_ball_balanced_2d.png
new file mode 100644 (file)
index 0000000..4d41bd6
Binary files /dev/null and b/doc/doxygen/images/hyper_ball_balanced_2d.png differ

diff --git a/doc/doxygen/images/hyper_ball_balanced_3d.png b/doc/doxygen/images/hyper_ball_balanced_3d.png
new file mode 100644 (file)
index 0000000..8a5d9a4
Binary files /dev/null and b/doc/doxygen/images/hyper_ball_balanced_3d.png differ

diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h
index 0684b4403395d7161abbd9c71b3174b75fc8d6dc..5ce7c04e2b55afa2c68ab8f0b81998129fe128eb 100644 (file)
--- a/include/deal.II/grid/grid_generator.h
+++ b/include/deal.II/grid/grid_generator.h
@@ -738,6 +738,38 @@ namespace GridGenerator
              const double        radius = 1.,
              const bool attach_spherical_manifold_on_boundary_cells = false);
 
+  /**
+   * This is an alternative to hyper_ball with 12 cells in 2d and 32 cells in
+   * 3d, which provides a better balance between the size of the cells around
+   * the outer curved boundaries and the cell in the interior. The mesh is
+   * based on the cells used by GridGenerator::quarter_hyper_ball() with
+   * appropriate copies and rotations to fill the whole ball.
+   *
+   * The following pictures show the resulting mesh in 2D (left) and 3D:
+   * <table align="center" class="doxtable">
+   *   <tr>
+   *     <td>
+   *       <img src="hyper_ball_balanced_2d.png" alt="" width="40%">
+   *     </td>
+   *     <td>
+   *       <img src="hyper_ball_balanced_3d.png" alt="" width="40%">
+   *     </td>
+   *   </tr>
+   * </table>
+   *
+   * By default, the manifold_id is set to 0 on the boundary faces, 1 on the
+   * boundary cells, and types::flat_manifold_id on the central cell and on
+   * internal faces.
+   *
+   * @pre The triangulation passed as argument needs to be empty when calling
+   * this function.
+   */
+  template <int dim>
+  void
+  hyper_ball_balanced(Triangulation<dim> &tria,
+                      const Point<dim> &  center = Point<dim>(),
+                      const double        radius = 1.);
+
   /**
    * Creates a hyper sphere, i.e., a surface of a ball in @p spacedim
    * dimensions. This function only exists for dim+1=spacedim in 2 and 3 space

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index c31e8151157670b65cd38cde00d14098ad2ea90c..c96b628a1d4ba18780b956ef84860a63d65431a7 100644 (file)
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -5333,6 +5333,99 @@ namespace GridGenerator
   }
 
 
+
+  template <int dim>
+  void
+  hyper_ball_balanced(Triangulation<dim> &tria,
+                      const Point<dim> &  p,
+                      const double        radius)
+  {
+    // We create the ball by duplicating the information in each dimension at
+    // a time by appropriate rotations, starting from the quarter ball. The
+    // rotations make sure we do not generate inverted cells that would appear
+    // if we tried the slightly simpler approach to simply mirror the cells.
+
+    Triangulation<dim> tria_piece;
+    GridGenerator::quarter_hyper_ball(tria_piece, p, radius);
+
+    for (unsigned int round = 0; round < dim; ++round)
+      {
+        Triangulation<dim> tria_copy;
+        tria_copy.copy_triangulation(tria_piece);
+        tria_piece.clear();
+        std::vector<Point<dim>> new_points(tria_copy.n_vertices());
+        if (round == 0)
+          for (unsigned int v = 0; v < tria_copy.n_vertices(); ++v)
+            {
+              // rotate by 90 degrees counterclockwise
+              new_points[v][0] = -tria_copy.get_vertices()[v][1];
+              new_points[v][1] = tria_copy.get_vertices()[v][0];
+              if (dim == 3)
+                new_points[v][2] = tria_copy.get_vertices()[v][2];
+            }
+        else if (round == 1)
+          {
+            for (unsigned int v = 0; v < tria_copy.n_vertices(); ++v)
+              {
+                // rotate by 180 degrees along the xy plane
+                new_points[v][0] = -tria_copy.get_vertices()[v][0];
+                new_points[v][1] = -tria_copy.get_vertices()[v][1];
+                if (dim == 3)
+                  new_points[v][2] = tria_copy.get_vertices()[v][2];
+              }
+          }
+        else if (round == 2)
+          for (unsigned int v = 0; v < tria_copy.n_vertices(); ++v)
+            {
+              // rotate by 180 degrees along the xz plane
+              Assert(dim == 3, ExcInternalError());
+              new_points[v][0] = -tria_copy.get_vertices()[v][0];
+              new_points[v][1] = tria_copy.get_vertices()[v][1];
+              new_points[v][2] = -tria_copy.get_vertices()[v][2];
+            }
+        else
+          Assert(false, ExcInternalError());
+
+
+        // the cell data is exactly the same as before
+        std::vector<CellData<dim>> cells;
+        cells.reserve(tria_copy.n_cells());
+        for (const auto &cell : tria_copy.cell_iterators())
+          {
+            CellData<dim> data;
+            for (unsigned int v : GeometryInfo<dim>::vertex_indices())
+              data.vertices[v] = cell->vertex_index(v);
+            data.material_id = cell->material_id();
+            data.manifold_id = cell->manifold_id();
+            cells.push_back(data);
+          }
+
+        Triangulation<dim> rotated_tria;
+        rotated_tria.create_triangulation(new_points, cells, SubCellData());
+
+        // merge the triangulations - this will make sure that the duplicate
+        // vertices in the interior are absorbed
+        if (round == dim - 1)
+          merge_triangulations(tria_copy, rotated_tria, tria, 1e-12 * radius);
+        else
+          merge_triangulations(tria_copy,
+                               rotated_tria,
+                               tria_piece,
+                               1e-12 * radius);
+      }
+
+    for (const auto &cell : tria.cell_iterators())
+      if (cell->center().norm_square() > 0.4 * radius)
+        cell->set_manifold_id(1);
+      else
+        cell->set_all_manifold_ids(numbers::flat_manifold_id);
+
+    tria.set_all_manifold_ids_on_boundary(0);
+    tria.set_manifold(0, SphericalManifold<dim>(p));
+  }
+
+
+
   template <>
   void hyper_shell(Triangulation<3> & tria,
                    const Point<3> &   p,

diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in
index 6db7ae06c4fa7962afa4628f1e3f7a70e310faf9..b12ceb858b1050da2b9cde03241ed00f64c450e6 100644 (file)
--- a/source/grid/grid_generator.inst.in
+++ b/source/grid/grid_generator.inst.in
@@ -223,6 +223,12 @@ for (deal_II_dimension : DIMENSIONS)
                             const double,
                             const unsigned int);
 
+#if deal_II_dimension >= 2
+      template void
+      hyper_ball_balanced<deal_II_dimension>(Triangulation<deal_II_dimension> &,
+                                             const Point<deal_II_dimension> &,
+                                             const double);
+#endif
     \}
   }