From: Martin Kronbichler Date: Thu, 2 Apr 2020 19:33:01 +0000 (+0200) Subject: Add alternative GridGenerator::hyper_ball with better balance of cells X-Git-Tag: v9.2.0-rc1~44^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4725ecee36eef2877aec16c7edcba398e015dd3e;p=dealii.git Add alternative GridGenerator::hyper_ball with better balance of cells --- diff --git a/doc/doxygen/images/hyper_ball_balanced_2d.png b/doc/doxygen/images/hyper_ball_balanced_2d.png new file mode 100644 index 0000000000..4d41bd679f Binary files /dev/null and b/doc/doxygen/images/hyper_ball_balanced_2d.png differ diff --git a/doc/doxygen/images/hyper_ball_balanced_3d.png b/doc/doxygen/images/hyper_ball_balanced_3d.png new file mode 100644 index 0000000000..8a5d9a443a Binary files /dev/null and b/doc/doxygen/images/hyper_ball_balanced_3d.png differ diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 0684b44033..5ce7c04e2b 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -738,6 +738,38 @@ namespace GridGenerator const double radius = 1., const bool attach_spherical_manifold_on_boundary_cells = false); + /** + * This is an alternative to hyper_ball with 12 cells in 2d and 32 cells in + * 3d, which provides a better balance between the size of the cells around + * the outer curved boundaries and the cell in the interior. The mesh is + * based on the cells used by GridGenerator::quarter_hyper_ball() with + * appropriate copies and rotations to fill the whole ball. + * + * The following pictures show the resulting mesh in 2D (left) and 3D: + * + * + * + * + * + *
+ * + * + * + *
+ * + * By default, the manifold_id is set to 0 on the boundary faces, 1 on the + * boundary cells, and types::flat_manifold_id on the central cell and on + * internal faces. + * + * @pre The triangulation passed as argument needs to be empty when calling + * this function. + */ + template + void + hyper_ball_balanced(Triangulation &tria, + const Point & center = Point(), + const double radius = 1.); + /** * Creates a hyper sphere, i.e., a surface of a ball in @p spacedim * dimensions. This function only exists for dim+1=spacedim in 2 and 3 space diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index c31e815115..c96b628a1d 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -5333,6 +5333,99 @@ namespace GridGenerator } + + template + void + hyper_ball_balanced(Triangulation &tria, + const Point & p, + const double radius) + { + // We create the ball by duplicating the information in each dimension at + // a time by appropriate rotations, starting from the quarter ball. The + // rotations make sure we do not generate inverted cells that would appear + // if we tried the slightly simpler approach to simply mirror the cells. + + Triangulation tria_piece; + GridGenerator::quarter_hyper_ball(tria_piece, p, radius); + + for (unsigned int round = 0; round < dim; ++round) + { + Triangulation tria_copy; + tria_copy.copy_triangulation(tria_piece); + tria_piece.clear(); + std::vector> new_points(tria_copy.n_vertices()); + if (round == 0) + for (unsigned int v = 0; v < tria_copy.n_vertices(); ++v) + { + // rotate by 90 degrees counterclockwise + new_points[v][0] = -tria_copy.get_vertices()[v][1]; + new_points[v][1] = tria_copy.get_vertices()[v][0]; + if (dim == 3) + new_points[v][2] = tria_copy.get_vertices()[v][2]; + } + else if (round == 1) + { + for (unsigned int v = 0; v < tria_copy.n_vertices(); ++v) + { + // rotate by 180 degrees along the xy plane + new_points[v][0] = -tria_copy.get_vertices()[v][0]; + new_points[v][1] = -tria_copy.get_vertices()[v][1]; + if (dim == 3) + new_points[v][2] = tria_copy.get_vertices()[v][2]; + } + } + else if (round == 2) + for (unsigned int v = 0; v < tria_copy.n_vertices(); ++v) + { + // rotate by 180 degrees along the xz plane + Assert(dim == 3, ExcInternalError()); + new_points[v][0] = -tria_copy.get_vertices()[v][0]; + new_points[v][1] = tria_copy.get_vertices()[v][1]; + new_points[v][2] = -tria_copy.get_vertices()[v][2]; + } + else + Assert(false, ExcInternalError()); + + + // the cell data is exactly the same as before + std::vector> cells; + cells.reserve(tria_copy.n_cells()); + for (const auto &cell : tria_copy.cell_iterators()) + { + CellData data; + for (unsigned int v : GeometryInfo::vertex_indices()) + data.vertices[v] = cell->vertex_index(v); + data.material_id = cell->material_id(); + data.manifold_id = cell->manifold_id(); + cells.push_back(data); + } + + Triangulation rotated_tria; + rotated_tria.create_triangulation(new_points, cells, SubCellData()); + + // merge the triangulations - this will make sure that the duplicate + // vertices in the interior are absorbed + if (round == dim - 1) + merge_triangulations(tria_copy, rotated_tria, tria, 1e-12 * radius); + else + merge_triangulations(tria_copy, + rotated_tria, + tria_piece, + 1e-12 * radius); + } + + for (const auto &cell : tria.cell_iterators()) + if (cell->center().norm_square() > 0.4 * radius) + cell->set_manifold_id(1); + else + cell->set_all_manifold_ids(numbers::flat_manifold_id); + + tria.set_all_manifold_ids_on_boundary(0); + tria.set_manifold(0, SphericalManifold(p)); + } + + + template <> void hyper_shell(Triangulation<3> & tria, const Point<3> & p, diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in index 6db7ae06c4..b12ceb858b 100644 --- a/source/grid/grid_generator.inst.in +++ b/source/grid/grid_generator.inst.in @@ -223,6 +223,12 @@ for (deal_II_dimension : DIMENSIONS) const double, const unsigned int); +#if deal_II_dimension >= 2 + template void + hyper_ball_balanced(Triangulation &, + const Point &, + const double); +#endif \} }