From aac755cb33b2dc01e9d91cce41d4d60e24bada06 Mon Sep 17 00:00:00 2001 From: Bruno Blais Date: Mon, 9 Jun 2025 20:47:39 -0400 Subject: [PATCH] WIP grid generator for the cylinder --- include/deal.II/grid/grid_generator.h | 85 +++++++++++- source/grid/grid_generator.cc | 186 ++++++++++++-------------- 2 files changed, 168 insertions(+), 103 deletions(-) diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index ae410ef739..3098805707 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -546,12 +546,91 @@ namespace GridGenerator const unsigned int n_shells = 2, const double skewness = 2.0, const bool colorize = false); + + /** + * Generate a grid consisting of a channel with a cylinder where the length + * and the height of the channel can be defined by the user. This generator is + * a generalized version of GridGenerator::channel_with_cylinder. It can be used + * for benchmarking Navier-Stokes solvers for various flows around a cylinder cases + * in 2D or 3D. The main limitation of this generator is that the diameter of the cylindre is + * fixed at one and that the dimensions of the channel must be an integer multiple of this diameter. + * Consequently, the length before the cylinder + * @length_pre, the length after the cylinder @length_post and the half height of the channel @half_height + * must be integer values. The geometry consists of a channel + * of size $[-L_{pre}, -H] \times [L_{post}, H] \times [0, W] $ (where the $z$ + * dimension is omitted in 2d) with a cylinder, parallel to the $z$ axis + * with diameter $1$, centered at $(0, 0, 0)$. The channel has three + * distinct regions: + *
    + *
  1. If @p n_shells is greater than zero, then there are that many shells + * centered around the cylinder,
  2. + *
  3. a blending region between the shells and the rest of the + * triangulation, and
  4. + *
  5. a bulk region consisting of Cartesian cells.
  6. + *
+ * Here is the grid (without additional global refinement) in 2d: + * + * @image html channel_with_cylinder_2d.png + * + * and in 3d: + * + * @image html channel_with_cylinder_3d.png + * + * The resulting Triangulation uses three manifolds: a PolarManifold (in 2d) + * or CylindricalManifold (in 3d) with manifold id $0$, a + * TransfiniteInterpolationManifold with manifold id $1$, and a FlatManifold + * everywhere else. For more information on this topic see + * @ref GlossManifoldIndicator "the glossary entry on manifold indicators". + * The + * cell faces on the cylinder and surrounding shells have manifold ids of + * $0$, while the cell volumes adjacent to the shells (or, if they do not + * exist, the cylinder) have a manifold id of $1$. Put another way: this + * grid uses TransfiniteInterpolationManifold to smoothly transition from + * the shells (generated with GridGenerator::concentric_hyper_shells) to the + * bulk region. All other cell volumes and faces have manifold id + * numbers::flat_manifold_id and use FlatManifold. All cells with id + * numbers::flat_manifold_id are rectangular prisms aligned with the + * coordinate axes. + * + * + * @param tria Triangulation to be created. Must be empty upon calling this + * function. + * + * @param half_height The half height of the channel (y- to 0 or 0 to y+). + * + * @param length_pre The length of the channel from the left side (x-) to the center of the cylinder (0) + * + * @param length_post The length of the channel from the cylinder (0) to the right side (x+) + * + * @param shell_region_radius Radius of the layer of shells around the cylinder. + * This value should be between larger than 0.5 (the radius of the cylinder) and smaller than 1 (the half-length of the box around the cylinder). + * + * @param n_shells Number of shells to use in the shell layer. + * + * @param skewness Parameter controlling how close the shells are + * to the cylinder: see the mathematical definition given in + * GridGenerator::concentric_hyper_shells. + * + * @param colorize If `true`, then assign different boundary ids to + * different parts of the boundary. For more + * information on boundary indicators see + * @ref GlossBoundaryIndicator "this glossary entry". + * The left boundary (at $x = -L_{pre}$) is assigned an id of $0$, the right + * boundary (at $x = L_{post}$) is assigned an id of $1$; the boundary of + * the obstacle in the middle (i.e., the circle in 2d or the cylinder + * walls in 3d) is assigned an id of $2$, the bottom wall (at $y=-H$) is assigned and id of $/$, the top wall (at $y=H$) + * is assigned an id of $4$. In 3D, the front wall ($z=0$) is assigned an id of $5$ and the back wall ($z=W$) is assigned + * an id of $6$. + */ template void custom_channel_with_cylinder(Triangulation &tria, - const double half_height, - const double length_pre, - const double length_post, + const unsigned int half_height, + const unsigned int length_pre, + const unsigned int length_post, + const double depth=1, + const unsigned int depth_division=1, + const double shell_region_radius = 0.75, const unsigned int n_shells = 2, const double skewness = 2.0, const bool colorize = false); diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 22d635efd9..7c494ada78 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -3832,39 +3832,41 @@ namespace GridGenerator template <> void custom_channel_with_cylinder(Triangulation<2> &tria, - const double half_height, - const double length_pre, - const double length_post, + const unsigned int half_height, + const unsigned int length_pre, + const unsigned int length_post, + [[maybe_unused]] const double depth, + [[maybe_unused]] unsigned int depth_division, + const double shell_region_radius, const unsigned int n_shells, const double skewness, const bool colorize) { - const types::manifold_id polar_manifold_id = 0; + const types::manifold_id polar_manifold_id = 0; const types::manifold_id tfi_manifold_id = 1; // The radius of the cylinder is 0.5, so the diameter is 1. const double radius = 0.5; - const double shell_region_width =radius; - - // We assume that the cylinder, centered at (0,0), with a diameter of 1 - // is the characteristic length of the channel. - // The number of repetitions is chosen to ensure that the cylinder - // occupies four cells. + const double box_radius = 1; + // We assume that the cylinder is centered at (0,0) and has a diameter of 1. + // We use the cylinder diameter as the characteristic length of the channel. + // The number of repetitions is chosen to ensure that the cylinder + // occupies four cells. - const unsigned int length_repetitions = - 2*static_cast(std::ceil((length_pre + length_post))); - const unsigned int height_repetitions = - 2*static_cast(std::ceil(2.0 * half_height)); + const unsigned int length_repetitions =length_pre+length_post; + const unsigned int height_repetitions =2*half_height; + const double x_length = -double(length_pre); - // We begin by setting up a grid that is 4 by 22 cells. While not - // squares, these have pretty good aspect ratios. + // We begin by setting up a grid that is length_repetition by height_repetitions cells. + // These cells are all square Triangulation<2> bulk_tria; GridGenerator::subdivided_hyper_rectangle(bulk_tria, - {length_repetitions, height_repetitions}, - Point<2>(length_pre, -half_height), - Point<2>(length_pre+length_post, half_height)); + {(length_repetitions), height_repetitions}, + Point<2>(-double(length_pre), -double(half_height)), + Point<2>(double(length_post), double(half_height))); + // bulk_tria now looks like this: // // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+ @@ -3878,42 +3880,26 @@ namespace GridGenerator // +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+ // // The next step is to remove the cells marked with XXs: we will place - // the grid around the cylinder there later. The next loop does two - // things: - // 1. Determines which cells need to be removed from the Triangulation + // the grid around the cylinder there later. The following loop determines which cells need to be removed from the Triangulation // (i.e., find the cells marked with XX in the picture). - // 2. Finds the location of the vertex marked with 'O' and uses that to - // calculate the shift vector for aligning cylinder_tria with - // tria_without_cylinder. std::set::active_cell_iterator> cells_to_remove; - Tensor<1, 2> cylinder_triangulation_offset; for (const auto &cell : bulk_tria.active_cell_iterators()) { - if ((cell->center() - Point<2>(0., 0.)).norm() < radius) + if ((cell->center() - Point<2>(0., 0.)).norm() < 1.1*box_radius) cells_to_remove.insert(cell); } + Triangulation<2> tria_without_cylinder; GridGenerator::create_triangulation_with_removed_cells( bulk_tria, cells_to_remove, tria_without_cylinder); - // set up the cylinder triangulation. Note that this function sets the + // Set up the cylinder triangulation. Note that this function sets the // manifold ids of the interior boundary cells to 0 // (polar_manifold_id). Triangulation<2> cylinder_tria; GridGenerator::hyper_cube_with_cylindrical_hole(cylinder_tria, - 0.5 + 0.5, - 0.41 / 4.0); - // The bulk cells are not quite squares, so we need to move the left - // and right sides of cylinder_tria inwards so that it fits in - // bulk_tria: - for (const auto &cell : cylinder_tria.active_cell_iterators()) - for (const unsigned int vertex_n : GeometryInfo<2>::vertex_indices()) - { - if (std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1e-10) - cell->vertex(vertex_n)[0] = -0.1; - else if (std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1e-10) - cell->vertex(vertex_n)[0] = 0.1; - } + shell_region_radius, + box_radius); // Assign interior manifold ids to be the TFI id. for (const auto &cell : cylinder_tria.active_cell_iterators()) @@ -3923,7 +3909,9 @@ namespace GridGenerator if (!cell->face(face_n)->at_boundary()) cell->face(face_n)->set_manifold_id(tfi_manifold_id); } - if (0.0 < shell_region_width) + + // The shell region should have a radius that is larger than the radius of the cylinder + if (radius < shell_region_radius ) { Assert(0 < n_shells, ExcMessage("If the shell region has positive width then " @@ -3932,7 +3920,7 @@ namespace GridGenerator GridGenerator::concentric_hyper_shells(shell_tria, Point<2>(), radius, - radius + radius, + shell_region_radius, n_shells, skewness, 8); @@ -3950,7 +3938,6 @@ namespace GridGenerator shell_tria, cylinder_tria, temp, vertex_tolerance, true); cylinder_tria = std::move(temp); } - GridTools::shift(cylinder_triangulation_offset, cylinder_tria); // Compute the tolerance again, since the shells may be very close to // each-other: @@ -3958,34 +3945,10 @@ namespace GridGenerator std::min(internal::minimal_vertex_distance(tria_without_cylinder), internal::minimal_vertex_distance(cylinder_tria)) / 10; + GridGenerator::merge_triangulations( tria_without_cylinder, cylinder_tria, tria, vertex_tolerance, true); - // Move the vertices in the middle of the faces of cylinder_tria slightly - // to give a better mesh quality. We have to balance the quality of these - // cells with the quality of the outer cells (initially rectangles). For - // constant radial distance, we would place them at the distance 0.1 * - // sqrt(2.) from the center. In case the shell region width is more than - // 0.1/6., we choose to place them at 0.1 * 4./3. from the center, which - // ensures that the shortest edge of the outer cells is 2./3. of the - // original length. If the shell region width is less, we make the edge - // length of the inner part and outer part (in the shorter x direction) - // the same. - { - const double shift = - std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.); - for (const auto &cell : tria.active_cell_iterators()) - for (const unsigned int v : GeometryInfo<2>::vertex_indices()) - if (cell->vertex(v).distance(Point<2>(0.1, 0.205)) < 1e-10) - cell->vertex(v) = Point<2>(0.2 - shift, 0.205); - else if (cell->vertex(v).distance(Point<2>(0.3, 0.205)) < 1e-10) - cell->vertex(v) = Point<2>(0.2 + shift, 0.205); - else if (cell->vertex(v).distance(Point<2>(0.2, 0.1025)) < 1e-10) - cell->vertex(v) = Point<2>(0.2, 0.2 - shift); - else if (cell->vertex(v).distance(Point<2>(0.2, 0.3075)) < 1e-10) - cell->vertex(v) = Point<2>(0.2, 0.2 + shift); - } - // Ensure that all manifold ids on a polar cell really are set to the // polar manifold id: for (const auto &cell : tria.active_cell_iterators()) @@ -3999,32 +3962,8 @@ namespace GridGenerator cell->manifold_id() != tfi_manifold_id) cell->set_all_manifold_ids(numbers::flat_manifold_id); - // We need to calculate the current center so that we can move it later: - // to start get a unique list of (points to) vertices on the cylinder - std::vector *> cylinder_pointers; - for (const auto &face : tria.active_face_iterators()) - if (face->manifold_id() == polar_manifold_id) - { - cylinder_pointers.push_back(&face->vertex(0)); - cylinder_pointers.push_back(&face->vertex(1)); - } - // de-duplicate - std::sort(cylinder_pointers.begin(), cylinder_pointers.end()); - cylinder_pointers.erase(std::unique(cylinder_pointers.begin(), - cylinder_pointers.end()), - cylinder_pointers.end()); - - // find the current center... - Point<2> center; - for (const Point<2> *const ptr : cylinder_pointers) - center += *ptr / double(cylinder_pointers.size()); - - // and recenter at (0.2, 0.2) - for (Point<2> *const ptr : cylinder_pointers) - *ptr += Point<2>(0.2, 0.2) - center; - // attach manifolds - PolarManifold<2> polar_manifold(Point<2>(0.2, 0.2)); + PolarManifold<2> polar_manifold(Point<2>(0., 0.)); tria.set_manifold(polar_manifold_id, polar_manifold); tria.set_manifold(tfi_manifold_id, FlatManifold<2>()); @@ -4038,25 +3977,72 @@ namespace GridGenerator { const Point<2> center = face->center(); // left side - if (std::abs(center[0] - 0.0) < 1e-10) + if (std::abs(center[0] - length_pre) < 1e-10) face->set_boundary_id(0); // right side - else if (std::abs(center[0] - 2.2) < 1e-10) + else if (std::abs(center[0] - length_post) < 1e-10) face->set_boundary_id(1); // cylinder boundary else if (face->manifold_id() == polar_manifold_id) face->set_boundary_id(2); - // sides of channel + // Top and bottom of channel else { - Assert(std::abs(center[1] - 0.00) < 1.0e-10 || - std::abs(center[1] - 0.41) < 1.0e-10, - ExcInternalError()); face->set_boundary_id(3); } } } + template <> + void + custom_channel_with_cylinder(Triangulation<3> &tria, + const unsigned int half_height, + const unsigned int length_pre, + const unsigned int length_post, + const double depth, + unsigned int depth_division, + const double shell_region_radius, + const unsigned int n_shells, + const double skewness, + const bool colorize) { + + Triangulation<2> tria_2; + custom_channel_with_cylinder(tria_2, + half_height, + length_pre, + length_post, + depth, + depth_division, + shell_region_radius, + n_shells, + skewness, + colorize); + + // extrude to 3d + extrude_triangulation(tria_2, depth_division, depth, tria, true); + + // set up the new 3d manifolds + const types::manifold_id cylindrical_manifold_id = 0; + const types::manifold_id tfi_manifold_id = 1; + const PolarManifold<2> *const m_ptr = + dynamic_cast *>( + &tria_2.get_manifold(cylindrical_manifold_id)); + Assert(m_ptr != nullptr, ExcInternalError()); + const Point<3> axial_point(m_ptr->center[0], m_ptr->center[1], 0.0); + const Tensor<1, 3> direction{{0.0, 0.0, 1.0}}; + + tria.set_manifold(cylindrical_manifold_id, FlatManifold<3>()); + tria.set_manifold(tfi_manifold_id, FlatManifold<3>()); + const CylindricalManifold<3> cylindrical_manifold(direction, axial_point); + TransfiniteInterpolationManifold<3> inner_manifold; + inner_manifold.initialize(tria); + tria.set_manifold(cylindrical_manifold_id, cylindrical_manifold); + tria.set_manifold(tfi_manifold_id, inner_manifold); + + // From extrude_triangulation: since the maximum boundary id of tria_2 was + // 4, the front boundary id is 4 and the back is 5. They remain unchanged. + } + template void hyper_cross(Triangulation &tria, -- 2.39.5