From 5cbff7eb9a425e870517a2b5ce584bc2b412c1f4 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 3 Apr 2020 11:11:28 +0200 Subject: [PATCH] Add additional cell options for 3D hyper shell --- include/deal.II/grid/grid_generator.h | 23 ++- source/grid/grid_generator.cc | 204 ++++++++++++++------------ 2 files changed, 125 insertions(+), 102 deletions(-) diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 0684b44033..7210e43516 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -1050,12 +1050,23 @@ namespace GridGenerator * is zero (as is the default), then it is computed adaptively such that the * resulting elements have the least aspect ratio. * - * In 3d, only certain numbers are allowed, 6 (or the default 0) for a - * surface based on a hexahedron (i.e. 6 panels on the inner sphere extruded - * in radial direction to form 6 cells), 12 for the rhombic dodecahedron, - * and 96. This choice dates from an older version of deal.II before the - * Manifold classes were implemented: today all three choices are roughly - * equivalent (after performing global refinement, of course). + * In 3d, only certain numbers are allowed + * + * The versions with 24 and 48 cells are useful if the shell is thin and the + * radial lengths should be made more similar to the circumferential lengths. * * The grids with 12 and 96 cells are plotted below: * diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index c31e815115..0d02fa7aad 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -5361,111 +5361,123 @@ namespace GridGenerator {-1, +1, +1}, // {+1, +1, +1}}}; - // Start with the shell bounded by - // two nested cubes - if (n == 6) + switch (n) { - for (unsigned int i = 0; i < 8; ++i) - vertices.push_back(p + hexahedron[i] * irad); - for (unsigned int i = 0; i < 8; ++i) - vertices.push_back(p + hexahedron[i] * orad); - - const unsigned int n_cells = 6; - const int cell_vertices[n_cells][8] = { - {8, 9, 10, 11, 0, 1, 2, 3}, // bottom - {9, 11, 1, 3, 13, 15, 5, 7}, // right - {12, 13, 4, 5, 14, 15, 6, 7}, // top - {8, 0, 10, 2, 12, 4, 14, 6}, // left - {8, 9, 0, 1, 12, 13, 4, 5}, // front - {10, 2, 11, 3, 14, 6, 15, 7}}; // back - - cells.resize(n_cells, CellData<3>()); - - for (unsigned int i = 0; i < n_cells; ++i) + case 6: { - for (const unsigned int j : GeometryInfo<3>::vertex_indices()) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - } - - tria.create_triangulation(vertices, cells, SubCellData()); - } - // A more regular subdivision can - // be obtained by two nested - // rhombic dodecahedra + // Start with the shell bounded by two nested cubes + for (unsigned int i = 0; i < 8; ++i) + vertices.push_back(p + hexahedron[i] * irad); + for (unsigned int i = 0; i < 8; ++i) + vertices.push_back(p + hexahedron[i] * orad); + + const unsigned int n_cells = 6; + const int cell_vertices[n_cells][8] = { + {8, 9, 10, 11, 0, 1, 2, 3}, // bottom + {9, 11, 1, 3, 13, 15, 5, 7}, // right + {12, 13, 4, 5, 14, 15, 6, 7}, // top + {8, 0, 10, 2, 12, 4, 14, 6}, // left + {8, 9, 0, 1, 12, 13, 4, 5}, // front + {10, 2, 11, 3, 14, 6, 15, 7}}; // back + + cells.resize(n_cells, CellData<3>()); + + for (unsigned int i = 0; i < n_cells; ++i) + { + for (const unsigned int j : GeometryInfo<3>::vertex_indices()) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } - else if (n == 12) - { - // Octahedron inscribed in the cube [-1,1]^3 - static const std::array, 6> octahedron = {{{-1, 0, 0}, // - {1, 0, 0}, // - {0, -1, 0}, // - {0, 1, 0}, // - {0, 0, -1}, // - {0, 0, 1}}}; - - for (unsigned int i = 0; i < 8; ++i) - vertices.push_back(p + hexahedron[i] * irad); - for (unsigned int i = 0; i < 6; ++i) - vertices.push_back(p + octahedron[i] * inner_radius); - for (unsigned int i = 0; i < 8; ++i) - vertices.push_back(p + hexahedron[i] * orad); - for (unsigned int i = 0; i < 6; ++i) - vertices.push_back(p + octahedron[i] * outer_radius); - - const unsigned int n_cells = 12; - const unsigned int rhombi[n_cells][4] = {{10, 4, 0, 8}, - {4, 13, 8, 6}, - {10, 5, 4, 13}, - {1, 9, 10, 5}, - {9, 7, 5, 13}, - {7, 11, 13, 6}, - {9, 3, 7, 11}, - {1, 12, 9, 3}, - {12, 2, 3, 11}, - {2, 8, 11, 6}, - {12, 0, 2, 8}, - {1, 10, 12, 0}}; - - cells.resize(n_cells, CellData<3>()); - - for (unsigned int i = 0; i < n_cells; ++i) + tria.create_triangulation(vertices, cells, SubCellData()); + break; + } + case 12: { - for (unsigned int j = 0; j < 4; ++j) + // A more regular subdivision can be obtained by two nested rhombic + // dodecahedra + // + // Octahedron inscribed in the cube [-1,1]^3 + static const std::array, 6> octahedron = {{{-1, 0, 0}, // + {1, 0, 0}, // + {0, -1, 0}, // + {0, 1, 0}, // + {0, 0, -1}, // + {0, 0, 1}}}; + + for (unsigned int i = 0; i < 8; ++i) + vertices.push_back(p + hexahedron[i] * irad); + for (unsigned int i = 0; i < 6; ++i) + vertices.push_back(p + octahedron[i] * inner_radius); + for (unsigned int i = 0; i < 8; ++i) + vertices.push_back(p + hexahedron[i] * orad); + for (unsigned int i = 0; i < 6; ++i) + vertices.push_back(p + octahedron[i] * outer_radius); + + const unsigned int n_cells = 12; + const unsigned int rhombi[n_cells][4] = {{10, 4, 0, 8}, + {4, 13, 8, 6}, + {10, 5, 4, 13}, + {1, 9, 10, 5}, + {9, 7, 5, 13}, + {7, 11, 13, 6}, + {9, 3, 7, 11}, + {1, 12, 9, 3}, + {12, 2, 3, 11}, + {2, 8, 11, 6}, + {12, 0, 2, 8}, + {1, 10, 12, 0}}; + + cells.resize(n_cells, CellData<3>()); + + for (unsigned int i = 0; i < n_cells; ++i) { - cells[i].vertices[j] = rhombi[i][j]; - cells[i].vertices[j + 4] = rhombi[i][j] + 14; + for (unsigned int j = 0; j < 4; ++j) + { + cells[i].vertices[j] = rhombi[i][j]; + cells[i].vertices[j + 4] = rhombi[i][j] + 14; + } + cells[i].material_id = 0; } - cells[i].material_id = 0; - } - tria.create_triangulation(vertices, cells, SubCellData()); - } - else if (n == 96) - { - // create a triangulation based on the 12-cell version. This function - // was needed before SphericalManifold was written: it manually - // adjusted the interior vertices to lie along concentric - // spheres. Nowadays we can just refine globally: - Triangulation<3> tmp; - hyper_shell(tmp, p, inner_radius, outer_radius, 12); - tmp.refine_global(1); - - // now copy the resulting level 1 cells into the new triangulation, - cells.resize(tmp.n_active_cells(), CellData<3>()); - for (const auto &cell : tmp.active_cell_iterators()) + tria.create_triangulation(vertices, cells, SubCellData()); + break; + } + case 24: + case 48: { - const unsigned int cell_index = cell->active_cell_index(); - for (const unsigned int v : GeometryInfo<3>::vertex_indices()) - cells[cell_index].vertices[v] = cell->vertex_index(v); - cells[cell_index].material_id = 0; + // These two meshes are created by first creating a mesh of the + // 6-cell/12-cell version, refining globally once, and finally + // removing the outer half of the cells + Triangulation<3> tmp; + hyper_shell( + tmp, p, inner_radius, 2 * outer_radius - inner_radius, n / 4); + tmp.refine_global(1); + std::set::active_cell_iterator> cells_to_remove; + for (const auto &cell : tmp.active_cell_iterators()) + if (cell->center(true).norm_square() > + outer_radius * outer_radius) + cells_to_remove.insert(cell); + AssertDimension(cells_to_remove.size(), n); + create_triangulation_with_removed_cells(tmp, cells_to_remove, tria); + break; + } + case 96: + { + // create a triangulation based on the 12-cell version. This + // function was needed before SphericalManifold was written: it + // manually adjusted the interior vertices to lie along concentric + // spheres. Nowadays we can just refine globally: + Triangulation<3> tmp; + hyper_shell(tmp, p, inner_radius, outer_radius, 12); + tmp.refine_global(1); + flatten_triangulation(tmp, tria); + break; + } + default: + { + Assert(false, ExcMessage("Invalid number of coarse mesh cells.")); } - - tria.create_triangulation(tmp.get_vertices(), cells, SubCellData()); - } - else - { - Assert(false, ExcMessage("Invalid number of coarse mesh cells.")); } if (colorize) -- 2.39.5