From 3a0dc2fd4a10d2386baf2a868287961b0525fbc4 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 4 Jan 2018 17:08:21 +0100 Subject: [PATCH] Implement CylindricalManifold::push_forward_gradient --- include/deal.II/grid/manifold_lib.h | 9 ++++++++ source/grid/manifold_lib.cc | 36 +++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index feac819df8..27f44282f1 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -381,6 +381,15 @@ public: virtual Point push_forward(const Point<3> &chart_point) const override; + /** + * Compute the derivatives of the mapping from cylindrical coordinates + * $(r, \phi, \lambda)$ to cartesian coordinates where $r$ denotes the + * distance from the axis, $\phi$ the angle between the given point and the + * computed normal direction and $\lambda$ the axial position. + */ + virtual DerivativeForm<1, 3, spacedim> + push_forward_gradient(const Point<3> &chart_point) const override; + /** * Compute new points on the CylindricalManifold. See the documentation of * the base class for a detailed description of what this function does. diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 72dabc7656..52e8361038 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -906,6 +906,42 @@ CylindricalManifold::push_forward(const Point<3> &chart_point) co +template +DerivativeForm<1, 3, spacedim> +CylindricalManifold::push_forward_gradient(const Point<3> &chart_point) const +{ + Tensor<2, 3> derivatives; + + // Rotate the orthogonal direction by the given angle. + // Formula from Section 5.2 in + // http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/ + // simplified assuming normal_direction and direction are orthogonal + // and unit vectors. + const double sine = std::sin(chart_point(1)); + const double cosine = std::cos(chart_point(1)); + const Tensor<1,3> dxn = cross_product_3d(direction, normal_direction); + const Tensor<1,3> intermediate = normal_direction*cosine+dxn*sine; + + // derivative w.r.t the radius + derivatives[0][0] = intermediate[0]; + derivatives[1][0] = intermediate[1]; + derivatives[2][0] = intermediate[2]; + + // derivatives w.r.t the angle + derivatives[0][1] = -normal_direction[0]*sine + dxn[0]*cosine; + derivatives[1][1] = -normal_direction[1]*sine + dxn[1]*cosine; + derivatives[2][1] = -normal_direction[2]*sine + dxn[2]*cosine; + + // derivatives w.r.t the direction of the axis + derivatives[0][2] = direction[0]; + derivatives[1][2] = direction[1]; + derivatives[2][2] = direction[2]; + + return derivatives; +} + + + // ============================================================ // FunctionManifold // ============================================================ -- 2.39.5