Implement additional options for 3d GridGenerator::hyper_shell

diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h
index 7210e43516af0ddaa6c920d94b4033599d53a115..c412b03205b06dc7bc6f7c733ca5940a88e6a725 100644 (file)
--- a/include/deal.II/grid/grid_generator.h
+++ b/include/deal.II/grid/grid_generator.h
@@ -1050,25 +1050,32 @@ namespace GridGenerator
    * is zero (as is the default), then it is computed adaptively such that the
    * resulting elements have the least aspect ratio.
    *
-   * In 3d, only certain numbers are allowed
+   * In 3d, only certain numbers are allowed:
    * <ul>
    * <li> 6 (or the default 0) for a surface based on a hexahedron (i.e. 6
    *      panels on the inner sphere extruded in radial direction to form 6
    *      cells),
    * <li> 12 for the rhombic dodecahedron,
-   * <li> 24 for the hexahedron-based surface refined once in the spherical
+   * <li> 24 for the hexahedron-based surface refined once in the azimuthal
    *      directions but not in the radial direction,
-   * <li> 48 for the rhombic dodecahedron refined once in the spherical
+   * <li> 48 for the rhombic dodecahedron refined once in the azimuthal
    *      directions but not in the radial direction,
    * <li> 96 for the rhombic dodecahedron refined once. This choice dates from
    *      an older version of deal.II before the Manifold classes were
    *      implemented: today this choce is equivalent to the rhombic
    *      dodecahedron after performing one global refinement.
+   * <li> Numbers of the kind $192\times 2^m$ with $m\geq 0$ integer. This
+   *      choice is similar to the 24 and 48 cell cases, but provides
+   *      additional refinements in azimuthal direction combined with a single
+   *      layer in radial direction. The base mesh is either the 6 or 12 cell
+   *      version, depending on whether $m$ in the power is odd or even,
+   *      respectively.
    * </ul>
-   * The versions with 24 and 48 cells are useful if the shell is thin and the
-   * radial lengths should be made more similar to the circumferential lengths.
+   * The versions with 24, 48, and $2^m 192$ cells are useful if the shell is
+   * thin and the radial lengths should be made more similar to the
+   * circumferential lengths.
    *
-   * The grids with 12 and 96 cells are plotted below:
+   * The 3d grids with 12 and 96 cells are plotted below:
    *
    * @image html hypershell3d-12.png
    * @image html hypershell3d-96.png

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 0d02fa7aad9a286ee086a09fdef43cc81aa25b12..2ddd40cab1766a5a21c35c12c295f44869d2eed7 100644 (file)
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -5344,7 +5344,22 @@ namespace GridGenerator
     Assert((inner_radius > 0) && (inner_radius < outer_radius),
            ExcInvalidRadii());
 
-    const unsigned int n = (n_cells == 0) ? 6 : n_cells;
+    unsigned int n_refinement_steps = 0;
+    unsigned int n_cells_coarsened  = n_cells;
+    if (n_cells != 96 && n_cells > 12)
+      while (n_cells_coarsened > 12 && n_cells_coarsened % 4 == 0)
+        {
+          ++n_refinement_steps;
+          n_cells_coarsened /= 4;
+        }
+    Assert(n_cells == 0 || n_cells == 6 || n_cells == 12 || n_cells == 96 ||
+             (n_refinement_steps > 0 &&
+              (n_cells_coarsened == 6 || n_cells_coarsened == 12)),
+           ExcMessage("Invalid number of coarse mesh cells"));
+
+    const unsigned int n = n_refinement_steps > 0 ?
+                             4 * n_cells_coarsened :
+                             ((n_cells == 0) ? 6 : n_cells);
 
     const double             irad = inner_radius / std::sqrt(3.0);
     const double             orad = outer_radius / std::sqrt(3.0);
@@ -5447,19 +5462,60 @@ namespace GridGenerator
         case 48:
           {
             // These two meshes are created by first creating a mesh of the
-            // 6-cell/12-cell version, refining globally once, and finally
-            // removing the outer half of the cells
-            Triangulation<3> tmp;
-            hyper_shell(
-              tmp, p, inner_radius, 2 * outer_radius - inner_radius, n / 4);
-            tmp.refine_global(1);
-            std::set<Triangulation<3>::active_cell_iterator> cells_to_remove;
-            for (const auto &cell : tmp.active_cell_iterators())
-              if (cell->center(true).norm_square() >
-                  outer_radius * outer_radius)
-                cells_to_remove.insert(cell);
-            AssertDimension(cells_to_remove.size(), n);
-            create_triangulation_with_removed_cells(tmp, cells_to_remove, tria);
+            // 6-cell/12-cell version, refining globally, and removing the
+            // outer half of the cells. For 192 and more cells, we do this
+            // iteratively several times, always refining and removing the
+            // outer half. Thus, the outer radius for the start is larger and
+            // set as 2^n_refinement_steps such that it exactly gives the
+            // desired radius in the end. It would have been slightly less
+            // code to treat refinement steps recursively for 192 cells or
+            // beyond, but unfortunately we could end up with the 96 cell case
+            // which is not what we want. Thus, we need to implement a loop
+            // manually here.
+            Triangulation<3>   tmp;
+            const unsigned int outer_radius_factor = 1 << n_refinement_steps;
+            hyper_shell(tmp,
+                        p,
+                        inner_radius,
+                        outer_radius_factor * outer_radius -
+                          (outer_radius_factor - 1) * inner_radius,
+                        n / 4);
+            for (unsigned int r = 0; r < n_refinement_steps; ++r)
+              {
+                tmp.refine_global(1);
+                std::set<Triangulation<3>::active_cell_iterator>
+                  cells_to_remove;
+
+                // We remove all cells which do not have exactly four vertices
+                // at the inner radius (plus some tolerance).
+                for (const auto &cell : tmp.active_cell_iterators())
+                  {
+                    unsigned int n_vertices_inside = 0;
+                    for (const auto v : GeometryInfo<3>::vertex_indices())
+                      if ((cell->vertex(v) - p).norm_square() <
+                          inner_radius * inner_radius * (1 + 1e-12))
+                        ++n_vertices_inside;
+                    if (n_vertices_inside < 4)
+                      cells_to_remove.insert(cell);
+                  }
+
+                AssertDimension(cells_to_remove.size(),
+                                tmp.n_active_cells() / 2);
+                if (r == n_refinement_steps - 1)
+                  create_triangulation_with_removed_cells(tmp,
+                                                          cells_to_remove,
+                                                          tria);
+                else
+                  {
+                    Triangulation<3> copy;
+                    create_triangulation_with_removed_cells(tmp,
+                                                            cells_to_remove,
+                                                            copy);
+                    tmp = std::move(copy);
+                    tmp.set_all_manifold_ids(0);
+                    tmp.set_manifold(0, SphericalManifold<3>(p));
+                  }
+              }
             break;
           }
         case 96:
@@ -5480,6 +5536,9 @@ namespace GridGenerator
           }
       }
 
+    if (n_cells > 0)
+      AssertDimension(tria.n_global_active_cells(), n_cells);
+
     if (colorize)
       colorize_hyper_shell(tria, p, inner_radius, outer_radius);
     tria.set_all_manifold_ids(0);