From 3b3e75c77f71399c9cb6e9794aaf819177052520 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 18 Jun 2025 11:57:14 -0600 Subject: [PATCH] Fix quarter hyper shell colorization for thin shells. --- doc/news/changes/minor/20250618Gassmoeller | 8 ++ source/grid/grid_generator.cc | 108 +++++++++--------- tests/grid/grid_generator_hyper_shell.cc | 110 +++++++++++++++++++ tests/grid/grid_generator_hyper_shell.output | 21 ++++ 4 files changed, 189 insertions(+), 58 deletions(-) create mode 100644 doc/news/changes/minor/20250618Gassmoeller create mode 100644 tests/grid/grid_generator_hyper_shell.cc create mode 100644 tests/grid/grid_generator_hyper_shell.output diff --git a/doc/news/changes/minor/20250618Gassmoeller b/doc/news/changes/minor/20250618Gassmoeller new file mode 100644 index 0000000000..8539a662ac --- /dev/null +++ b/doc/news/changes/minor/20250618Gassmoeller @@ -0,0 +1,8 @@ +Fixed: The function GridGenerator::quarter_hyper_shell would colorize +the boundary ids of the triangulation incorrectly if the created shell +was too thin. In particular the top and bottom boundary ids would +not be set. This did not matter for the bottom boundary since it is 0, +which is also the default value, but the top boundary would not be +marked correctly for these geometries. This is fixed now. +
+(Rene Gassmoeller, 2025/06/18) diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 3514905d2d..2d244b404b 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -1454,74 +1454,66 @@ namespace GridGenerator if (tria.n_cells() != 3) AssertThrow(false, ExcNotImplemented()); - double middle = (outer_radius - inner_radius) / 2e0 + inner_radius; - double eps = 1e-3 * middle; - Triangulation<3>::cell_iterator cell = tria.begin(); + const double middle_radius = + (outer_radius - inner_radius) / 2e0 + inner_radius; + const double eps = 1e-3 * middle_radius; - for (; cell != tria.end(); ++cell) - for (const unsigned int f : GeometryInfo<3>::face_indices()) + for (const auto &cell : tria.cell_iterators()) + for (const unsigned int f : cell->face_indices()) { - if (!cell->face(f)->at_boundary()) + const auto face = cell->face(f); + if (!face->at_boundary()) continue; - double radius = cell->face(f)->center().norm() - center.norm(); - if (std::fabs(cell->face(f)->center()[0]) < - eps) // x = 0 set boundary 2 + const double face_center_radius = + (face->center(true) - center).norm(); + + if (std::fabs(face->center()[0]) < eps) // x = 0 set boundary 2 { - cell->face(f)->set_boundary_id(2); - for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; - ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - - cell->face(f)->line(j)->vertex(1).norm()) > - eps) - cell->face(f)->line(j)->set_boundary_id(2); + face->set_boundary_id(2); + for (const unsigned int line_no : face->line_indices()) + if (face->line(line_no)->at_boundary()) + if (std::fabs(face->line(line_no)->vertex(0).norm() - + face->line(line_no)->vertex(1).norm()) > eps) + face->line(line_no)->set_boundary_id(2); } - else if (std::fabs(cell->face(f)->center()[1]) < - eps) // y = 0 set boundary 3 + else if (std::fabs(face->center()[1]) < eps) // y = 0 set boundary 3 { - cell->face(f)->set_boundary_id(3); - for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; - ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - - cell->face(f)->line(j)->vertex(1).norm()) > - eps) - cell->face(f)->line(j)->set_boundary_id(3); + face->set_boundary_id(3); + for (const unsigned int line_no : face->line_indices()) + if (face->line(line_no)->at_boundary()) + if (std::fabs(face->line(line_no)->vertex(0).norm() - + face->line(line_no)->vertex(1).norm()) > eps) + face->line(line_no)->set_boundary_id(3); } - else if (std::fabs(cell->face(f)->center()[2]) < - eps) // z = 0 set boundary 4 + else if (std::fabs(face->center()[2]) < eps) // z = 0 set boundary 4 { - cell->face(f)->set_boundary_id(4); - for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; - ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - - cell->face(f)->line(j)->vertex(1).norm()) > - eps) - cell->face(f)->line(j)->set_boundary_id(4); + face->set_boundary_id(4); + for (const unsigned int line_no : face->line_indices()) + if (face->line(line_no)->at_boundary()) + if (std::fabs(face->line(line_no)->vertex(0).norm() - + face->line(line_no)->vertex(1).norm()) > eps) + face->line(line_no)->set_boundary_id(4); } - else if (radius < middle) // inner radius set boundary 0 + else if (face_center_radius < + middle_radius) // inner radius set boundary 0 { - cell->face(f)->set_boundary_id(0); - for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; - ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - - cell->face(f)->line(j)->vertex(1).norm()) < - eps) - cell->face(f)->line(j)->set_boundary_id(0); + face->set_boundary_id(0); + for (const unsigned int line_no : face->line_indices()) + if (face->line(line_no)->at_boundary()) + if (std::fabs(face->line(line_no)->vertex(0).norm() - + face->line(line_no)->vertex(1).norm()) < eps) + face->line(line_no)->set_boundary_id(0); } - else if (radius > middle) // outer radius set boundary 1 + else if (face_center_radius > + middle_radius) // outer radius set boundary 1 { - cell->face(f)->set_boundary_id(1); - for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; - ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - - cell->face(f)->line(j)->vertex(1).norm()) < - eps) - cell->face(f)->line(j)->set_boundary_id(1); + face->set_boundary_id(1); + for (const unsigned int line_no : face->line_indices()) + if (face->line(line_no)->at_boundary()) + if (std::fabs(face->line(line_no)->vertex(0).norm() - + face->line(line_no)->vertex(1).norm()) < eps) + face->line(line_no)->set_boundary_id(1); } else DEAL_II_ASSERT_UNREACHABLE(); @@ -6408,11 +6400,11 @@ namespace GridGenerator AssertThrow(false, ExcNotImplemented()); } - if (colorize) - colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius); - tria.set_all_manifold_ids(0); tria.set_manifold(0, SphericalManifold<3>(center)); + + if (colorize) + colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius); } diff --git a/tests/grid/grid_generator_hyper_shell.cc b/tests/grid/grid_generator_hyper_shell.cc new file mode 100644 index 0000000000..df30dfe431 --- /dev/null +++ b/tests/grid/grid_generator_hyper_shell.cc @@ -0,0 +1,110 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2013 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// Test if colorizing a thin quarter_hyper_shell yields correct +// boundary indicators. The function used to not colorize +// boundaries correctly if the hyper shell was too thin. + +#include + +#include +#include +#include +#include + +#include "../tests.h" + + + +template +void +test(const Point origin = Point()) +{ + SphericalManifold manifold(origin); + { + deallog << "quarter_hyper_shell with origin: " << origin << std::endl; + Triangulation tr; + GridGenerator::quarter_hyper_shell(tr, origin, 1460800, 1560800, 0, true); + tr.set_manifold(0, manifold); + + bool found_top_boundary = false; + const types::boundary_id top_id = 1; + deallog << std::endl << "Boundary ids of boundary faces:"; + + for (const auto &cell : tr.active_cell_iterators()) + for (const unsigned int face_no : cell->face_indices()) + { + const typename Triangulation::face_iterator face = + cell->face(face_no); + + if (face->at_boundary()) + { + const types::boundary_id boundary_id = face->boundary_id(); + deallog << " " << boundary_id; + + if (boundary_id == top_id) + { + found_top_boundary = true; + } + } + } + + if (dim > 2) + { + deallog << std::endl << std::endl << "Boundary ids of boundary lines:"; + + for (const auto &cell : tr.active_cell_iterators()) + for (const unsigned int face_no : cell->face_indices()) + { + const typename Triangulation::face_iterator face = + cell->face(face_no); + + if (face->at_boundary()) + { + for (const unsigned l : face->line_indices()) + deallog << " " << face->line(l)->boundary_id(); + } + } + } + + deallog << std::endl; + + if (found_top_boundary) + deallog << "Ok. Found top boundary." << std::endl; + else + deallog << "Error. Failed to find top boundary." << std::endl; + } +} + + +int +main() +{ + initlog(); + + deallog.get_file_stream() << std::setprecision(9); + + deallog.push("2d"); + test<2>(); + test<2>(Point<2>(1.0, 0.0)); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + test<3>(Point<3>(1.0, 0.0, 0.0)); + + deallog.pop(); +} diff --git a/tests/grid/grid_generator_hyper_shell.output b/tests/grid/grid_generator_hyper_shell.output new file mode 100644 index 0000000000..a8d86594fd --- /dev/null +++ b/tests/grid/grid_generator_hyper_shell.output @@ -0,0 +1,21 @@ + +DEAL:2d::quarter_hyper_shell with origin: 0.00000 0.00000 +DEAL:2d:: +DEAL:2d::Boundary ids of boundary faces: 3 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 2 1 0 +DEAL:2d::Ok. Found top boundary. +DEAL:2d::quarter_hyper_shell with origin: 1.00000 0.00000 +DEAL:2d:: +DEAL:2d::Boundary ids of boundary faces: 3 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 2 1 0 +DEAL:2d::Ok. Found top boundary. +DEAL:3d::quarter_hyper_shell with origin: 0.00000 0.00000 0.00000 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary faces: 2 0 1 4 3 0 1 4 2 0 1 3 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary lines: 0 1 4 2 0 0 0 0 1 1 1 1 4 4 0 1 0 1 4 3 0 0 0 0 1 1 1 1 4 4 0 1 0 1 2 3 0 0 0 0 1 1 1 1 3 3 0 1 +DEAL:3d::Ok. Found top boundary. +DEAL:3d::quarter_hyper_shell with origin: 1.00000 0.00000 0.00000 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary faces: 2 0 1 4 3 0 1 4 2 0 1 3 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary lines: 0 1 4 2 0 0 0 0 1 1 1 1 4 4 0 1 0 1 4 3 0 0 0 0 1 1 1 1 4 4 0 1 0 1 2 3 0 0 0 0 1 1 1 1 3 3 0 1 +DEAL:3d::Ok. Found top boundary. -- 2.39.5