From: Niklas Fehn Date: Wed, 16 Oct 2019 16:15:24 +0000 (+0200) Subject: extend GridGenerator::torus<3,3>() so that one can generate an open torus with angle... X-Git-Tag: v9.2.0-rc1~975^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3a36fead050155630c4d0f232b06b7a37be9b586;p=dealii.git extend GridGenerator::torus<3,3>() so that one can generate an open torus with angle 0 < phi <= 2*pi --- diff --git a/doc/news/changes/minor/20191016NiklasFehn b/doc/news/changes/minor/20191016NiklasFehn new file mode 100644 index 0000000000..ac34c61ef6 --- /dev/null +++ b/doc/news/changes/minor/20191016NiklasFehn @@ -0,0 +1,4 @@ +New: Extend function GridGenerator::torus<3,3>() so that one can generate +an open torus with an angle of 0 < phi <= 2*pi. +
+(Niklas Fehn, 2019/10/16) diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index e3033fffc8..616b3eaa48 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -1188,6 +1188,13 @@ namespace GridGenerator * @param n_cells_toroidal Optional argument to set the number of cell * layers in toroidal direction. The default is 6 cell layers. * + * @param phi Optional argument to generate an open torus with angle + * $0 < \varphi <= 2 \pi$. The default value is $2 \pi$, + * in which case a closed torus is generated. If the torus is open, + * the torus is cut at two planes perpendicular to the torus centerline. + * The center of these two planes are located at $(x_1, y_1, z_1) = (R, 0, 0)$ + * and $(x_2, y_2, z_2) = (R \cos(\varphi), 0, R \sin(\varphi))$. + * * @note Implemented for Triangulation<2,3> and Triangulation<3,3>. */ template @@ -1195,7 +1202,8 @@ namespace GridGenerator torus(Triangulation &tria, const double R, const double r, - const unsigned int n_cells_toroidal = 6); + const unsigned int n_cells_toroidal = 6, + const double phi = 2.0 * numbers::PI); /** * This function produces a square in the xy-plane with a cylindrical diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index a335abc88a..3523c89bd0 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -592,7 +592,8 @@ namespace GridGenerator void torus<2, 3>(Triangulation<2, 3> &tria, const double R, const double r, - const unsigned int) + const unsigned int, + const double) { Assert(R > r, ExcMessage("Outer radius R must be greater than the inner " @@ -734,7 +735,8 @@ namespace GridGenerator void torus<3, 3>(Triangulation<3, 3> &tria, const double R, const double r, - const unsigned int n_cells_toroidal) + const unsigned int n_cells_toroidal, + const double phi) { Assert(R > r, ExcMessage("Outer radius R must be greater than the inner " @@ -743,11 +745,16 @@ namespace GridGenerator Assert(n_cells_toroidal > 2, ExcMessage("Number of cells in toroidal direction has " "to be at least 3.")); + AssertThrow(phi > 0.0 && phi < 2.0 * numbers::PI + 1.0e-15, + ExcMessage("Invalid angle phi specified.")); // the first 8 vertices are in the x-y-plane - Point<3> const p = Point<3>(R, 0.0, 0.0); - double const a = 1. / (1 + std::sqrt(2.0)); - std::vector> vertices(8 * n_cells_toroidal); + Point<3> const p = Point<3>(R, 0.0, 0.0); + double const a = 1. / (1 + std::sqrt(2.0)); + unsigned int additional_layer = 0; // torus is closed (angle of 2*pi) + if (phi < 2.0 * numbers::PI - 1.0e-15) + additional_layer = 1; // indicates "open" torus with angle < 2*pi + std::vector> vertices(8 * (n_cells_toroidal + additional_layer)); vertices[0] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0)), vertices[1] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0)), vertices[2] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0) * a), @@ -759,8 +766,8 @@ namespace GridGenerator // create remaining vertices by rotating around negative y-axis (the // direction is to ensure positive cell measures) - double phi_cell = 2.0 * numbers::PI / n_cells_toroidal; - for (unsigned int c = 1; c < n_cells_toroidal; ++c) + double const phi_cell = phi / n_cells_toroidal; + for (unsigned int c = 1; c < n_cells_toroidal + additional_layer; ++c) { for (unsigned int v = 0; v < 8; ++v) { @@ -777,13 +784,15 @@ namespace GridGenerator { for (unsigned int j = 0; j < 2; ++j) { - unsigned int offset = (8 * (c + j)) % (8 * n_cells_toroidal); + unsigned int offset = + (8 * (c + j)) % (8 * (n_cells_toroidal + additional_layer)); + // cell 0 in x-y-plane cells[5 * c].vertices[0 + j * 4] = offset + 0; cells[5 * c].vertices[1 + j * 4] = offset + 1; cells[5 * c].vertices[2 + j * 4] = offset + 2; cells[5 * c].vertices[3 + j * 4] = offset + 3; - // cell 1 in x-y-plane + // cell 1 in x-y-plane (cell on torus centerline) cells[5 * c + 1].vertices[0 + j * 4] = offset + 2; cells[5 * c + 1].vertices[1 + j * 4] = offset + 3; cells[5 * c + 1].vertices[2 + j * 4] = offset + 4; @@ -805,8 +814,9 @@ namespace GridGenerator cells[5 * c + 4].vertices[3 + j * 4] = offset + 7; } - cells[5 * c].material_id = 0; - cells[5 * c + 1].material_id = 0; + cells[5 * c].material_id = 0; + // cell on torus centerline + cells[5 * c + 1].material_id = 1; cells[5 * c + 2].material_id = 0; cells[5 * c + 3].material_id = 0; cells[5 * c + 4].material_id = 0; @@ -819,14 +829,26 @@ namespace GridGenerator for (auto &cell : tria.cell_iterators()) { - bool cell_at_boundary = false; + // identify faces on torus surface and set manifold to 1 for (unsigned int f = 0; f < GeometryInfo<3>::faces_per_cell; ++f) - if (cell->at_boundary(f)) - cell_at_boundary = true; - if (cell_at_boundary == false) - cell->set_all_manifold_ids(2); + { + // faces 4 and 5 are those with normal vector aligned with torus + // centerline + if (cell->face(f)->at_boundary() && f != 4 && f != 5) + { + cell->face(f)->set_all_manifold_ids(1); + } + } + + // set manifold id to 2 for those cells that are on the torus centerline + if (cell->material_id() == 1) + { + cell->set_all_manifold_ids(2); + // reset to 0 + cell->set_material_id(0); + } } - tria.set_all_manifold_ids_on_boundary(1); + tria.set_manifold(1, TorusManifold<3>(R, r)); tria.set_manifold(2, CylindricalManifold<3>(Tensor<1, 3>({0., 1., 0.}), diff --git a/source/grid/grid_generator_from_name.cc b/source/grid/grid_generator_from_name.cc index 55256c090f..321824e268 100644 --- a/source/grid/grid_generator_from_name.cc +++ b/source/grid/grid_generator_from_name.cc @@ -275,9 +275,9 @@ namespace GridGenerator parse_and_create<3, 3, unsigned int, unsigned int, double, double>( moebius, arguments, tria); else if (name == "torus") - parse_and_create<3, 3, double, double, unsigned int>(torus, - arguments, - tria); + parse_and_create<3, 3, double, double, unsigned int, double>(torus, + arguments, + tria); else { return false; @@ -296,9 +296,9 @@ namespace GridGenerator Triangulation<2, 3> &tria) { if (name == "torus") - parse_and_create<2, 3, double, double, unsigned int>(torus, - arguments, - tria); + parse_and_create<2, 3, double, double, unsigned int, double>(torus, + arguments, + tria); else { return false;