]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix compilation with clang-13. 14009/head
authorDavid Wells <drwells@email.unc.edu>
Fri, 17 Jun 2022 16:46:30 +0000 (12:46 -0400)
committerDavid Wells <drwells@email.unc.edu>
Fri, 17 Jun 2022 21:34:53 +0000 (17:34 -0400)
source/grid/grid_generator.cc

index d1e491dcc6c3d1a50e677c84c4589d218462e305..ebaaae4f0be0c86a398fc5c72b6dd4fbae2db268 100644 (file)
@@ -5927,6 +5927,212 @@ namespace GridGenerator
     tria.set_manifold(0, SphericalManifold<dim>(p));
   }
 
+  // To work around an internal clang-13 error we need to split up the
+  // individual hyper shell functions. This has the added bonus of making the
+  // control flow easier to follow - some hyper shell functions call others.
+  namespace internal
+  {
+    namespace
+    {
+      void
+      hyper_shell_6(Triangulation<3> &tria,
+                    const Point<3> &  p,
+                    const double      inner_radius,
+                    const double      outer_radius)
+      {
+        std::vector<Point<3>>    vertices;
+        std::vector<CellData<3>> cells;
+
+        const double irad = inner_radius / std::sqrt(3.0);
+        const double orad = outer_radius / std::sqrt(3.0);
+
+        // Corner points of the cube [-1,1]^3
+        static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1}, //
+                                                            {+1, -1, -1}, //
+                                                            {-1, +1, -1}, //
+                                                            {+1, +1, -1}, //
+                                                            {-1, -1, +1}, //
+                                                            {+1, -1, +1}, //
+                                                            {-1, +1, +1}, //
+                                                            {+1, +1, +1}}};
+
+        // 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;
+          }
+
+        tria.create_triangulation(vertices, cells, SubCellData());
+        tria.set_all_manifold_ids(0);
+        tria.set_manifold(0, SphericalManifold<3>(p));
+      }
+
+      void
+      hyper_shell_12(Triangulation<3> &tria,
+                     const Point<3> &  p,
+                     const double      inner_radius,
+                     const double      outer_radius)
+      {
+        std::vector<Point<3>>    vertices;
+        std::vector<CellData<3>> cells;
+
+        const double irad = inner_radius / std::sqrt(3.0);
+        const double orad = outer_radius / std::sqrt(3.0);
+
+        // 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}}};
+
+        // Corner points of the cube [-1,1]^3
+        static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1}, //
+                                                            {+1, -1, -1}, //
+                                                            {-1, +1, -1}, //
+                                                            {+1, +1, -1}, //
+                                                            {-1, -1, +1}, //
+                                                            {+1, -1, +1}, //
+                                                            {-1, +1, +1}, //
+                                                            {+1, +1, +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)
+          {
+            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;
+          }
+
+        tria.create_triangulation(vertices, cells, SubCellData());
+        tria.set_all_manifold_ids(0);
+        tria.set_manifold(0, SphericalManifold<3>(p));
+      }
+
+      void
+      hyper_shell_24_48(Triangulation<3> & tria,
+                        const unsigned int n,
+                        const unsigned int n_refinement_steps,
+                        const Point<3> &   p,
+                        const double       inner_radius,
+                        const double       outer_radius)
+      {
+        // These two meshes are created by first creating a mesh of the
+        // 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;
+        if (n == 24)
+          hyper_shell_6(tmp,
+                        p,
+                        inner_radius,
+                        outer_radius_factor * outer_radius -
+                          (outer_radius_factor - 1) * inner_radius);
+        else if (n == 48)
+          hyper_shell_12(tmp,
+                         p,
+                         inner_radius,
+                         outer_radius_factor * outer_radius -
+                           (outer_radius_factor - 1) * inner_radius);
+        else
+          Assert(n == 24 || n == 48, ExcInternalError());
+        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));
+              }
+          }
+        tria.set_all_manifold_ids(0);
+        tria.set_manifold(0, SphericalManifold<3>(p));
+      }
+
+    } // namespace
+  }   // namespace internal
+
 
 
   template <>
@@ -5958,163 +6164,19 @@ namespace GridGenerator
                              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);
-    std::vector<Point<3>>    vertices;
-    std::vector<CellData<3>> cells;
-
-    // Corner points of the cube [-1,1]^3
-    static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1}, //
-                                                        {+1, -1, -1}, //
-                                                        {-1, +1, -1}, //
-                                                        {+1, +1, -1}, //
-                                                        {-1, -1, +1}, //
-                                                        {+1, -1, +1}, //
-                                                        {-1, +1, +1}, //
-                                                        {+1, +1, +1}}};
-
     switch (n)
       {
         case 6:
-          {
-            // 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;
-              }
-
-            tria.create_triangulation(vertices, cells, SubCellData());
-            break;
-          }
+          internal::hyper_shell_6(tria, p, inner_radius, outer_radius);
+          break;
         case 12:
-          {
-            // 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)
-              {
-                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;
-              }
-
-            tria.create_triangulation(vertices, cells, SubCellData());
-            break;
-          }
+          internal::hyper_shell_12(tria, p, inner_radius, outer_radius);
+          break;
         case 24:
         case 48:
-          {
-            // These two meshes are created by first creating a mesh of the
-            // 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;
-          }
+          internal::hyper_shell_24_48(
+            tria, n, n_refinement_steps, p, inner_radius, outer_radius);
+          break;
         case 96:
           {
             // create a triangulation based on the 12-cell version. This
@@ -6122,9 +6184,11 @@ namespace GridGenerator
             // 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);
+            internal::hyper_shell_12(tmp, p, inner_radius, outer_radius);
             tmp.refine_global(1);
             flatten_triangulation(tmp, tria);
+            tria.set_all_manifold_ids(0);
+            tria.set_manifold(0, SphericalManifold<3>(p));
             break;
           }
         default:
@@ -6138,8 +6202,6 @@ namespace GridGenerator
 
     if (colorize)
       colorize_hyper_shell(tria, p, inner_radius, outer_radius);
-    tria.set_all_manifold_ids(0);
-    tria.set_manifold(0, SphericalManifold<3>(p));
   }
 
 

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.