]> https://gitweb.dealii.org/ - dealii.git/commitdiff
WIP working on new cylinder grid generator
authorBruno Blais <blais.bruno@gmail.com>
Tue, 3 Jun 2025 17:26:09 +0000 (13:26 -0400)
committerBruno Blais <blais.bruno@gmail.com>
Tue, 3 Jun 2025 17:26:09 +0000 (13:26 -0400)
include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc

index b41a125b76c8d45cb51c415f218b98a03c889356..ae410ef739a04fa6c5b4818fff5f7120ec46bcc6 100644 (file)
@@ -546,7 +546,16 @@ namespace GridGenerator
                         const unsigned int  n_shells           = 2,
                         const double        skewness           = 2.0,
                         const bool          colorize           = false);
-
+  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  n_shells           = 2,
+    const double        skewness           = 2.0,
+    const bool          colorize           = false);
+                      
   /**
    * A general @p dim -dimensional cell (a segment if dim is 1, a quadrilateral
    * if @p dim is 2, or a hexahedron if @p dim is 3) immersed in a
index ff80925162c961ef45529ab775a1e5e6444b5370..22d635efd9ac4cb5ec14d83d22794b8071d9ae19 100644 (file)
@@ -3829,7 +3829,233 @@ namespace GridGenerator
           face->set_boundary_id(3);
   }
 
+  template <>
+  void
+  custom_channel_with_cylinder(Triangulation<2>  &tria,
+                        const double       half_height,
+                        const double       length_pre,
+                        const double       length_post,
+                        const unsigned int n_shells,
+                        const double       skewness,
+                        const bool         colorize)
+  {
+    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 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));
+
+
+    // We begin by setting up a grid that is 4 by 22 cells. While not
+    // squares, these have pretty good aspect ratios.
+    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));
+    // bulk_tria now looks like this:
+    //
+    //   +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
+    //   |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
+    //   +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
+    //   |  |XX|XX|  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
+    //   +--+--O--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
+    //   |  |XX|XX|  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
+    //   +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
+    //   |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
+    //   +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
+    //
+    // 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
+    //    (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)
+          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
+    // 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;
+        }
+
+    // Assign interior manifold ids to be the TFI id.
+    for (const auto &cell : cylinder_tria.active_cell_iterators())
+      {
+        cell->set_manifold_id(tfi_manifold_id);
+        for (const unsigned int face_n : GeometryInfo<2>::face_indices())
+          if (!cell->face(face_n)->at_boundary())
+            cell->face(face_n)->set_manifold_id(tfi_manifold_id);
+      }
+    if (0.0 < shell_region_width)
+      {
+        Assert(0 < n_shells,
+               ExcMessage("If the shell region has positive width then "
+                          "there must be at least one shell."));
+        Triangulation<2> shell_tria;
+        GridGenerator::concentric_hyper_shells(shell_tria,
+                                               Point<2>(),
+                                               radius,
+                                               radius + radius,
+                                               n_shells,
+                                               skewness,
+                                               8);
+
+        // Make the tolerance as large as possible since these cells can
+        // be quite close together
+        const double vertex_tolerance =
+          std::min(internal::minimal_vertex_distance(shell_tria),
+                   internal::minimal_vertex_distance(cylinder_tria)) *
+          0.5;
+
+        shell_tria.set_all_manifold_ids(polar_manifold_id);
+        Triangulation<2> temp;
+        GridGenerator::merge_triangulations(
+          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:
+    const double vertex_tolerance =
+      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())
+      if (cell->manifold_id() == polar_manifold_id)
+        cell->set_all_manifold_ids(polar_manifold_id);
+
+    // Ensure that all other manifold ids (including the interior faces
+    // opposite the cylinder) are set to the flat manifold id:
+    for (const auto &cell : tria.active_cell_iterators())
+      if (cell->manifold_id() != polar_manifold_id &&
+          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));
+    tria.set_manifold(polar_manifold_id, polar_manifold);
+
+    tria.set_manifold(tfi_manifold_id, FlatManifold<2>());
+    TransfiniteInterpolationManifold<2> inner_manifold;
+    inner_manifold.initialize(tria);
+    tria.set_manifold(tfi_manifold_id, inner_manifold);
+
+    if (colorize)
+      for (const auto &face : tria.active_face_iterators())
+        if (face->at_boundary())
+          {
+            const Point<2> center = face->center();
+            // left side
+            if (std::abs(center[0] - 0.0) < 1e-10)
+              face->set_boundary_id(0);
+            // right side
+            else if (std::abs(center[0] - 2.2) < 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
+            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 <int dim, int spacedim>
   void

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.