From da85dbe30915285507047b8b53701d8994d77a90 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 9 Jul 2023 23:50:33 -0600 Subject: [PATCH] Make a function more robust. --- source/grid/grid_generator.cc | 57 ++++++++++++++++++----------------- 1 file changed, 30 insertions(+), 27 deletions(-) diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 2871cbbe39..b37029d32b 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -7488,46 +7488,49 @@ namespace GridGenerator Assert(L > 0, ExcMessage("Must give positive extension L")); Assert(Nz >= 1, ExcLowerRange(1, Nz)); + // Start with a cylinder shell with the correct inner and outer radius + // and as many layers as requested cylinder_shell(triangulation, L, inner_radius, outer_radius, 8, Nz); triangulation.set_all_manifold_ids(numbers::flat_manifold_id); + // Then loop over all vertices that are at the boundary (by looping + // over all cells, their faces, and if the face is at the boundary, + // their vertices. If we haven't touched that vertex yet, see if + // we need to move it from its cylinder mantle position to the + // outer boundary of the box. std::vector treated_vertices(triangulation.n_vertices(), false); for (const auto &cell : triangulation.active_cell_iterators()) { - for (auto f : GeometryInfo::face_indices()) + for (const auto f : cell->face_indices()) if (cell->face(f)->at_boundary()) { - for (const unsigned int v : cell->vertex_indices()) + for (const unsigned int v : cell->face(f)->vertex_indices()) { const unsigned int vv = cell->face(f)->vertex_index(v); if (treated_vertices[vv] == false) { treated_vertices[vv] = true; - for (unsigned int i = 0; i <= Nz; ++i) - { - double d = i * L / Nz; - switch (vv - i * 16) - { - case 1: - cell->face(f)->vertex(v) = - Point(outer_radius, outer_radius, d); - break; - case 3: - cell->face(f)->vertex(v) = - Point(-outer_radius, outer_radius, d); - break; - case 5: - cell->face(f)->vertex(v) = - Point(-outer_radius, -outer_radius, d); - break; - case 7: - cell->face(f)->vertex(v) = - Point(outer_radius, -outer_radius, d); - break; - default: - break; - } - } + + // The vertices we have to treat are the ones that + // have x=y or x=-y and are at the outer ring -- that is, + // they are on the diagonal in the x-y plane and radius + // equal to outer_radius. These need to be pulled out to + // the corner point of the square, i.e., their x and y + // coordinates need to be multiplied by sqrt(2), + // whereas the z coordinate remains unchanged: + const Point vertex_location = + cell->face(f)->vertex(v); + if ((std::fabs(std::fabs(vertex_location[0]) - + std::fabs(vertex_location[1])) < + 1e-12 * outer_radius) && + (std::fabs(vertex_location[0] * vertex_location[0] + + vertex_location[1] * vertex_location[1] - + outer_radius * outer_radius) < + 1e-12 * outer_radius)) + cell->face(f)->vertex(v) = + Point<3>(vertex_location[0] * std::sqrt(2.0), + vertex_location[1] * std::sqrt(2.0), + vertex_location[2]); } } } -- 2.39.5