]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix quarter hyper shell colorization for thin shells.
authorRene Gassmoeller <rene.gassmoeller@mailbox.org>
Wed, 18 Jun 2025 17:57:14 +0000 (11:57 -0600)
committerRene Gassmoeller <rene.gassmoeller@mailbox.org>
Fri, 20 Jun 2025 00:04:44 +0000 (18:04 -0600)
doc/news/changes/minor/20250618Gassmoeller [new file with mode: 0644]
source/grid/grid_generator.cc
tests/grid/grid_generator_hyper_shell.cc [new file with mode: 0644]
tests/grid/grid_generator_hyper_shell.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20250618Gassmoeller b/doc/news/changes/minor/20250618Gassmoeller
new file mode 100644 (file)
index 0000000..8539a66
--- /dev/null
@@ -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.
+<br>
+(Rene Gassmoeller, 2025/06/18)
index 3514905d2d1f008696edb5ab91c3dd3cc319cc47..2d244b404b7c43679621b76d9f3183a89f0b6ebc 100644 (file)
@@ -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 (file)
index 0000000..df30dfe
--- /dev/null
@@ -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 <deal.II/base/tensor.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test(const Point<dim> origin = Point<dim>())
+{
+  SphericalManifold<dim> manifold(origin);
+  {
+    deallog << "quarter_hyper_shell with origin: " << origin << std::endl;
+    Triangulation<dim> 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<dim>::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<dim>::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 (file)
index 0000000..a8d8659
--- /dev/null
@@ -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.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.