From 92ff104980f391b664a0b53a079774095f4878ac Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 4 Jan 2018 17:01:41 +0100 Subject: [PATCH] Simplify CylindricalManifold::push_forward --- source/grid/manifold_lib.cc | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 3c79f1c550..72dabc7656 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -62,7 +62,7 @@ namespace internal // helper function to compute a vector orthogonal to a given one. Point<3> - compute_normal(const Tensor<1,3> &vector) + compute_normal(const Tensor<1,3> &vector, bool normalize=false) { Assert(vector.norm_square() != 0., ExcMessage("The direction parameter must not be zero!")); @@ -87,6 +87,8 @@ namespace internal normal[1]=-1.; normal[2]=(vector[0]+vector[1])/vector[2]; } + if (normalize) + normal /= normal.norm(); return normal; } } @@ -823,7 +825,7 @@ CylindricalManifold::CylindricalManifold(const Point &d const Point &point_on_axis_, const double tolerance) : ChartManifold(Tensor<1,3>({0,2.*numbers::PI,0})), - normal_direction(internal::compute_normal(direction_)), + normal_direction(internal::compute_normal(direction_, true)), direction (direction_/direction_.norm()), point_on_axis (point_on_axis_), tolerance(tolerance) @@ -891,25 +893,15 @@ CylindricalManifold::push_forward(const Point<3> &chart_point) co // 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. - const double sine = std::sin(chart_point(1)); - const double cosine = std::cos(chart_point(1)); - const double x = normal_direction[0]*cosine - -sine*(direction[2]*normal_direction[1] - -direction[1]*normal_direction[2]); - const double y = normal_direction[1]*cosine - -sine*(direction[0]*normal_direction[2] - -direction[2]*normal_direction[0]); - const double z = normal_direction[2]*cosine - -sine*(direction[1]*normal_direction[0] - -direction[0]*normal_direction[1]); - - // Rescale according to the given distance from the axis. - Point<3> intermediate (x,y,z); - intermediate *= chart_point(0)/std::sqrt(intermediate.square()); + // simplified assuming normal_direction and direction are orthogonal + // and unit vectors. + const double sine_r = std::sin(chart_point(1))*chart_point(0); + const double cosine_r = std::cos(chart_point(1))*chart_point(0); + const Tensor<1,3> dxn = cross_product_3d(direction, normal_direction); + const Tensor<1,3> intermediate = normal_direction*cosine_r+dxn*sine_r; // Finally, put everything together. - return intermediate+point_on_axis+direction*chart_point(2); + return point_on_axis+direction*chart_point(2)+intermediate; } -- 2.39.5