Made CylindricalManifold derived from Manifold instead of FlatManifold.

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index 8397511a05301d72c1d00688252f75f58ccaca58..fda040d710996feb36b7df766411d155f32c38cf 100644 (file)
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -124,26 +124,30 @@ private:
  * @author Luca Heltai, 2014
  */
 template <int dim, int spacedim = dim>
-class CylindricalManifold : public FlatManifold<dim,spacedim>
+class CylindricalManifold : public Manifold<dim,spacedim>
 {
 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.
+   * z-axis, respectively. The tolerance value is used to determine
+   * if a point is on the axis.
    */
-  CylindricalManifold (const unsigned int axis = 0);
+  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.
+   * any point on the axis. The tolerance value is used to determine
+   * if a point is on the axis. 
    */
   CylindricalManifold (const Point<spacedim> &direction,
-                      const Point<spacedim> &point_on_axis);
+                      const Point<spacedim> &point_on_axis, 
+                      const double tolerance = 1e-10);
 
  /**
    * Compute new points on the CylindricalManifold. See the documentation
@@ -166,10 +170,14 @@ protected:
 
 private:
   /**
-   * Given a number for the axis, return a vector that denotes this
-   * direction.
+   * Helper FlatManifold to compute temptative midpoints.
    */
-  static Point<spacedim> get_axis_vector (const unsigned int axis);
+  FlatManifold<dim,spacedim> flat_manifold;
+  
+  /**
+   * Relative tolerance to measure zero distances.
+   */
+  double tolerance;
 };
 
 DEAL_II_NAMESPACE_CLOSE

diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index efc765ff1768427ebd3778cb65fbafe268f9cc39..bd479a125568aa972f74841d4a121608c718fabc 100644 (file)
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -131,22 +131,12 @@ SphericalManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) c
 // CylindricalManifold
 // ============================================================
 
-
-template <int dim, int spacedim>
-Point<spacedim>
-CylindricalManifold<dim,spacedim>::get_axis_vector (const unsigned int axis)
-{
-  Assert (axis < spacedim, ExcIndexRange (axis, 0, spacedim));
-  Point<spacedim> axis_vector;
-  axis_vector[axis] = 1;
-  return axis_vector;
-}
-
-
 template <int dim, int spacedim>
-CylindricalManifold<dim,spacedim>::CylindricalManifold(const unsigned int axis) :
-  direction (get_axis_vector (axis)),
-  point_on_axis (Point<spacedim>())
+CylindricalManifold<dim,spacedim>::CylindricalManifold(const unsigned int axis, 
+                                                      const double tolerance) :
+  direction (Point<spacedim>::unit_vector(axis)),
+  point_on_axis (Point<spacedim>()),
+  tolerance(tolerance)
 {
   Assert(spacedim > 1, ExcImpossibleInDim(1));
 }
@@ -154,9 +144,11 @@ CylindricalManifold<dim,spacedim>::CylindricalManifold(const unsigned int axis)
 
 template <int dim, int spacedim>
 CylindricalManifold<dim,spacedim>::CylindricalManifold(const Point<spacedim> &direction,
-                                                      const Point<spacedim> &point_on_axis) :
+                                                      const Point<spacedim> &point_on_axis, 
+                                                      const double tolerance) :
   direction (direction),
-  point_on_axis (point_on_axis)
+  point_on_axis (point_on_axis),
+  tolerance(tolerance)
 {
   Assert(spacedim > 2, ExcImpossibleInDim(spacedim));
 }
@@ -173,7 +165,7 @@ get_new_point (const Quadrature<spacedim> &quad) const
   const std::vector<double> &weights = quad.get_weights();
 
   // compute a proposed new point  
-  Point<spacedim> middle = FlatManifold<dim,spacedim>::get_new_point(quad);
+  Point<spacedim> middle = flat_manifold.get_new_point(quad);
 
   double radius = 0;
   Point<spacedim> on_plane;
@@ -193,7 +185,7 @@ get_new_point (const Quadrature<spacedim> &quad) const
 
   // scale it to the desired length and put everything back together,
   // unless we have a point on the axis
-  if (vector_from_axis.norm() <= this->tolerance * middle.norm())
+  if (vector_from_axis.norm() <= tolerance * middle.norm())
     return middle;
 
   else