From: Bruno Blais Date: Tue, 3 Jun 2025 17:26:09 +0000 (-0400) Subject: WIP working on new cylinder grid generator X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=411270e84ee640af680fa142ae8241af7c25e96e;p=dealii.git WIP working on new cylinder grid generator --- diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index b41a125b76..ae410ef739 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -546,7 +546,16 @@ namespace GridGenerator const unsigned int n_shells = 2, const double skewness = 2.0, const bool colorize = false); - + template + void + custom_channel_with_cylinder(Triangulation &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 diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index ff80925162..22d635efd9 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -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(std::ceil((length_pre + length_post))); + const unsigned int height_repetitions = + 2*static_cast(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::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 *> 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 void