]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Correctly use manifold in half hyper shell colorization 18574/head
authorRene Gassmoeller <rene.gassmoeller@mailbox.org>
Thu, 19 Jun 2025 23:15:51 +0000 (17:15 -0600)
committerRene Gassmoeller <rene.gassmoeller@mailbox.org>
Fri, 20 Jun 2025 17:19:06 +0000 (11:19 -0600)
doc/news/changes/minor/20250619Gassmoeller [new file with mode: 0644]
source/grid/grid_generator.cc
tests/grid/grid_generator_half_hyper_shell.output
tests/grid/grid_generator_half_hyper_shell_02.cc [new file with mode: 0644]
tests/grid/grid_generator_half_hyper_shell_02.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20250619Gassmoeller b/doc/news/changes/minor/20250619Gassmoeller
new file mode 100644 (file)
index 0000000..ef3e029
--- /dev/null
@@ -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.
+<br>
+(Rene Gassmoeller, 2025/06/19)
index 3514905d2d1f008696edb5ab91c3dd3cc319cc47..c061b92e06f605355aa309f917b97847d2713153 100644 (file)
@@ -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));
   }
 
 
index 9cc562e6695fec6abf63a0b5f6c057b56272f10a..fa9618d58c03c883b76f43598f082adf2dee0de4 100644 (file)
@@ -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 (file)
index 0000000..7d40812
--- /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 half_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 << "half_hyper_shell with origin: " << origin << std::endl;
+    Triangulation<dim> 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<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_half_hyper_shell_02.output b/tests/grid/grid_generator_half_hyper_shell_02.output
new file mode 100644 (file)
index 0000000..40a29dc
--- /dev/null
@@ -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.

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.