From b09287ccd19b4d44883fae964833729fd94637e8 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Thu, 19 Jun 2025 17:15:51 -0600 Subject: [PATCH] Correctly use manifold in half hyper shell colorization --- doc/news/changes/minor/20250619Gassmoeller | 6 + source/grid/grid_generator.cc | 28 ++--- .../grid_generator_half_hyper_shell.output | 8 +- .../grid_generator_half_hyper_shell_02.cc | 110 ++++++++++++++++++ .../grid_generator_half_hyper_shell_02.output | 21 ++++ 5 files changed, 155 insertions(+), 18 deletions(-) create mode 100644 doc/news/changes/minor/20250619Gassmoeller create mode 100644 tests/grid/grid_generator_half_hyper_shell_02.cc create mode 100644 tests/grid/grid_generator_half_hyper_shell_02.output diff --git a/doc/news/changes/minor/20250619Gassmoeller b/doc/news/changes/minor/20250619Gassmoeller new file mode 100644 index 0000000000..ef3e029bc4 --- /dev/null +++ b/doc/news/changes/minor/20250619Gassmoeller @@ -0,0 +1,6 @@ +Fixed: The function GridGenerator::half_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 +be set incorrectly. This is fixed now. +
+(Rene Gassmoeller, 2025/06/19) diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 3514905d2d..c061b92e06 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -6309,28 +6309,30 @@ namespace GridGenerator cells, SubCellData()); // no boundary information + tria.set_all_manifold_ids(0); + tria.set_manifold(0, SphericalManifold<3>(center)); + if (colorize) { - // We want to use a standard boundary description where - // the boundary is not curved. Hence set boundary id 2 to + // We want to use a flat boundary description where + // the boundary is not curved. Hence set boundary id 2 // to all faces in a first step. - Triangulation<3>::cell_iterator cell = tria.begin(); - for (; cell != tria.end(); ++cell) - for (const unsigned int i : GeometryInfo<3>::face_indices()) - if (cell->at_boundary(i)) - cell->face(i)->set_all_boundary_ids(2); + for (const auto &cell : tria.cell_iterators()) + for (const unsigned int f : cell->face_indices()) + if (cell->at_boundary(f)) + cell->face(f)->set_all_boundary_ids(2); // Next look for the curved boundaries. If the x value of the // center of the face is not equal to center(0), we're on a curved // boundary. Then decide whether the center is nearer to the inner // or outer boundary to set the correct boundary id. - for (cell = tria.begin(); cell != tria.end(); ++cell) - for (const unsigned int i : GeometryInfo<3>::face_indices()) - if (cell->at_boundary(i)) + for (const auto &cell : tria.cell_iterators()) + for (const unsigned int f : cell->face_indices()) + if (cell->at_boundary(f)) { - const Triangulation<3>::face_iterator face = cell->face(i); + const Triangulation<3>::face_iterator face = cell->face(f); + const Point<3> face_center(face->center(true)); - const Point<3> face_center(face->center()); if (std::abs(face_center[0] - center[0]) > 1.e-6 * face_center.norm()) { @@ -6342,8 +6344,6 @@ namespace GridGenerator } } } - tria.set_all_manifold_ids(0); - tria.set_manifold(0, SphericalManifold<3>(center)); } diff --git a/tests/grid/grid_generator_half_hyper_shell.output b/tests/grid/grid_generator_half_hyper_shell.output index 9cc562e669..fa9618d58c 100644 --- a/tests/grid/grid_generator_half_hyper_shell.output +++ b/tests/grid/grid_generator_half_hyper_shell.output @@ -12,15 +12,15 @@ DEAL::half_hyper_shell<2>:: cell 4, face 1, id 2 DEAL::half_hyper_shell<2>:: cell 4, face 2, id 1 DEAL::half_hyper_shell<2>:: cell 4, face 3, id 0 DEAL::half_hyper_shell<3>:: cell 0, face 2, id 2 -DEAL::half_hyper_shell<3>:: cell 0, face 4, id 0 +DEAL::half_hyper_shell<3>:: cell 0, face 4, id 1 DEAL::half_hyper_shell<3>:: cell 0, face 5, id 0 -DEAL::half_hyper_shell<3>:: cell 1, face 0, id 0 +DEAL::half_hyper_shell<3>:: cell 1, face 0, id 1 DEAL::half_hyper_shell<3>:: cell 1, face 1, id 0 DEAL::half_hyper_shell<3>:: cell 1, face 2, id 2 DEAL::half_hyper_shell<3>:: cell 2, face 2, id 2 -DEAL::half_hyper_shell<3>:: cell 2, face 4, id 0 +DEAL::half_hyper_shell<3>:: cell 2, face 4, id 1 DEAL::half_hyper_shell<3>:: cell 2, face 5, id 0 -DEAL::half_hyper_shell<3>:: cell 3, face 0, id 0 +DEAL::half_hyper_shell<3>:: cell 3, face 0, id 1 DEAL::half_hyper_shell<3>:: cell 3, face 1, id 0 DEAL::half_hyper_shell<3>:: cell 3, face 2, id 2 DEAL::half_hyper_shell<3>:: cell 4, face 0, id 1 diff --git a/tests/grid/grid_generator_half_hyper_shell_02.cc b/tests/grid/grid_generator_half_hyper_shell_02.cc new file mode 100644 index 0000000000..7d40812f9b --- /dev/null +++ b/tests/grid/grid_generator_half_hyper_shell_02.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 half_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 << "half_hyper_shell with origin: " << origin << std::endl; + Triangulation tr; + GridGenerator::half_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_half_hyper_shell_02.output b/tests/grid/grid_generator_half_hyper_shell_02.output new file mode 100644 index 0000000000..40a29dc9af --- /dev/null +++ b/tests/grid/grid_generator_half_hyper_shell_02.output @@ -0,0 +1,21 @@ + +DEAL:2d::half_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 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 1 0 2 1 0 +DEAL:2d::Ok. Found top boundary. +DEAL:2d::half_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 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 1 0 2 1 0 +DEAL:2d::Ok. Found top boundary. +DEAL:3d::half_hyper_shell with origin: 0.00000 0.00000 0.00000 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary faces: 2 1 0 1 0 2 2 1 0 1 0 2 1 0 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary lines: 1 0 2 2 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 2 2 1 0 1 0 2 2 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 2 2 1 0 1 1 1 1 0 0 0 0 +DEAL:3d::Ok. Found top boundary. +DEAL:3d::half_hyper_shell with origin: 1.00000 0.00000 0.00000 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary faces: 2 1 0 1 0 2 2 1 0 1 0 2 1 0 +DEAL:3d:: +DEAL:3d::Boundary ids of boundary lines: 1 0 2 2 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 2 2 1 0 1 0 2 2 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 2 2 1 0 1 1 1 1 0 0 0 0 +DEAL:3d::Ok. Found top boundary. -- 2.39.5