From: Martin Kronbichler Date: Tue, 7 Apr 2020 16:19:09 +0000 (+0200) Subject: Implement additional options for 3d GridGenerator::hyper_shell X-Git-Tag: v9.2.0-rc1~19^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bfa024ac3ad9a7344423571430bc5ba3f105b5b9;p=dealii.git Implement additional options for 3d GridGenerator::hyper_shell --- diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 7210e43516..c412b03205 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -1050,25 +1050,32 @@ 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 + * In 3d, only certain numbers are allowed: * - * 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 versions with 24, 48, and $2^m 192$ 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: + * The 3d grids with 12 and 96 cells are plotted below: * * @image html hypershell3d-12.png * @image html hypershell3d-96.png diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 0d02fa7aad..2ddd40cab1 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -5344,7 +5344,22 @@ namespace GridGenerator Assert((inner_radius > 0) && (inner_radius < outer_radius), ExcInvalidRadii()); - const unsigned int n = (n_cells == 0) ? 6 : n_cells; + unsigned int n_refinement_steps = 0; + unsigned int n_cells_coarsened = n_cells; + if (n_cells != 96 && n_cells > 12) + while (n_cells_coarsened > 12 && n_cells_coarsened % 4 == 0) + { + ++n_refinement_steps; + n_cells_coarsened /= 4; + } + Assert(n_cells == 0 || n_cells == 6 || n_cells == 12 || n_cells == 96 || + (n_refinement_steps > 0 && + (n_cells_coarsened == 6 || n_cells_coarsened == 12)), + ExcMessage("Invalid number of coarse mesh cells")); + + const unsigned int n = n_refinement_steps > 0 ? + 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); @@ -5447,19 +5462,60 @@ namespace GridGenerator case 48: { // 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::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); + // 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; } case 96: @@ -5480,6 +5536,9 @@ namespace GridGenerator } } + if (n_cells > 0) + AssertDimension(tria.n_global_active_cells(), n_cells); + if (colorize) colorize_hyper_shell(tria, p, inner_radius, outer_radius); tria.set_all_manifold_ids(0);