]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove the hidden CylindricalManifold<dim, 3> code. 5852/head
authorDavid Wells <wellsd2@rpi.edu>
Sat, 3 Feb 2018 18:13:45 +0000 (13:13 -0500)
committerDavid Wells <wellsd2@rpi.edu>
Sat, 3 Feb 2018 18:13:45 +0000 (13:13 -0500)
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
source/grid/manifold_lib.cc

index 1392896a3c8c5d2e785a45a6833cad2c97ead387..3313575111692eb8a6d8d49726efa37fd7157b6f 100644 (file)
@@ -936,92 +936,6 @@ private:
   FlatManifold<dim> chart_manifold;
 };
 
-#ifndef DOXYGEN
-/**
- * Specialization for the only properly implemented spacedim parameter.
- */
-template <int dim>
-class CylindricalManifold<dim,3> : public ChartManifold<dim,3,3>
-{
-public:
-  /**
-   * Constructor. Using default values for the constructor arguments yields a
-   * cylinder along the x-axis (<tt>axis=0</tt>). Choose <tt>axis=1</tt> or
-   * <tt>axis=2</tt> 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<const Point<3>> &surrounding_points,
-                 const ArrayView<const double>   &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
index c17484212ba11cd79b01d08131796e711eb0d86e..255c06681966e33620da94dc0fdb4b9e05e82b3b 100644 (file)
@@ -67,6 +67,15 @@ namespace internal
   }
 
   // helper function to compute a vector orthogonal to a given one.
+  // does nothing unless spacedim == 3.
+  template <int spacedim>
+  Point<spacedim>
+  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<spacedim> &p1,
 
 
 
-template<int dim, int spacedim>
+template <int dim, int spacedim>
 Tensor<1, spacedim>
 SphericalManifold<dim, spacedim>::
 normal_vector (const typename Triangulation<dim, spacedim >::face_iterator &face,
@@ -905,24 +914,19 @@ get_new_point (const ArrayView<const Tensor<1,3>> &directions,
 // ============================================================
 // CylindricalManifold
 // ============================================================
-
 template <int dim, int spacedim>
 CylindricalManifold<dim, spacedim>::CylindricalManifold(const unsigned int axis,
-                                                        const double tolerance) :
+                                                        const double tolerance)
+  :
   CylindricalManifold<dim, spacedim>(Point<spacedim>::unit_vector(axis),
                                      Point<spacedim>(),
                                      tolerance)
-{}
-
-
-
-template <int dim>
-CylindricalManifold<dim, 3>::CylindricalManifold(const unsigned int axis,
-                                                 const double tolerance) :
-  CylindricalManifold<dim, 3>(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<dim, spacedim>::CylindricalManifold(const Point<spacedim> &d
                                                         const Point<spacedim> &point_on_axis_,
                                                         const double tolerance) :
   ChartManifold<dim,spacedim,3>(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<dim, spacedim>::CylindricalManifold(const Point<spacedim> &d
 
 
 
-template <int dim>
-CylindricalManifold<dim, 3>::CylindricalManifold(const Point<3> &direction_,
-                                                 const Point<3> &point_on_axis_,
-                                                 const double tolerance) :
-  ChartManifold<dim,3,3>(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 <int dim, int spacedim>
 Point<spacedim>
 CylindricalManifold<dim,spacedim>::
-get_new_point (const ArrayView<const Point<spacedim>> &/*surrounding_points*/,
-               const ArrayView<const double>          &/*weights*/) const
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+               const ArrayView<const double>          &weights) const
 {
   Assert (spacedim==3,
           ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
-  return Point<spacedim>();
-}
-
 
-
-template <int dim>
-Point<3>
-CylindricalManifold<dim,3>::
-get_new_point (const ArrayView<const Point<3>> &surrounding_points,
-               const ArrayView<const double>          &weights) const
-{
   // First check if the average in space lies on the axis.
-  Point<3> middle;
+  Point<spacedim> middle;
   double average_length = 0.;
   for (unsigned int i=0; i<surrounding_points.size(); ++i)
     {
@@ -988,34 +969,26 @@ get_new_point (const ArrayView<const Point<3>> &surrounding_points,
   const double lambda = middle*direction;
 
   if ((middle-direction*lambda).square() < tolerance*average_length)
-    return Point<3>()+direction*lambda;
+    return Point<spacedim>()+direction*lambda;
   else // If not, using the ChartManifold should yield valid results.
-    return ChartManifold<dim, 3, 3>::get_new_point(surrounding_points,
-                                                   weights);
+    return ChartManifold<dim, spacedim, 3>::get_new_point(surrounding_points,
+                                                          weights);
 }
 
 
 
 template <int dim, int spacedim>
 Point<3>
-CylindricalManifold<dim, spacedim>::pull_back(const Point<spacedim> &/*space_point*/) const
+CylindricalManifold<dim, spacedim>::pull_back(const Point<spacedim> &space_point) const
 {
   Assert (spacedim==3,
           ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
-  return Point<3>();
-}
 
-
-
-template <int dim>
-Point<3>
-CylindricalManifold<dim, 3>::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<spacedim> 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<dim, 3>::pull_back(const Point<3> &space_point) const
 
 template <int dim, int spacedim>
 Point<spacedim>
-CylindricalManifold<dim, spacedim>::push_forward(const Point<3> &/*chart_point*/) const
+CylindricalManifold<dim, spacedim>::push_forward(const Point<3> &chart_point) const
 {
   Assert (spacedim==3,
           ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
-  return Point<spacedim>();
-}
 
-
-
-template<int dim>
-Point<3>
-CylindricalManifold<dim, 3>::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<dim, 3>::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<dim, 3>::push_forward(const Point<3> &chart_point) const
 
 
 
-template<int dim, int spacedim>
+template <int dim, int spacedim>
 DerivativeForm<1, 3, spacedim>
-CylindricalManifold<dim, spacedim>::push_forward_gradient(const Point<3> &/*chart_point*/) const
+CylindricalManifold<dim, spacedim>::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<int dim>
-DerivativeForm<1, 3, 3>
-CylindricalManifold<dim, 3>::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<dim, 3>::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];

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.