From 93802238df352053f362e771d0b059aa0d90010b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 23 Jan 2025 15:48:37 -0700 Subject: [PATCH] Make GridGenerator::hyper_ball() work in dim=2, spacedim=3. --- include/deal.II/grid/grid_generator.h | 21 ++- source/grid/grid_generator.cc | 230 +++++++++++++------------- source/grid/grid_generator.inst.in | 8 +- 3 files changed, 132 insertions(+), 127 deletions(-) diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 90f9c9fa98..e574a82fbd 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -730,12 +730,12 @@ namespace GridGenerator /** * Initialize the given triangulation with several * @ref GlossCoarseMesh "coarse mesh cells" - * that cover a hyperball, i.e. a circle in 2d or a - * ball in 3d, around @p center with given @p radius. The function is + * that cover a hyperball, i.e. a circular disk if `dim==2` or a + * ball if `dim==3`, around @p center with given @p radius. The function is * used in step-6. * - * In order to avoid degenerate cells at the boundaries, the circle is - * triangulated by five cells, whereas in 3d the ball is subdivided by + * In order to avoid degenerate cells at the boundaries, the circular disk is + * triangulated by five cells, whereas in 3d the ball is subdivided into * seven cells. Specifically, these * cells are one cell in the center plus one "cap" cell on each of the faces * of this center cell. This ensures that under repeated refinement, none @@ -747,7 +747,10 @@ namespace GridGenerator * after one refinement is optimized. * * This function is declared to exist for triangulations of all space - * dimensions, but throws an error if called in 1d. + * dimensions, but throws an error if called if `dim==1`. If `spacedim>dim` + * (which is only possible if `dim==2, spacedim==3`) then the function + * creates a circular disk for which all vertices have a $z$ value equal + * to the $z$ coordinate of the center point provided by the second argument. * * By default, the manifold_id is set to 0 on the boundary faces, 1 on the * boundary cells, and numbers::flat_manifold_id on the central cell and on @@ -806,11 +809,11 @@ namespace GridGenerator * @note This function is available through the python interface as * `triangulation.generate_hyper_ball(point, radius = 1.)`. */ - template + template void - hyper_ball(Triangulation &tria, - const Point ¢er = Point(), - const double radius = 1., + hyper_ball(Triangulation &tria, + const Point ¢er = {}, + const double radius = 1., const bool attach_spherical_manifold_on_boundary_cells = false); /** diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index f15ecb478d..e1b25bc712 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -3929,15 +3929,6 @@ namespace GridGenerator - template <> - void - hyper_ball(Triangulation<1> &, const Point<1> &, const double, const bool) - { - DEAL_II_NOT_IMPLEMENTED(); - } - - - template <> void hyper_ball_balanced(Triangulation<1> &, const Point<1> &, const double) @@ -4293,48 +4284,125 @@ namespace GridGenerator - // Implementation for 2d only - template <> + template void - hyper_ball(Triangulation<2> &tria, - const Point<2> &p, - const double radius, - const bool internal_manifolds) + hyper_ball(Triangulation &tria, + const Point &p, + const double radius, + const bool internal_manifolds) { - // equilibrate cell sizes at - // transition from the inner part - // to the radial cells - const double a = 1. / (1 + std::sqrt(2.0)); - const Point<2> vertices[8] = { - p + Point<2>(-1, -1) * (radius / std::sqrt(2.0)), - p + Point<2>(+1, -1) * (radius / std::sqrt(2.0)), - p + Point<2>(-1, -1) * (radius / std::sqrt(2.0) * a), - p + Point<2>(+1, -1) * (radius / std::sqrt(2.0) * a), - p + Point<2>(-1, +1) * (radius / std::sqrt(2.0) * a), - p + Point<2>(+1, +1) * (radius / std::sqrt(2.0) * a), - p + Point<2>(-1, +1) * (radius / std::sqrt(2.0)), - p + Point<2>(+1, +1) * (radius / std::sqrt(2.0))}; + if constexpr (dim == 2) + { + const auto embed_point = [](const double x, + const double y) -> Point { + if constexpr (spacedim == 2) + return {x, y}; + else if constexpr (spacedim == 3) + return {x, y, 0}; + else + DEAL_II_NOT_IMPLEMENTED(); + }; - std::vector> cells(5, CellData<2>()); - for (unsigned int i = 0; i < 5; ++i) - { - for (unsigned int j = 0; j < 4; ++j) - cells[i].vertices[j] = circle_cell_vertices[i][j]; - cells[i].material_id = 0; - cells[i].manifold_id = i == 2 ? numbers::flat_manifold_id : 1; - } + // Equilibrate cell sizes at transition from the inner part + // to the radial cells + const double a = 1. / (1 + std::sqrt(2.0)); + const Point vertices[8] = { + p + embed_point(-1, -1) * (radius / std::sqrt(2.0)), + p + embed_point(+1, -1) * (radius / std::sqrt(2.0)), + p + embed_point(-1, -1) * (radius / std::sqrt(2.0) * a), + p + embed_point(+1, -1) * (radius / std::sqrt(2.0) * a), + p + embed_point(-1, +1) * (radius / std::sqrt(2.0) * a), + p + embed_point(+1, +1) * (radius / std::sqrt(2.0) * a), + p + embed_point(-1, +1) * (radius / std::sqrt(2.0)), + p + embed_point(+1, +1) * (radius / std::sqrt(2.0))}; - tria.create_triangulation(std::vector>(std::begin(vertices), - std::end(vertices)), - cells, - SubCellData()); // no boundary information - tria.set_all_manifold_ids_on_boundary(0); - tria.set_manifold(0, SphericalManifold<2>(p)); - if (internal_manifolds) - tria.set_manifold(1, SphericalManifold<2>(p)); + std::vector> cells(5, CellData<2>()); + + for (unsigned int i = 0; i < 5; ++i) + { + for (unsigned int j = 0; j < 4; ++j) + cells[i].vertices[j] = circle_cell_vertices[i][j]; + cells[i].material_id = 0; + cells[i].manifold_id = i == 2 ? numbers::flat_manifold_id : 1; + } + + tria.create_triangulation(std::vector>( + std::begin(vertices), std::end(vertices)), + cells, + SubCellData()); // no boundary information + tria.set_all_manifold_ids_on_boundary(0); + tria.set_manifold(0, SphericalManifold(p)); + if (internal_manifolds) + tria.set_manifold(1, SphericalManifold(p)); + else + tria.set_manifold(1, FlatManifold()); + } + else if constexpr (dim == 3) + { + const double a = + 1. / (1 + std::sqrt(3.0)); // equilibrate cell sizes at transition + // from the inner part to the radial + // cells + const unsigned int n_vertices = 16; + const Point<3> vertices[n_vertices] = { + // first the vertices of the inner + // cell + p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0) * a), + p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0) * a), + p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0) * a), + p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0) * a), + p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0) * a), + p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0) * a), + p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0) * a), + p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0) * a), + // now the eight vertices at + // the outer sphere + p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0)), + p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0)), + p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0)), + p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0)), + p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0)), + p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0)), + p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0)), + p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0)), + }; + + // one needs to draw the seven cubes to + // understand what's going on here + const unsigned int n_cells = 7; + const int cell_vertices[n_cells][8] = { + {0, 1, 4, 5, 3, 2, 7, 6}, // center + {8, 9, 12, 13, 0, 1, 4, 5}, // bottom + {9, 13, 1, 5, 10, 14, 2, 6}, // right + {11, 10, 3, 2, 15, 14, 7, 6}, // top + {8, 0, 12, 4, 11, 3, 15, 7}, // left + {8, 9, 0, 1, 11, 10, 3, 2}, // front + {12, 4, 13, 5, 15, 7, 14, 6}}; // back + + std::vector> cells(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; + cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1; + } + + tria.create_triangulation(std::vector>(std::begin(vertices), + std::end(vertices)), + cells, + SubCellData()); // no boundary information + tria.set_all_manifold_ids_on_boundary(0); + tria.set_manifold(0, SphericalManifold<3>(p)); + if (internal_manifolds) + tria.set_manifold(1, SphericalManifold<3>(p)); + else + tria.set_manifold(1, FlatManifold<3>()); + } else - tria.set_manifold(1, FlatManifold<2>()); + DEAL_II_NOT_IMPLEMENTED(); } @@ -5079,78 +5147,6 @@ namespace GridGenerator - // Implementation for 3d only - template <> - void - hyper_ball(Triangulation<3> &tria, - const Point<3> &p, - const double radius, - const bool internal_manifold) - { - const double a = - 1. / (1 + std::sqrt(3.0)); // equilibrate cell sizes at transition - // from the inner part to the radial - // cells - const unsigned int n_vertices = 16; - const Point<3> vertices[n_vertices] = { - // first the vertices of the inner - // cell - p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0) * a), - p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0) * a), - p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0) * a), - p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0) * a), - p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0) * a), - p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0) * a), - p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0) * a), - p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0) * a), - // now the eight vertices at - // the outer sphere - p + Point<3>(-1, -1, -1) * (radius / std::sqrt(3.0)), - p + Point<3>(+1, -1, -1) * (radius / std::sqrt(3.0)), - p + Point<3>(+1, -1, +1) * (radius / std::sqrt(3.0)), - p + Point<3>(-1, -1, +1) * (radius / std::sqrt(3.0)), - p + Point<3>(-1, +1, -1) * (radius / std::sqrt(3.0)), - p + Point<3>(+1, +1, -1) * (radius / std::sqrt(3.0)), - p + Point<3>(+1, +1, +1) * (radius / std::sqrt(3.0)), - p + Point<3>(-1, +1, +1) * (radius / std::sqrt(3.0)), - }; - - // one needs to draw the seven cubes to - // understand what's going on here - const unsigned int n_cells = 7; - const int cell_vertices[n_cells][8] = { - {0, 1, 4, 5, 3, 2, 7, 6}, // center - {8, 9, 12, 13, 0, 1, 4, 5}, // bottom - {9, 13, 1, 5, 10, 14, 2, 6}, // right - {11, 10, 3, 2, 15, 14, 7, 6}, // top - {8, 0, 12, 4, 11, 3, 15, 7}, // left - {8, 9, 0, 1, 11, 10, 3, 2}, // front - {12, 4, 13, 5, 15, 7, 14, 6}}; // back - - std::vector> cells(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; - cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1; - } - - tria.create_triangulation(std::vector>(std::begin(vertices), - std::end(vertices)), - cells, - SubCellData()); // no boundary information - tria.set_all_manifold_ids_on_boundary(0); - tria.set_manifold(0, SphericalManifold<3>(p)); - if (internal_manifold) - tria.set_manifold(1, SphericalManifold<3>(p)); - else - tria.set_manifold(1, FlatManifold<3>()); - } - - - void non_standard_orientation_mesh(Triangulation<2> &tria, const unsigned int n_rotate_middle_square) diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in index 75a7d06b57..189ce57a84 100644 --- a/source/grid/grid_generator.inst.in +++ b/source/grid/grid_generator.inst.in @@ -40,7 +40,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const double, const bool); - template void subdivided_hyper_rectangle( Triangulation &, @@ -72,6 +71,13 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const Point &, const std::vector &); + template void + hyper_ball( + Triangulation &, + const Point &, + const double, + const bool); + template void cheese( Triangulation &, -- 2.39.5