Add additional cell options for 3D hyper shell

diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h
index 0684b4403395d7161abbd9c71b3174b75fc8d6dc..7210e43516af0ddaa6c920d94b4033599d53a115 100644 (file)
--- a/include/deal.II/grid/grid_generator.h
+++ b/include/deal.II/grid/grid_generator.h
@@ -1050,12 +1050,23 @@ 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, 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), 12 for the rhombic dodecahedron,
-   * and 96. This choice dates from an older version of deal.II before the
-   * Manifold classes were implemented: today all three choices are roughly
-   * equivalent (after performing global refinement, of course).
+   * 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
+   *      directions but not in the radial direction,
+   * <li> 48 for the rhombic dodecahedron refined once in the spherical
+   *      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.
+   * </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 grids with 12 and 96 cells are plotted below:
    *

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index c31e8151157670b65cd38cde00d14098ad2ea90c..0d02fa7aad9a286ee086a09fdef43cc81aa25b12 100644 (file)
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -5361,111 +5361,123 @@ namespace GridGenerator
                                                         {-1, +1, +1}, //
                                                         {+1, +1, +1}}};
 
-    // Start with the shell bounded by
-    // two nested cubes
-    if (n == 6)
+    switch (n)
       {
-        for (unsigned int i = 0; i < 8; ++i)
-          vertices.push_back(p + hexahedron[i] * irad);
-        for (unsigned int i = 0; i < 8; ++i)
-          vertices.push_back(p + hexahedron[i] * orad);
-
-        const unsigned int n_cells                   = 6;
-        const int          cell_vertices[n_cells][8] = {
-          {8, 9, 10, 11, 0, 1, 2, 3},    // bottom
-          {9, 11, 1, 3, 13, 15, 5, 7},   // right
-          {12, 13, 4, 5, 14, 15, 6, 7},  // top
-          {8, 0, 10, 2, 12, 4, 14, 6},   // left
-          {8, 9, 0, 1, 12, 13, 4, 5},    // front
-          {10, 2, 11, 3, 14, 6, 15, 7}}; // back
-
-        cells.resize(n_cells, CellData<3>());
-
-        for (unsigned int i = 0; i < n_cells; ++i)
+        case 6:
           {
-            for (const unsigned int j : GeometryInfo<3>::vertex_indices())
-              cells[i].vertices[j] = cell_vertices[i][j];
-            cells[i].material_id = 0;
-          }
-
-        tria.create_triangulation(vertices, cells, SubCellData());
-      }
-    // A more regular subdivision can
-    // be obtained by two nested
-    // rhombic dodecahedra
+            // Start with the shell bounded by two nested cubes
+            for (unsigned int i = 0; i < 8; ++i)
+              vertices.push_back(p + hexahedron[i] * irad);
+            for (unsigned int i = 0; i < 8; ++i)
+              vertices.push_back(p + hexahedron[i] * orad);
+
+            const unsigned int n_cells                   = 6;
+            const int          cell_vertices[n_cells][8] = {
+              {8, 9, 10, 11, 0, 1, 2, 3},    // bottom
+              {9, 11, 1, 3, 13, 15, 5, 7},   // right
+              {12, 13, 4, 5, 14, 15, 6, 7},  // top
+              {8, 0, 10, 2, 12, 4, 14, 6},   // left
+              {8, 9, 0, 1, 12, 13, 4, 5},    // front
+              {10, 2, 11, 3, 14, 6, 15, 7}}; // back
+
+            cells.resize(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;
+              }
 
-    else if (n == 12)
-      {
-        // Octahedron inscribed in the cube [-1,1]^3
-        static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0}, //
-                                                            {1, 0, 0},  //
-                                                            {0, -1, 0}, //
-                                                            {0, 1, 0},  //
-                                                            {0, 0, -1}, //
-                                                            {0, 0, 1}}};
-
-        for (unsigned int i = 0; i < 8; ++i)
-          vertices.push_back(p + hexahedron[i] * irad);
-        for (unsigned int i = 0; i < 6; ++i)
-          vertices.push_back(p + octahedron[i] * inner_radius);
-        for (unsigned int i = 0; i < 8; ++i)
-          vertices.push_back(p + hexahedron[i] * orad);
-        for (unsigned int i = 0; i < 6; ++i)
-          vertices.push_back(p + octahedron[i] * outer_radius);
-
-        const unsigned int n_cells            = 12;
-        const unsigned int rhombi[n_cells][4] = {{10, 4, 0, 8},
-                                                 {4, 13, 8, 6},
-                                                 {10, 5, 4, 13},
-                                                 {1, 9, 10, 5},
-                                                 {9, 7, 5, 13},
-                                                 {7, 11, 13, 6},
-                                                 {9, 3, 7, 11},
-                                                 {1, 12, 9, 3},
-                                                 {12, 2, 3, 11},
-                                                 {2, 8, 11, 6},
-                                                 {12, 0, 2, 8},
-                                                 {1, 10, 12, 0}};
-
-        cells.resize(n_cells, CellData<3>());
-
-        for (unsigned int i = 0; i < n_cells; ++i)
+            tria.create_triangulation(vertices, cells, SubCellData());
+            break;
+          }
+        case 12:
           {
-            for (unsigned int j = 0; j < 4; ++j)
+            // A more regular subdivision can be obtained by two nested rhombic
+            // dodecahedra
+            //
+            // Octahedron inscribed in the cube [-1,1]^3
+            static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0}, //
+                                                                {1, 0, 0},  //
+                                                                {0, -1, 0}, //
+                                                                {0, 1, 0},  //
+                                                                {0, 0, -1}, //
+                                                                {0, 0, 1}}};
+
+            for (unsigned int i = 0; i < 8; ++i)
+              vertices.push_back(p + hexahedron[i] * irad);
+            for (unsigned int i = 0; i < 6; ++i)
+              vertices.push_back(p + octahedron[i] * inner_radius);
+            for (unsigned int i = 0; i < 8; ++i)
+              vertices.push_back(p + hexahedron[i] * orad);
+            for (unsigned int i = 0; i < 6; ++i)
+              vertices.push_back(p + octahedron[i] * outer_radius);
+
+            const unsigned int n_cells            = 12;
+            const unsigned int rhombi[n_cells][4] = {{10, 4, 0, 8},
+                                                     {4, 13, 8, 6},
+                                                     {10, 5, 4, 13},
+                                                     {1, 9, 10, 5},
+                                                     {9, 7, 5, 13},
+                                                     {7, 11, 13, 6},
+                                                     {9, 3, 7, 11},
+                                                     {1, 12, 9, 3},
+                                                     {12, 2, 3, 11},
+                                                     {2, 8, 11, 6},
+                                                     {12, 0, 2, 8},
+                                                     {1, 10, 12, 0}};
+
+            cells.resize(n_cells, CellData<3>());
+
+            for (unsigned int i = 0; i < n_cells; ++i)
               {
-                cells[i].vertices[j]     = rhombi[i][j];
-                cells[i].vertices[j + 4] = rhombi[i][j] + 14;
+                for (unsigned int j = 0; j < 4; ++j)
+                  {
+                    cells[i].vertices[j]     = rhombi[i][j];
+                    cells[i].vertices[j + 4] = rhombi[i][j] + 14;
+                  }
+                cells[i].material_id = 0;
               }
-            cells[i].material_id = 0;
-          }
 
-        tria.create_triangulation(vertices, cells, SubCellData());
-      }
-    else if (n == 96)
-      {
-        // create a triangulation based on the 12-cell version. This function
-        // was needed before SphericalManifold was written: it manually
-        // adjusted the interior vertices to lie along concentric
-        // spheres. Nowadays we can just refine globally:
-        Triangulation<3> tmp;
-        hyper_shell(tmp, p, inner_radius, outer_radius, 12);
-        tmp.refine_global(1);
-
-        // now copy the resulting level 1 cells into the new triangulation,
-        cells.resize(tmp.n_active_cells(), CellData<3>());
-        for (const auto &cell : tmp.active_cell_iterators())
+            tria.create_triangulation(vertices, cells, SubCellData());
+            break;
+          }
+        case 24:
+        case 48:
           {
-            const unsigned int cell_index = cell->active_cell_index();
-            for (const unsigned int v : GeometryInfo<3>::vertex_indices())
-              cells[cell_index].vertices[v] = cell->vertex_index(v);
-            cells[cell_index].material_id = 0;
+            // 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);
+            break;
+          }
+        case 96:
+          {
+            // create a triangulation based on the 12-cell version. This
+            // function was needed before SphericalManifold was written: it
+            // manually adjusted the interior vertices to lie along concentric
+            // spheres. Nowadays we can just refine globally:
+            Triangulation<3> tmp;
+            hyper_shell(tmp, p, inner_radius, outer_radius, 12);
+            tmp.refine_global(1);
+            flatten_triangulation(tmp, tria);
+            break;
+          }
+        default:
+          {
+            Assert(false, ExcMessage("Invalid number of coarse mesh cells."));
           }
-
-        tria.create_triangulation(tmp.get_vertices(), cells, SubCellData());
-      }
-    else
-      {
-        Assert(false, ExcMessage("Invalid number of coarse mesh cells."));
       }
 
     if (colorize)