From 13dfeb9c83c704eebc4b89b77165d0532d592e31 Mon Sep 17 00:00:00 2001 From: Stefano Dominici Date: Fri, 30 Nov 2018 17:56:26 +0100 Subject: [PATCH] Changelog entry added as minor/20181130StefanoDominici. Minor modifications: * Formatted code to add three spaces after each function. * Extended the documentation in test elliptical_manifold_01.cc. * Added @copydoc to function declarations that were not showing up in the class documentation. * Fixed grammar. --- .../changes/minor/20181130StefanoDominici | 5 ++ include/deal.II/grid/manifold_lib.h | 16 +++- source/grid/manifold_lib.cc | 87 ++++++++++++------- tests/manifold/elliptical_manifold_01.cc | 55 ++++++++++++ tests/manifold/elliptical_manifold_02.cc | 11 +++ 5 files changed, 138 insertions(+), 36 deletions(-) create mode 100644 doc/news/changes/minor/20181130StefanoDominici diff --git a/doc/news/changes/minor/20181130StefanoDominici b/doc/news/changes/minor/20181130StefanoDominici new file mode 100644 index 0000000000..1c4caab11b --- /dev/null +++ b/doc/news/changes/minor/20181130StefanoDominici @@ -0,0 +1,5 @@ +New: Class EllipticalManifold derived from ChartManifold, +valid only for dim=2 and spacedim=2. It maps points from a system of +cartesian coordinates to a system of elliptical coordinates and vice-versa. +
+(Stefano Dominici, 2018/11/30) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 4449d25445..dfcc581858 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -499,17 +499,25 @@ public: EllipticalManifold(const Point & center, const Tensor<1, spacedim> &major_axis_direction, const double c_parameter); - virtual ~EllipticalManifold() = default; virtual std::unique_ptr> clone() const override; + /** + * @copydoc ChartManifold::pull_back() + */ virtual Point pull_back(const Point &space_point) const override; + /** + * @copydoc ChartManifold::push_forward() + */ virtual Point push_forward(const Point &chart_point) const override; + /** + * @copydoc ChartManifold::push_forward_gradient() + */ virtual DerivativeForm<1, spacedim, spacedim> push_forward_gradient(const Point &chart_point) const override; @@ -530,10 +538,10 @@ protected: private: /** - * Return the periodicity of the coordinate variables of the manifold. + * @copydoc ChartManifold::get_periodicity() * - * For $dim=2$ and $spacedim=2$, first coordinate is non-periodic, while - * second coordinate has a periodicity of $2\pi$. + * For $\text{dim}=2$ and $\text{spacedim}=2$, the first coordinate is + * non-periodic, while the second coordinate has a periodicity of $2\pi$. */ static Tensor<1, spacedim> get_periodicity(); diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 957503a218..e809bcf27e 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -1148,6 +1148,7 @@ CylindricalManifold::push_forward_gradient( } + // ============================================================ // EllipticalManifold // ============================================================ @@ -1174,6 +1175,8 @@ EllipticalManifold::EllipticalManifold( direction /= direction_norm; } + + template std::unique_ptr> EllipticalManifold::clone() const @@ -1183,17 +1186,21 @@ EllipticalManifold::clone() const c_parameter); } + + template Tensor<1, spacedim> EllipticalManifold::get_periodicity() { Tensor<1, spacedim> periodicity; // The second elliptical coordinate is periodic, while the first is not. - // Enforce periodicity on last variable. + // Enforce periodicity on the last variable. periodicity[spacedim - 1] = 2.0 * numbers::PI; return periodicity; } + + template Point EllipticalManifold::push_forward(const Point &) const @@ -1201,6 +1208,9 @@ EllipticalManifold::push_forward(const Point &) const Assert(false, ExcNotImplemented()); return Point(); } + + + template <> Point<2> EllipticalManifold<2, 2>::push_forward(const Point<2> &chart_point) const @@ -1218,6 +1228,9 @@ EllipticalManifold<2, 2>::push_forward(const Point<2> &chart_point) const direction[1] * x + direction[0] * y); return p + center; } + + + template Point EllipticalManifold::pull_back(const Point &) const @@ -1225,6 +1238,9 @@ EllipticalManifold::pull_back(const Point &) const Assert(false, ExcNotImplemented()); return Point(); } + + + template <> Point<2> EllipticalManifold<2, 2>::pull_back(const Point<2> &space_point) const @@ -1235,11 +1251,12 @@ EllipticalManifold<2, 2>::pull_back(const Point<2> &space_point) const const double y0 = space_point[1] - center[1]; const double x = direction[0] * x0 + direction[1] * y0; const double y = -direction[1] * x0 + direction[0] * y0; - // From here we try to find solutions to equation - // x^2/(1-q)-y^2/q = c^2 for q, and q = -sinh^2(pt[0]), - // and to equation - // x^2/(1-p)-y^2/p = c^2 for p, and p = sin^2(pt[1]). - // Note that p and q are solutions to the same quadratic equation. + // From here we try to find solutions of two equations: + // x^2/(1-q)-y^2/q = c^2 for q, given q = -sinh^2(pt[0]), + // and + // x^2/(1-p)-y^2/p = c^2 for p, given p = sin^2(pt[1]). + // Note that, in the end, p and q are solutions to the same quadratic + // equation. const double x2 = x * x; const double y2 = y * y; // @@ -1252,50 +1269,56 @@ EllipticalManifold<2, 2>::pull_back(const Point<2> &space_point) const p = 0.0; if (p > 1.0) p = 1.0; - const double q = (-b - srdelta) / (2.0 * c2); - const bool x_pos = std::signbit(x) == 0; - const bool y_pos = std::signbit(y) == 0; - const double eta0 = std::asin(std::sqrt(p)); + const double q = (-b - srdelta) / (2.0 * c2); + const bool x_is_positive = !std::signbit(x); + const bool y_is_positive = !std::signbit(y); + const double eta0 = std::asin(std::sqrt(p)); // Given q = -sinh^2(pt[0]), // pt[0] is calculated by straight inversion since q <= 0 for any // pt[0]. Point<2> pt(std::log(std::sqrt(-q) + std::sqrt(1.0 - q)), eta0); // Unfolding pt[1] according to the quadrant. - if (x_pos && !y_pos) + if (x_is_positive && !y_is_positive) pt[1] = 2.0 * numbers::PI - eta0; - else if (!x_pos && y_pos) + else if (!x_is_positive && y_is_positive) pt[1] = numbers::PI - eta0; - else if (!x_pos && !y_pos) + else if (!x_is_positive && !y_is_positive) pt[1] = numbers::PI + eta0; return pt; } + + template DerivativeForm<1, spacedim, spacedim> EllipticalManifold::push_forward_gradient( - const Point &chart_point) const + const Point &) const { - DerivativeForm<1, spacedim, spacedim> dX; - const double ch = std::cosh(chart_point[0]); - const double sh = std::sinh(chart_point[0]); - const double cs = std::cos(chart_point[1]); - const double sn = std::sin(chart_point[1]); - switch (spacedim) - { - case 2: - { - dX[0][0] = c_parameter * sh * cs; - dX[0][1] = -c_parameter * ch * sn; - dX[1][0] = c_parameter * ch * sn; - dX[1][1] = c_parameter * sh * cs; - } - break; - default: - Assert(false, ExcNotImplemented()); - } + Assert(false, ExcNotImplemented()); + return DerivativeForm<1, spacedim, spacedim>(); +} + + + +template <> +DerivativeForm<1, 2, 2> +EllipticalManifold<2, 2>::push_forward_gradient( + const Point<2> &chart_point) const +{ + const double ch = std::cosh(chart_point[0]); + const double sh = std::sinh(chart_point[0]); + const double cs = std::cos(chart_point[1]); + const double sn = std::sin(chart_point[1]); + DerivativeForm<1, 2, 2> dX; + dX[0][0] = c_parameter * sh * cs; + dX[0][1] = -c_parameter * ch * sn; + dX[1][0] = c_parameter * ch * sn; + dX[1][1] = c_parameter * sh * cs; return dX; } + + // ============================================================ // FunctionManifold // ============================================================ diff --git a/tests/manifold/elliptical_manifold_01.cc b/tests/manifold/elliptical_manifold_01.cc index 5e58bde5e6..a9a6b0affd 100644 --- a/tests/manifold/elliptical_manifold_01.cc +++ b/tests/manifold/elliptical_manifold_01.cc @@ -41,6 +41,35 @@ // hyper_shell made of 8 elements. namespace { + /* The following system of equations + // / + // |x = ellPar*cosh(u)*cos(v) + // |y = ellPar*sinh(u)*sin(v) + // \ + // allows to transform a point (u,v) in the chart space to a + // point (x,y) in the cartesian space. + // + // Setting v=0, the system of equations becomes + // / + // |x = x_0 = ellPar*cosh(u) + // |y = 0 + // \ + // which represents the crossing point of an ellipsis, non-rotated and + // centered at the origin of the cartesian system, with its major axis. + // + // For the sake of completeness, by setting v=pi/2, one finds y_0 = + // ellPar*sinh(u) as the crossing point on the minor axis. + // + // This function takes in input a point pt={x_0,v} + // and transforms it into (x,y) in accordance with the first system of + // equations above. + // + // On a final note, the condition x_0 >= ellPar has to be satisfied. This is + // required because min(cosh(u)) = 1 is obtained for u = 0, and imposing v = 0 + // the minimum value that x can assume is x = ellPar. + // For the sake of completeness, it's worth to notice that there would be no + // restrictions on the input if we passed y_0, instead of x_0, as pt[0]. + */ Point<2> chart_to_cartesian(const Point<2> &pt, const double ellPar) { @@ -51,6 +80,9 @@ namespace c[1] = ellPar * sh * std::sin(pt[1]); return c; } + + + std::vector> generate_shell_points(const Point<2> ¢er, const double radius0, @@ -78,6 +110,21 @@ namespace } return points; } + + + + // Generate an hyper_shell over an EllipticalManifold having an abitrary + // center, and the major axis oriented in the direction of the x-axis. + // + // inner_radius and outer_radius parameters correspond to the + // distances from the center of the manifold of the + // crossing points of two concentric ellipsis with their major axis. + // + // Input parameters must respect the following constrains: + // c_parameter < inner_radius < outer_radius. + // + // To understand the constrain on the c_parameter refer to the + // documentation of function chart_to_cartesian() whithin this namespace. void build_simple_hyper_shell(Triangulation<2, 2> &grid, const Point<2> & center, const double inner_radius, @@ -107,7 +154,13 @@ namespace grid.set_manifold(0, EllipticalManifold<2, 2>(center, axis, c_param)); grid.set_all_manifold_ids(0); } + + + } // namespace + + + // Helper function // Generate a simple hyper_shell over an elliptical manifold centered at the // origin. Major axis is the x-axis. @@ -134,6 +187,8 @@ test(unsigned int ref = 1) gridout.write_msh(tria, deallog.get_file_stream()); } + + int main() { diff --git a/tests/manifold/elliptical_manifold_02.cc b/tests/manifold/elliptical_manifold_02.cc index 00d9e8b7e9..17b8966091 100644 --- a/tests/manifold/elliptical_manifold_02.cc +++ b/tests/manifold/elliptical_manifold_02.cc @@ -30,6 +30,8 @@ #define TOLERANCE 1e-12 + + // Helper functions template Point @@ -46,6 +48,9 @@ test_push(const Tensor<1, spacedim> &axis, deallog << chart_point << std::endl; return chart_point; } + + + template Point test_pull(const Tensor<1, spacedim> &axis, @@ -62,6 +67,8 @@ test_pull(const Tensor<1, spacedim> &axis, return space_point; } + + // Function that tests pull_back() and push_forward(). void local_test(const Tensor<1, 2> &axis, @@ -80,6 +87,8 @@ local_test(const Tensor<1, 2> &axis, } } + + void test() { @@ -103,6 +112,8 @@ test() } } + + int main() { -- 2.39.5