From dd165c62628f81ee11adb6788bda9094036856b0 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 6 Jul 2024 09:28:33 -0400 Subject: [PATCH] TorusManifold: rename R and r to centerline_radius and inner_radius. These names are self-describing. --- include/deal.II/grid/grid_generator.h | 12 ++-- include/deal.II/grid/manifold_lib.h | 12 ++-- source/grid/grid_generator.cc | 86 ++++++++++++++------------- source/grid/manifold_lib.cc | 52 +++++++++------- 4 files changed, 87 insertions(+), 75 deletions(-) diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 1d7bd02481..bb58a41af4 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -1648,10 +1648,12 @@ namespace GridGenerator * * @param tria The triangulation to be filled. * - * @param R The radius of the circle, which forms the middle line of the - * torus containing the loop of cells. Must be greater than @p r. + * @param centerline_radius The radius of the circle which forms the center + * line of the torus containing the loop of cells. Must be greater than @p + * inner_radius. * - * @param r The inner radius of the torus. + * @param inner_radius The distance between the inner edge of the torus and + * origin. * * @param n_cells_toroidal Optional argument to set the number of cell * layers in toroidal direction. The default is 6 cell layers. @@ -1668,8 +1670,8 @@ namespace GridGenerator template void torus(Triangulation &tria, - const double R, - const double r, + const double centerline_radius, + const double inner_radius, const unsigned int n_cells_toroidal = 6, const double phi = 2.0 * numbers::PI); diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 479bf52dfd..1cc44a64f1 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -786,11 +786,11 @@ public: static const int spacedim = 3; /** - * Constructor. Specify the radius of the centerline @p R and the radius - * of the torus itself (@p r). The variables have the same meaning as - * the parameters in GridGenerator::torus(). + * Constructor. Specify the radius of the centerline @p centerline_radius and + * the radius of the torus' inner circle (@p inner_radius). The variables have + * the same meaning as the parameters in GridGenerator::torus(). */ - TorusManifold(const double R, const double r); + TorusManifold(const double centerline_radius, const double inner_radius); /** * Make a clone of this Manifold object. @@ -817,7 +817,9 @@ public: push_forward_gradient(const Point<3> &chart_point) const override; private: - double r, R; + double centerline_radius; + + double inner_radius; }; diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 53bdfcfb50..3ef230e19a 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -1848,36 +1848,37 @@ namespace GridGenerator template <> void torus<2, 3>(Triangulation<2, 3> &tria, - const double R, - const double r, + const double centerline_radius, + const double inner_radius, const unsigned int, const double) { - Assert(R > r, - ExcMessage("Outer radius R must be greater than the inner " - "radius r.")); - Assert(r > 0.0, ExcMessage("The inner radius r must be positive.")); + Assert(centerline_radius > inner_radius, + ExcMessage("The centerline radius must be greater than the " + "inner radius.")); + Assert(inner_radius > 0.0, + ExcMessage("The inner radius must be positive.")); const unsigned int dim = 2; const unsigned int spacedim = 3; std::vector> vertices(16); - vertices[0] = Point(R - r, 0, 0); - vertices[1] = Point(R, -r, 0); - vertices[2] = Point(R + r, 0, 0); - vertices[3] = Point(R, r, 0); - vertices[4] = Point(0, 0, R - r); - vertices[5] = Point(0, -r, R); - vertices[6] = Point(0, 0, R + r); - vertices[7] = Point(0, r, R); - vertices[8] = Point(-(R - r), 0, 0); - vertices[9] = Point(-R, -r, 0); - vertices[10] = Point(-(R + r), 0, 0); - vertices[11] = Point(-R, r, 0); - vertices[12] = Point(0, 0, -(R - r)); - vertices[13] = Point(0, -r, -R); - vertices[14] = Point(0, 0, -(R + r)); - vertices[15] = Point(0, r, -R); + vertices[0] = Point(centerline_radius - inner_radius, 0, 0); + vertices[1] = Point(centerline_radius, -inner_radius, 0); + vertices[2] = Point(centerline_radius + inner_radius, 0, 0); + vertices[3] = Point(centerline_radius, inner_radius, 0); + vertices[4] = Point(0, 0, centerline_radius - inner_radius); + vertices[5] = Point(0, -inner_radius, centerline_radius); + vertices[6] = Point(0, 0, centerline_radius + inner_radius); + vertices[7] = Point(0, inner_radius, centerline_radius); + vertices[8] = Point(-(centerline_radius - inner_radius), 0, 0); + vertices[9] = Point(-centerline_radius, -inner_radius, 0); + vertices[10] = Point(-(centerline_radius + inner_radius), 0, 0); + vertices[11] = Point(-centerline_radius, inner_radius, 0); + vertices[12] = Point(0, 0, -(centerline_radius - inner_radius)); + vertices[13] = Point(0, -inner_radius, -centerline_radius); + vertices[14] = Point(0, 0, -(centerline_radius + inner_radius)); + vertices[15] = Point(0, inner_radius, -centerline_radius); std::vector> cells(16); // Right Hand Orientation @@ -1981,7 +1982,7 @@ namespace GridGenerator tria.create_triangulation(vertices, cells, SubCellData()); tria.set_all_manifold_ids(0); - tria.set_manifold(0, TorusManifold<2>(R, r)); + tria.set_manifold(0, TorusManifold<2>(centerline_radius, inner_radius)); } @@ -2000,15 +2001,16 @@ namespace GridGenerator template <> void torus<3, 3>(Triangulation<3, 3> &tria, - const double R, - const double r, + const double centerline_radius, + const double inner_radius, const unsigned int n_cells_toroidal, const double phi) { - Assert(R > r, - ExcMessage("Outer radius R must be greater than the inner " - "radius r.")); - Assert(r > 0.0, ExcMessage("The inner radius r must be positive.")); + Assert(centerline_radius > inner_radius, + ExcMessage("The centerline radius must be greater than the " + "inner radius.")); + Assert(inner_radius > 0.0, + ExcMessage("The inner radius must be positive.")); Assert(n_cells_toroidal > static_cast(phi / numbers::PI), ExcMessage("Number of cells in toroidal direction has " "to be at least 3 for a torus of polar extent 2*pi.")); @@ -2016,7 +2018,7 @@ namespace GridGenerator ExcMessage("Invalid angle phi specified.")); // the first 8 vertices are in the x-y-plane - const Point<3> p = Point<3>(R, 0.0, 0.0); + const Point<3> p = Point<3>(centerline_radius, 0.0, 0.0); const double a = 1. / (1 + std::sqrt(2.0)); // A value of 1 indicates "open" torus with angle < 2*pi, which // means that we need an additional layer of vertices @@ -2027,14 +2029,14 @@ namespace GridGenerator const unsigned int n_point_layers_toroidal = n_cells_toroidal + additional_layer; std::vector> vertices(8 * n_point_layers_toroidal); - 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), - vertices[3] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0) * a), - vertices[4] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0) * a), - vertices[5] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0) * a), - vertices[6] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0)), - vertices[7] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0)); + vertices[0] = p + Point<3>(-1, -1, 0) * (inner_radius / std::sqrt(2.0)), + vertices[1] = p + Point<3>(+1, -1, 0) * (inner_radius / std::sqrt(2.0)), + vertices[2] = p + Point<3>(-1, -1, 0) * (inner_radius / std::sqrt(2.0) * a), + vertices[3] = p + Point<3>(+1, -1, 0) * (inner_radius / std::sqrt(2.0) * a), + vertices[4] = p + Point<3>(-1, +1, 0) * (inner_radius / std::sqrt(2.0) * a), + vertices[5] = p + Point<3>(+1, +1, 0) * (inner_radius / std::sqrt(2.0) * a), + vertices[6] = p + Point<3>(-1, +1, 0) * (inner_radius / std::sqrt(2.0)), + vertices[7] = p + Point<3>(+1, +1, 0) * (inner_radius / std::sqrt(2.0)); // create remaining vertices by rotating around negative y-axis (the // direction is to ensure positive cell measures) @@ -2043,10 +2045,10 @@ namespace GridGenerator { for (unsigned int v = 0; v < 8; ++v) { - const double r_2d = vertices[v][0]; - vertices[8 * c + v][0] = r_2d * std::cos(phi_cell * c); + const double inner_radius_2d = vertices[v][0]; + vertices[8 * c + v][0] = inner_radius_2d * std::cos(phi_cell * c); vertices[8 * c + v][1] = vertices[v][1]; - vertices[8 * c + v][2] = r_2d * std::sin(phi_cell * c); + vertices[8 * c + v][2] = inner_radius_2d * std::sin(phi_cell * c); } } @@ -2101,7 +2103,7 @@ namespace GridGenerator } } - tria.set_manifold(1, TorusManifold<3>(R, r)); + tria.set_manifold(1, TorusManifold<3>(centerline_radius, inner_radius)); tria.set_manifold(2, CylindricalManifold<3>(Tensor<1, 3>({0., 1., 0.}), Point<3>())); diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index ddc375cb60..a7b4541939 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -1521,11 +1521,12 @@ TorusManifold::pull_back(const Point<3> &p) const double z = p[1]; double y = p[2]; double phi = std::atan2(y, x); - double theta = std::atan2(z, std::sqrt(x * x + y * y) - R); + double theta = std::atan2(z, std::sqrt(x * x + y * y) - centerline_radius); double w = - std::sqrt(Utilities::fixed_power<2>(y - std::sin(phi) * R) + - Utilities::fixed_power<2>(x - std::cos(phi) * R) + z * z) / - r; + std::sqrt(Utilities::fixed_power<2>(y - std::sin(phi) * centerline_radius) + + Utilities::fixed_power<2>(x - std::cos(phi) * centerline_radius) + + z * z) / + inner_radius; return {phi, theta, w}; } @@ -1539,23 +1540,26 @@ TorusManifold::push_forward(const Point<3> &chart_point) const double theta = chart_point[1]; double w = chart_point[2]; - return {std::cos(phi) * R + r * w * std::cos(theta) * std::cos(phi), - r * w * std::sin(theta), - std::sin(phi) * R + r * w * std::cos(theta) * std::sin(phi)}; + return {std::cos(phi) * centerline_radius + + inner_radius * w * std::cos(theta) * std::cos(phi), + inner_radius * w * std::sin(theta), + std::sin(phi) * centerline_radius + + inner_radius * w * std::cos(theta) * std::sin(phi)}; } template -TorusManifold::TorusManifold(const double R, const double r) +TorusManifold::TorusManifold(const double centerline_radius, + const double inner_radius) : ChartManifold(Point<3>(2 * numbers::PI, 2 * numbers::PI, 0.0)) - , r(r) - , R(R) + , centerline_radius(centerline_radius) + , inner_radius(inner_radius) { - Assert(R > r, - ExcMessage("Outer radius R must be greater than the inner " - "radius r.")); - Assert(r > 0.0, ExcMessage("inner radius must be positive.")); + Assert(centerline_radius > inner_radius, + ExcMessage("The centerline radius must be greater than the " + "inner radius.")); + Assert(inner_radius > 0.0, ExcMessage("The inner radius must be positive.")); } @@ -1564,7 +1568,7 @@ template std::unique_ptr> TorusManifold::clone() const { - return std::make_unique>(R, r); + return std::make_unique>(centerline_radius, inner_radius); } @@ -1579,17 +1583,19 @@ TorusManifold::push_forward_gradient(const Point<3> &chart_point) const double theta = chart_point[1]; double w = chart_point[2]; - DX[0][0] = -std::sin(phi) * R - r * w * std::cos(theta) * std::sin(phi); - DX[0][1] = -r * w * std::sin(theta) * std::cos(phi); - DX[0][2] = r * std::cos(theta) * std::cos(phi); + DX[0][0] = -std::sin(phi) * centerline_radius - + inner_radius * w * std::cos(theta) * std::sin(phi); + DX[0][1] = -inner_radius * w * std::sin(theta) * std::cos(phi); + DX[0][2] = inner_radius * std::cos(theta) * std::cos(phi); DX[1][0] = 0; - DX[1][1] = r * w * std::cos(theta); - DX[1][2] = r * std::sin(theta); + DX[1][1] = inner_radius * w * std::cos(theta); + DX[1][2] = inner_radius * std::sin(theta); - DX[2][0] = std::cos(phi) * R + r * w * std::cos(theta) * std::cos(phi); - DX[2][1] = -r * w * std::sin(theta) * std::sin(phi); - DX[2][2] = r * std::cos(theta) * std::sin(phi); + DX[2][0] = std::cos(phi) * centerline_radius + + inner_radius * w * std::cos(theta) * std::cos(phi); + DX[2][1] = -inner_radius * w * std::sin(theta) * std::sin(phi); + DX[2][2] = inner_radius * std::cos(theta) * std::sin(phi); return DX; } -- 2.39.5