/**
- * Manifold description for a polar space coordinate system. Given any
- * two points in space, tranform them in polar coordinate centered at
- * center, compute intermediate points in polar coordinate systems,
- * and returns its space representation.
+ * Manifold description for a polar space coordinate system. Given a
+ * set of points in space, this class computes their weighted average
+ * and the weighted average of their distance from the center. The
+ * average of the points is then shifted in the radial direction, to
+ * match the weighted average of the distances.
+ *
+ * You can use this Manifold object to describe any sphere, circle,
+ * hypersphere or hyperdisc in two or three dimensions, both as a
+ * co-dimension one manifold descriptor or as co-dimension zero
+ * manifold descriptor.
+ *
+ * @ingroup manifold
+ *
+ * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007, Luca Heltai,
+ * 2013
*/
template <int spacedim>
class PolarManifold : public FlatManifold<spacedim>
{
public:
+ /**
+ * The Constructor takes the
+ * center of the polar
+ * coordinates system.
+ */
PolarManifold(const Point<spacedim> center = Point<spacedim>());
+ /**
+ * Given a set of points in
+ * space, this function computes
+ * their weighted average and the
+ * weighted average of their
+ * distance from the center. The
+ * average of the points is then
+ * shifted in the radial
+ * direction, to match the
+ * weighted average of the
+ * distances.
+ */
virtual Point<spacedim>
get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &weights) const;
private:
- Point<spacedim> center;
+ /**
+ * The center of the polar
+ * coordinate system.
+ */
+ const Point<spacedim> center;
+};
+
+/**
+ * Manifold description for a cylindrical space coordinate system. In
+ * three dimensions, points are projected on a cylindrical tube along
+ * the <tt>x-</tt>, <tt>y-</tt> or <tt>z</tt>-axis (when using the
+ * first constructor of this class), or an arbitrarily oriented
+ * cylinder described by the direction of its axis and a point located
+ * on the axis. Similar to PolarManifold, new points are projected by
+ * taking the weighetd averages of the old points and adjusting the
+ * radius from the axis.
+ *
+ * This class was developed to be used in conjunction with the @p
+ * cylinder function of GridGenerator. Its use is discussed in detail
+ * in the results section of step-49.
+ *
+ * @ingroup manifold
+ *
+ * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007, Luca Heltai,
+ * 2013
+ */
+
+template <int spacedim>
+class CylinderManifold : public FlatManifold<spacedim>
+{
+ public:
+ /**
+ * Constructor. Using default
+ * values for the constructor
+ * arguments yields a cylindrical
+ * manifold tube along the x-axis
+ * (<tt>axis=0</tt>). Choose
+ * <tt>axis=1</tt> or
+ * <tt>axis=2</tt> for a
+ * cylindrical manifold along the
+ * y- or z-axis, respectively.
+ */
+ CylinderManifold(const int axis=0);
+
+ /**
+ * Constructor. If constructed
+ * with this constructor, the
+ * boundary 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.
+ */
+ CylinderManifold(const Point<spacedim> &direction,
+ const Point<spacedim> &point_on_axis);
+
+ /**
+ * Refer to the general
+ * documentation of this class
+ * and the documentation of the
+ * base class.
+ */
+ virtual Point<spacedim>
+ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+
+
+ private:
+ /**
+ * Direction of the cylindrical
+ * coordinate system.
+ */
+ const Point<spacedim> direction;
+
+ /**
+ * One point on the axis of the
+ * cylindrical coordinate systems
+ */
+ const Point<spacedim> point_on_axis;
+
+ /**
+ * Returns the axis vector
+ * associated with the given axis
+ * number.
+ */
+ static Point<spacedim> get_axis_vector (const unsigned int axis);
};
-//
-// /**
-// * Boundary object for the hull of a cylinder. In three dimensions,
-// * points are projected on a circular tube along the <tt>x-</tt>,
-// * <tt>y-</tt> or <tt>z</tt>-axis (when using the first constructor of
-// * this class), or an arbitrarily oriented cylinder described by the
-// * direction of its axis and a point located on the axis. The radius
-// * of the tube can be given independently. Similar to
-// * HyperBallBoundary, new points are projected by dividing the
-// * straight line between the old two points and adjusting the radius
-// * from the axis.
-// *
-// * This class was developed to be used in conjunction with the
-// * @p cylinder function of GridGenerator. It should be used for
-// * the hull of the cylinder only (boundary indicator 0). Its use is
-// * discussed in detail in the results section of step-49.
-// *
-// * This class is derived from StraightBoundary rather than from
-// * Boundary, which would seem natural, since this way we can use the
-// * StraightBoundary::in_between() function.
-// *
-// * @ingroup boundary
-// *
-// * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007
-// */
-// template <int dim, int spacedim = dim>
-// class CylinderBoundary : public StraightBoundary<dim,spacedim>
-// {
-// public:
-// /**
-// * Constructor. Using default values for the constructor arguments yields a
-// * circular tube 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.
-// */
-// CylinderBoundary (const double radius = 1.0,
-// const unsigned int axis = 0);
-//
-// /**
-// * Constructor. If constructed
-// * with this constructor, the
-// * boundary 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.
-// */
-// CylinderBoundary (const double radius,
-// const Point<spacedim> &direction,
-// const Point<spacedim> &point_on_axis);
-//
-// /**
-// * Refer to the general documentation of
-// * this class and the documentation of the
-// * base class.
-// */
-// virtual Point<spacedim>
-// get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
-//
-// /**
-// * Refer to the general documentation of
-// * this class and the documentation of the
-// * base class.
-// */
-// virtual Point<spacedim>
-// get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
-//
-// /**
-// * Refer to the general
-// * documentation of this class
-// * and the documentation of the
-// * base class.
-// *
-// * Calls
-// * @p get_intermediate_points_between_points.
-// */
-// virtual void
-// get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
-// std::vector<Point<spacedim> > &points) const;
-//
-// /**
-// * Refer to the general
-// * documentation of this class
-// * and the documentation of the
-// * base class.
-// *
-// * Only implemented for <tt>dim=3</tt>
-// * and for <tt>points.size()==1</tt>.
-// */
-// virtual void
-// get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
-// std::vector<Point<spacedim> > &points) const;
-//
-// /**
-// * Compute the normals to the
-// * boundary at the vertices of
-// * the given face.
-// *
-// * Refer to the general
-// * documentation of this class
-// * and the documentation of the
-// * base class.
-// */
-// virtual void
-// get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
-// typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
-//
-// /**
-// * Return the radius of the cylinder.
-// */
-// double get_radius () const;
-//
-// /**
-// * Exception. Thrown by the
-// * @p get_radius if the
-// * @p compute_radius_automatically,
-// * see below, flag is set true.
-// */
-// DeclException0 (ExcRadiusNotSet);
-//
-//
-// protected:
-// /**
-// * Radius of the cylinder.
-// */
-// const double radius;
-//
-// /**
-// * The direction vector of the axis.
-// */
-// const Point<spacedim> direction;
-//
-// /**
-// * An arbitrary point on the axis.
-// */
-// const Point<spacedim> point_on_axis;
-//
-// private:
-//
-// /**
-// * Called by
-// * @p get_intermediate_points_on_line
-// * and by
-// * @p get_intermediate_points_on_quad.
-// *
-// * Refer to the general
-// * documentation of
-// * @p get_intermediate_points_on_line
-// * in the documentation of the
-// * base class.
-// */
-// void get_intermediate_points_between_points (const Point<spacedim> &p0, const Point<spacedim> &p1,
-// std::vector<Point<spacedim> > &points) const;
-//
-// /**
-// * Given a number for the axis,
-// * return a vector that denotes
-// * this direction.
-// */
-// static Point<spacedim> get_axis_vector (const unsigned int axis);
-// };
//
//
//
return p;
}
+// ============================================================
+// CylinderManifold
+// ============================================================
+template <int spacedim>
+CylinderManifold<spacedim>::CylinderManifold(const int axis):
+ direction (get_axis_vector (axis)),
+ point_on_axis (Point<spacedim>())
+{}
+
+template <int spacedim>
+CylinderManifold<spacedim>::CylinderManifold(const Point<spacedim> &direction,
+ const Point<spacedim> &point_on_axis):
+ direction (direction),
+ point_on_axis (point_on_axis)
+{}
+
+
+
+template <int spacedim>
+Point<spacedim>
+CylinderManifold<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 spacedim>
+Point<spacedim>
+CylinderManifold<spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const
+{
+ Assert(surrounding_points.size() == weights.size(),
+ ExcDimensionMismatch(surrounding_points.size(), weights.size()));
+
+ // compute a proposed new point
+ Point<spacedim> middle = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
+
+ double radius = 0;
+ Point<spacedim> on_plane;
+
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ {
+ on_plane = surrounding_points[i]-point_on_axis;
+ on_plane = on_plane - (on_plane*direction) * direction;
+ radius += weights[i]*on_plane.norm();
+ }
+
+ // we then have to project this
+ // point out to the given radius
+ // from the axis. to this end, we
+ // have to take into account the
+ // offset point_on_axis and the
+ // direction of the axis
+ const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
+ ((middle-point_on_axis) * direction) * direction;
+
+ // scale it to the desired length
+ // and put everything back
+ // together, unless we have a point
+ // on the axis
+ if (vector_from_axis.norm() <= 1e-10 * middle.norm())
+ return middle;
+
+ else
+ return (vector_from_axis / vector_from_axis.norm() * radius +
+ ((middle-point_on_axis) * direction) * direction +
+ point_on_axis);
+}
+
// explicit instantiations
#include "manifold_lib.inst"
projects / dealii-svn.git / commitdiff