From: Daniel Arndt Date: Tue, 2 Jan 2018 23:38:34 +0000 (+0100) Subject: Dummy implementation for CylindricalManifold for spacedim!=3 X-Git-Tag: v9.0.0-rc1~590^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=abf0f4c241f3130b9d77a678a62ae29eb1d60986;p=dealii.git Dummy implementation for CylindricalManifold for spacedim!=3 --- diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 27f44282f1..7686a0d7d0 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -920,6 +920,92 @@ 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 52e8361038..6840286dcd 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -813,19 +813,29 @@ get_new_point (const ArrayView> &directions, template CylindricalManifold::CylindricalManifold(const unsigned int axis, const double tolerance) : - CylindricalManifold(Point<3>::unit_vector(axis), - Point<3>(), + 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) +{} + + + template CylindricalManifold::CylindricalManifold(const Point &direction_, const Point &point_on_axis_, const double tolerance) : - ChartManifold(Tensor<1,3>({0,2.*numbers::PI,0})), - normal_direction(internal::compute_normal(direction_, true)), + ChartManifold(Tensor<1,3>({0,2.*numbers::PI,0})), + normal_direction(Tensor<1, spacedim>()), direction (direction_/direction_.norm()), point_on_axis (point_on_axis_), tolerance(tolerance) @@ -838,14 +848,40 @@ 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 +{ + 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 middle; + Point<3> 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()+direction*lambda; + return Point<3>()+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); } @@ -867,6 +903,17 @@ get_new_point (const ArrayView> &surrounding_points, template Point<3> 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; @@ -889,6 +936,17 @@ CylindricalManifold::pull_back(const Point &space_point template Point 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 @@ -909,6 +967,18 @@ CylindricalManifold::push_forward(const Point<3> &chart_point) co template DerivativeForm<1, 3, spacedim> 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; diff --git a/source/grid/manifold_lib.inst.in b/source/grid/manifold_lib.inst.in index d9d00a7c7b..7e181074c7 100644 --- a/source/grid/manifold_lib.inst.in +++ b/source/grid/manifold_lib.inst.in @@ -13,17 +13,13 @@ // // --------------------------------------------------------------------- -for (deal_II_dimension : DIMENSIONS) -{ - template class CylindricalManifold; -} - for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension template class PolarManifold; template class SphericalManifold; template class TransfiniteInterpolationManifold; + template class CylindricalManifold; #endif #if deal_II_dimension == deal_II_space_dimension template class TorusManifold;