From: Bruno Blais <blais.bruno@gmail.com>
Date: Sun, 25 Feb 2024 17:59:46 +0000 (-0500)
Subject: Fix colorization of BCs in cylinder shell
X-Git-Tag: v9.6.0-rc1~519^2~5
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c12828cc413e1e0ddde3158d28f80672a4744281;p=dealii.git

Fix colorization of BCs in cylinder shell
---

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index e417ce880f..514c8a9072 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -6480,8 +6480,41 @@ namespace GridGenerator
     // Top (Z+) part of the cylinder has boundary id 3
 
     // Define tolerance to help detect boundary conditions
-    double eps_z =
-      std::min(1e-3 * (outer_radius - inner_radius), 1e-3 * length);
+    // First we define the tolerance along the z axis to identify 
+    // bottom and top cells.
+    double eps_z = 1e-6 * length;
+
+
+
+      // Gather the inner radius from the faces instead of the argument, this is
+      // more robust for some aspect ratios. First initialize the outer to 0 and
+      // the inner to a large value
+      double inner_radius = DBL_MAX;
+      double outer_radius = 0.;
+
+      // Loop over the cells once to acquire the min and max radius at the face
+      // centers Otherwise, for some cell ratio, the center of the faces can be
+      // at a radius which is significantly different from the one prescribed.
+      for (const auto &cell : triangulation.active_cell_iterators())
+        for (const unsigned int f : GeometryInfo<3>::face_indices())
+          {
+            if (!cell->face(f)->at_boundary())
+              continue;
+
+            const auto   face_center = cell->face(f)->center();
+            const double z           = face_center[2];
+
+            if ((std::fabs(z) > eps_z) &&
+                (std::fabs(z - length) > eps_z)) // Not a zmin or zmax boundary
+              {
+                const double radius =
+                  std::sqrt(face_center[0] * face_center[0] +
+                            face_center[1] * face_center[1]);
+                inner_radius = std::min(inner_radius, radius);
+                outer_radius = std::max(outer_radius, radius);
+              }
+          }
+    
     double mid_radial_distance = 0.5 * (outer_radius - inner_radius);
 
     for (const auto &cell : tria.active_cell_iterators())