From 9fe262d9808a4e7146a03fa68d296254c5beaf26 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 3 Feb 2018 13:13:45 -0500 Subject: [PATCH] Remove the hidden CylindricalManifold code. It turns out that we have converted enough compile-time dimensionality checks to run-time dimensionality checks that we can use just one implementation for CylindricalManifold instead of two. --- include/deal.II/grid/manifold_lib.h | 86 --------------------- source/grid/manifold_lib.cc | 116 +++++++++------------------- 2 files changed, 36 insertions(+), 166 deletions(-) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 1392896a3c..3313575111 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -936,92 +936,6 @@ private: FlatManifold chart_manifold; }; -#ifndef DOXYGEN -/** - * Specialization for the only properly implemented spacedim parameter. - */ -template -class CylindricalManifold : public ChartManifold -{ -public: - /** - * Constructor. Using default values for the constructor arguments yields a - * cylinder along the x-axis (axis=0). Choose axis=1 or - * axis=2 for a tube along the y- or z-axis, respectively. The - * tolerance value is used to determine if a point is on the axis. - */ - CylindricalManifold (const unsigned int axis = 0, - const double tolerance = 1e-10); - - /** - * Constructor. If constructed with this constructor, the manifold described - * is a cylinder with an axis that points in direction #direction and goes - * through the given #point_on_axis. The direction may be arbitrarily - * scaled, and the given point may be any point on the axis. The tolerance - * value is used to determine if a point is on the axis. - */ - CylindricalManifold (const Point<3> &direction, - const Point<3> &point_on_axis, - const double tolerance = 1e-10); - - /** - * Compute the Cartesian coordinates for a point given in cylindrical - * coordinates. - */ - virtual Point<3> - pull_back(const Point<3> &space_point) const override; - - /** - * Compute the cylindrical coordinates $(r, \phi, \lambda)$ for the given - * point 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 Point<3> - 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, 3> - 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. - */ - virtual Point<3> - get_new_point (const ArrayView> &surrounding_points, - const ArrayView &weights) const override; - -protected: - /** - * A vector orthogonal to the normal direction. - */ - const Tensor<1,3> normal_direction; - - /** - * The direction vector of the axis. - */ - - const Tensor<1,3> direction; - /** - * An arbitrary point on the axis. - */ - const Point<3> point_on_axis; - -private: - /** - * Relative tolerance to measure zero distances. - */ - double tolerance; - -}; -#endif //DOXYGEN - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index c17484212b..255c066819 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -67,6 +67,15 @@ namespace internal } // helper function to compute a vector orthogonal to a given one. + // does nothing unless spacedim == 3. + template + Point + compute_normal(const Tensor<1, spacedim> &/*vector*/, + bool /*normalize*/=false) + { + return {}; + } + Point<3> compute_normal(const Tensor<1,3> &vector, bool normalize=false) { @@ -379,7 +388,7 @@ get_tangent_vector (const Point &p1, -template +template Tensor<1, spacedim> SphericalManifold:: normal_vector (const typename Triangulation::face_iterator &face, @@ -905,24 +914,19 @@ get_new_point (const ArrayView> &directions, // ============================================================ // CylindricalManifold // ============================================================ - template CylindricalManifold::CylindricalManifold(const unsigned int axis, - const double tolerance) : + const double tolerance) + : CylindricalManifold(Point::unit_vector(axis), Point(), tolerance) -{} - - - -template -CylindricalManifold::CylindricalManifold(const unsigned int axis, - const double tolerance) : - CylindricalManifold(Point<3>::unit_vector(axis), - Point<3>(), - tolerance) -{} +{ + // do not use static_assert to make dimension-independent programming + // easier. + Assert (spacedim==3, + ExcMessage("CylindricalManifold can only be used for spacedim==3!")); +} @@ -931,7 +935,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(Tensor<1, spacedim>()), + normal_direction(internal::compute_normal(direction_, true)), direction (direction_/direction_.norm()), point_on_axis (point_on_axis_), tolerance(tolerance) @@ -944,40 +948,17 @@ CylindricalManifold::CylindricalManifold(const Point &d -template -CylindricalManifold::CylindricalManifold(const Point<3> &direction_, - const Point<3> &point_on_axis_, - const double tolerance) : - ChartManifold(Tensor<1,3>({0,2.*numbers::PI,0})), - normal_direction(internal::compute_normal(direction_, true)), - direction (direction_/direction_.norm()), - point_on_axis (point_on_axis_), - tolerance(tolerance) -{} - - - template Point CylindricalManifold:: -get_new_point (const ArrayView> &/*surrounding_points*/, - const ArrayView &/*weights*/) const +get_new_point (const ArrayView> &surrounding_points, + const ArrayView &weights) const { Assert (spacedim==3, ExcMessage("CylindricalManifold can only be used for spacedim==3!")); - return Point(); -} - - -template -Point<3> -CylindricalManifold:: -get_new_point (const ArrayView> &surrounding_points, - const ArrayView &weights) const -{ // First check if the average in space lies on the axis. - Point<3> middle; + Point middle; double average_length = 0.; for (unsigned int i=0; i> &surrounding_points, const double lambda = middle*direction; if ((middle-direction*lambda).square() < tolerance*average_length) - return Point<3>()+direction*lambda; + return Point()+direction*lambda; else // If not, using the ChartManifold should yield valid results. - return ChartManifold::get_new_point(surrounding_points, - weights); + return ChartManifold::get_new_point(surrounding_points, + weights); } template Point<3> -CylindricalManifold::pull_back(const Point &/*space_point*/) const +CylindricalManifold::pull_back(const Point &space_point) const { Assert (spacedim==3, ExcMessage("CylindricalManifold can only be used for spacedim==3!")); - return Point<3>(); -} - - -template -Point<3> -CylindricalManifold::pull_back(const Point<3> &space_point) const -{ // First find the projection of the given point to the axis. - const Tensor<1,3> normalized_point = space_point-point_on_axis; + const Tensor<1,spacedim> normalized_point = space_point-point_on_axis; const double lambda = normalized_point*direction; - const Point<3> projection = point_on_axis + direction*lambda; - const Tensor<1,3> p_diff = space_point - projection; + const Point projection = point_on_axis + direction*lambda; + const Tensor<1,spacedim> p_diff = space_point - projection; // Then compute the angle between the projection direction and // another vector orthogonal to the direction vector. @@ -1031,19 +1004,11 @@ CylindricalManifold::pull_back(const Point<3> &space_point) const template Point -CylindricalManifold::push_forward(const Point<3> &/*chart_point*/) const +CylindricalManifold::push_forward(const Point<3> &chart_point) const { Assert (spacedim==3, ExcMessage("CylindricalManifold can only be used for spacedim==3!")); - return Point(); -} - - -template -Point<3> -CylindricalManifold::push_forward(const Point<3> &chart_point) const -{ // Rotate the orthogonal direction by the given angle. // Formula from Section 5.2 in // http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/ @@ -1051,8 +1016,8 @@ CylindricalManifold::push_forward(const Point<3> &chart_point) const // 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; + const Tensor<1, spacedim> dxn = cross_product_3d(direction, normal_direction); + const Tensor<1, spacedim> intermediate = normal_direction*cosine_r+dxn*sine_r; // Finally, put everything together. return point_on_axis+direction*chart_point(2)+intermediate; @@ -1060,22 +1025,13 @@ CylindricalManifold::push_forward(const Point<3> &chart_point) const -template +template DerivativeForm<1, 3, spacedim> -CylindricalManifold::push_forward_gradient(const Point<3> &/*chart_point*/) const +CylindricalManifold::push_forward_gradient(const Point<3> &chart_point) const { Assert (spacedim==3, ExcMessage("CylindricalManifold can only be used for spacedim==3!")); - return DerivativeForm<1,3,spacedim>(); -} - - - -template -DerivativeForm<1, 3, 3> -CylindricalManifold::push_forward_gradient(const Point<3> &chart_point) const -{ Tensor<2, 3> derivatives; // Rotate the orthogonal direction by the given angle. @@ -1085,8 +1041,8 @@ CylindricalManifold::push_forward_gradient(const Point<3> &chart_point) // 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; + const Tensor<1,spacedim> dxn = cross_product_3d(direction, normal_direction); + const Tensor<1,spacedim> intermediate = normal_direction*cosine+dxn*sine; // derivative w.r.t the radius derivatives[0][0] = intermediate[0]; -- 2.39.5