From dec8539cdd82da02e919e7196f3bd403b5f727d0 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 6 Jun 2017 15:14:14 +0200 Subject: [PATCH] Move computation of normal vector into anonymous namespace --- include/deal.II/grid/manifold_lib.h | 27 +++--- source/grid/manifold_lib.cc | 88 +++++++++---------- .../mappings/mapping_q_mixed_manifolds_01.cc | 24 ++++- .../mappings/mapping_q_mixed_manifolds_02.cc | 24 ++++- 4 files changed, 101 insertions(+), 62 deletions(-) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 53d59fed5d..60483e3af2 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -283,17 +283,17 @@ public: const double tolerance = 1e-10); /** - * Compute the cartesian coordinates for a point given in cylindrical + * Compute the Cartesian coordinates for a point given in cylindrical * coordinates. */ virtual Point<3> pull_back(const Point &space_point) const override; /** - * Compute the cylindrical coordinates \f$(r, \phi, \lambda)\f$ for the given - * point where \f$r\f$ denotes the distance from the axis, - * \f$\phi\f$ the angle between the given point and the computed normal - * direction and \f$\lambda\f$ the axial position. + * 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 push_forward(const Point<3> &chart_point) const override; @@ -306,13 +306,7 @@ public: get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const override; -private: - /** - * Compute a vector that is orthogonal to the given one. - * Used for initializing the member variable normal_direction. - */ - Point compute_normal(const Tensor<1,spacedim> &vector) const; - +protected: /** * A vector orthogonal to direcetion. */ @@ -328,10 +322,19 @@ private: */ const Point point_on_axis; +private: /** * Relative tolerance to measure zero distances. */ double tolerance; + + // explicitly check for sensible template arguments, but not on windows + // because MSVC creates bogus warnings during normal compilation +#ifndef DEAL_II_MSVC + static_assert (spacedim==3, + "CylindricalManifold can only be used for spacedim==3!"); +#endif + }; diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index f4d6e35663..3e8dba04a1 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -300,18 +300,49 @@ get_new_point (const std::vector > &vertices, // CylindricalManifold // ============================================================ +namespace +{ + // helper function to compute a vector orthogonal to a given one. + template + Point + compute_normal(const Tensor<1,spacedim> &vector) + { + AssertThrow(vector.norm() != 0., + ExcMessage("The direction parameter must not be zero!")); + Point<3> normal; + if (std::abs(vector[0]) >= std::abs(vector[1]) + && std::abs(vector[0]) >= std::abs(vector[2])) + { + normal[1]=-1.; + normal[2]=-1.; + normal[0]=(vector[1]+vector[2])/vector[0]; + } + else if (std::abs(vector[1]) >= std::abs(vector[0]) + && std::abs(vector[1]) >= std::abs(vector[2])) + { + normal[0]=-1.; + normal[2]=-1.; + normal[1]=(vector[0]+vector[2])/vector[1]; + } + else + { + normal[0]=-1.; + normal[1]=-1.; + normal[2]=(vector[0]+vector[1])/vector[2]; + } + return normal; + } +} + + + template CylindricalManifold::CylindricalManifold(const unsigned int axis, const double tolerance) : - ChartManifold(Tensor<1,3>({0,2.*numbers::PI,0})), - normal_direction(compute_normal(Point<3>::unit_vector(axis))), - direction (Point<3>::unit_vector(axis)), - point_on_axis (Point<3>()), - tolerance(tolerance) -{ - static_assert(spacedim==3, - "CylindricalManifold can only be used for spacedim==3!"); -} + CylindricalManifold(Point<3>::unit_vector(axis), + Point<3>(), + tolerance) +{} @@ -324,42 +355,7 @@ CylindricalManifold::CylindricalManifold(const Point &d direction (direction_/direction_.norm()), point_on_axis (point_on_axis_), tolerance(tolerance) -{ - static_assert(spacedim==3, - "CylindricalManifold can only be used for spacedim==3!"); -} - - - -template -Point -CylindricalManifold::compute_normal(const Tensor<1,spacedim> &vector) const -{ - AssertThrow(vector.norm() != 0., - ExcMessage("The direction parameter must not be zero!")); - Point<3> normal; - if (std::abs(vector[0]) >= std::abs(vector[1]) - && std::abs(vector[0]) >= std::abs(vector[2])) - { - normal[1]=-1.; - normal[2]=-1.; - normal[0]=(vector[1]+vector[2])/vector[0]; - } - else if (std::abs(vector[1]) >= std::abs(vector[0]) - && std::abs(vector[1]) >= std::abs(vector[2])) - { - normal[0]=-1.; - normal[2]=-1.; - normal[1]=(vector[0]+vector[2])/vector[1]; - } - else - { - normal[0]=-1.; - normal[1]=-1.; - normal[2]=(vector[0]+vector[1])/vector[2]; - } - return normal; -} +{} @@ -369,7 +365,7 @@ CylindricalManifold:: get_new_point (const std::vector > &surrounding_points, const std::vector &weights) const { - // First check if the the average in space lies on the axis. + // First check if the average in space lies on the axis. Point middle; double average_length = 0.; for (unsigned int i=0; i #include +template +struct ManifoldWrapper +{ + Manifold *operator()(const Point &direction, + const Point ¢er ) const; +}; +template <> +Manifold<2> * +ManifoldWrapper<2>::operator()(const Point<2> &/*direction*/, + const Point<2> ¢er ) const +{ + return new SphericalManifold<2>(center); +} + +template <> +Manifold<3> * +ManifoldWrapper<3>::operator()(const Point<3> &direction, + const Point<3> ¢er ) const +{ + return new CylindricalManifold<3>(direction, center); +} template void test() @@ -45,8 +66,7 @@ void test() direction[dim-1] = 1.; std::shared_ptr > cylinder_manifold - (dim == 2 ? static_cast*>(new SphericalManifold(Point())) : - static_cast*>(new CylindricalManifold(direction, Point()))); + (ManifoldWrapper()(direction, Point())); // create cube and shift it to position 0.7 Triangulation tria; diff --git a/tests/mappings/mapping_q_mixed_manifolds_02.cc b/tests/mappings/mapping_q_mixed_manifolds_02.cc index ec551b0559..9f525c88bd 100644 --- a/tests/mappings/mapping_q_mixed_manifolds_02.cc +++ b/tests/mappings/mapping_q_mixed_manifolds_02.cc @@ -150,7 +150,28 @@ void create_triangulation(Triangulation<3> &tria) } } +template +struct ManifoldWrapper +{ + Manifold *operator()(const Point &direction, + const Point ¢er ) const; +}; + +template <> +Manifold<2> * +ManifoldWrapper<2>::operator()(const Point<2> &/*direction*/, + const Point<2> ¢er ) const +{ + return new SphericalManifold<2>(center); +} +template <> +Manifold<3> * +ManifoldWrapper<3>::operator()(const Point<3> &direction, + const Point<3> ¢er ) const +{ + return new CylindricalManifold<3>(direction, center); +} template void test() @@ -162,8 +183,7 @@ void test() direction[dim-1] = 1.; std::shared_ptr > cylinder_manifold - (dim == 2 ? static_cast*>(new SphericalManifold(center)) : - static_cast*>(new CylindricalManifold(direction, center))); + (ManifoldWrapper()(direction, center)); Triangulation tria; create_triangulation(tria); tria.set_manifold(MANIFOLD_ID, *cylinder_manifold); -- 2.39.5