From c12828cc413e1e0ddde3158d28f80672a4744281 Mon Sep 17 00:00:00 2001 From: Bruno Blais Date: Sun, 25 Feb 2024 12:59:46 -0500 Subject: [PATCH] Fix colorization of BCs in cylinder shell --- source/grid/grid_generator.cc | 37 +++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) 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()) -- 2.39.5