Fix colorization of BCs in cylinder shell

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index e417ce880f76f2450d2e62dc0a21d12b6805ac0d..514c8a90720c760e08bdfcc4614eae28e6af99db 100644 (file)
--- 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())