From: David Wells Date: Fri, 17 Jun 2022 16:46:30 +0000 (-0400) Subject: Fix compilation with clang-13. X-Git-Tag: v9.4.0-rc2~5^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5ca98472f86aaf9704a2809929249266c1eb504a;p=dealii.git Fix compilation with clang-13. --- diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index d1e491dcc6..ebaaae4f0b 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -5927,6 +5927,212 @@ namespace GridGenerator tria.set_manifold(0, SphericalManifold(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> vertices; + std::vector> 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, 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> vertices; + std::vector> 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, 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, 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::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> vertices; - std::vector> cells; - - // Corner points of the cube [-1,1]^3 - static const std::array, 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, 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::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)); }