]> https://gitweb.dealii.org/ - dealii.git/commitdiff
WIP grid generator for the cylinder
authorBruno Blais <blais.bruno@gmail.com>
Tue, 10 Jun 2025 00:47:39 +0000 (20:47 -0400)
committerBruno Blais <blais.bruno@gmail.com>
Tue, 10 Jun 2025 00:47:39 +0000 (20:47 -0400)
include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc

index ae410ef739a04fa6c5b4818fff5f7120ec46bcc6..30988057078e5a7323b477fc56424047808e2172 100644 (file)
@@ -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:
+   * <ol>
+   *   <li>If @p n_shells is greater than zero, then there are that many shells
+   *   centered around the cylinder,</li>
+   *   <li>a blending region between the shells and the rest of the
+   *   triangulation, and</li>
+   *   <li>a bulk region consisting of Cartesian cells.</li>
+   * </ol>
+   * 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 <int dim>
   void
   custom_channel_with_cylinder(Triangulation<dim> &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);
index 22d635efd9ac4cb5ec14d83d22794b8071d9ae19..7c494ada788ccbbd491f5db065e16ba304ae4ea6 100644 (file)
@@ -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<unsigned int>(std::ceil((length_pre + length_post)));
-    const unsigned int height_repetitions =
-      2*static_cast<unsigned int>(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<Triangulation<2>::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<Point<2> *> 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<const PolarManifold<2> *>(
+        &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 <int dim, int spacedim>
   void
   hyper_cross(Triangulation<dim, spacedim>    &tria,

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.