// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
/**
* 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. The radius of the tube can be
- * set. Similar to HyperBallBoundary, new points are projected by
- * dividing the straight line between the old two points and adjusting
- * the radius in the @p yz-plane (xz-plane or xy-plane, respectively).
+ * <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
*
* @ingroup boundary
*
- * @author Guido Kanschat, 2001
+ * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007
*/
template <int dim>
class CylinderBoundary : public StraightBoundary<dim>
* along the y- or z-axis,
* respectively.
*/
- CylinderBoundary (const double radius = 1.0,
- const unsigned int axis = 0);
+ 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<dim> direction,
+ const Point<dim> point_on_axis);
/**
* Refer to the general documentation of
const double radius;
/**
- * Denotes the axis of the
- * circular
- * tube. <tt>axis=0</tt>,
- * <tt>axis=1</tt> and
- * <tt>axis=2</tt> denote the
- * <tt>x-axis</tt>,
- * <tt>y-axis</tt> and
- * <tt>z-axis</tt>, respectively.
+ * The direction vector of the axis.
+ */
+ const Point<dim> direction;
+
+ /**
+ * An arbitrary point on the axis.
*/
- const unsigned int axis;
+ const Point<dim> point_on_axis;
private:
* base class.
*/
void get_intermediate_points_between_points (const Point<dim> &p0, const Point<dim> &p1,
- std::vector<Point<dim> > &points) const;
+ std::vector<Point<dim> > &points) const;
+
+ /**
+ * Given a number for the axis,
+ * return a vector that denotes
+ * this direction.
+ */
+ static Point<dim> get_axis_vector (const unsigned int axis);
};
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
DEAL_II_NAMESPACE_OPEN
+
+
+
template <int dim>
CylinderBoundary<dim>::CylinderBoundary (const double radius,
- const unsigned int axis) :
+ const unsigned int axis)
+ :
radius(radius),
- axis(axis)
+ direction (get_axis_vector (axis)),
+ point_on_axis (Point<dim>())
{}
+template <int dim>
+CylinderBoundary<dim>::CylinderBoundary (const double radius,
+ const Point<dim> direction,
+ const Point<dim> point_on_axis)
+ :
+ radius(radius),
+ direction (direction / direction.norm()),
+ point_on_axis (point_on_axis)
+{}
+
template <int dim>
Point<dim>
-CylinderBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+CylinderBoundary<dim>::get_axis_vector (const unsigned int axis)
{
- Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
- // project to boundary
- if (dim>=3
- && line->vertex(0).square()-line->vertex(0)(axis)*line->vertex(0)(axis) >= radius*radius-1.e-10
- && line->vertex(1).square()-line->vertex(1)(axis)*line->vertex(1)(axis) >= radius*radius-1.e-10)
- {
- const double f = radius / std::sqrt(middle.square()-middle(axis)*middle(axis));
- for (unsigned int i=0; i<dim; ++i)
- if (i!=axis)
- middle(i) *= f;
- }
- return middle;
+ Assert (axis < dim, ExcIndexRange (axis, 0, dim));
+
+ Point<dim> axis_vector;
+ axis_vector[axis] = 1;
+ return axis_vector;
+}
+
+
+
+template <int dim>
+Point<dim>
+CylinderBoundary<dim>::
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+ // compute a proposed new point
+ const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+
+ // 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<dim> 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);
}
CylinderBoundary<3>::
get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
{
- Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+ const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
- // project to boundary
- if (quad->vertex(0).square()-quad->vertex(0)(axis)*quad->vertex(0)(axis) >= radius*radius-1.e-10
- && quad->vertex(1).square()-quad->vertex(1)(axis)*quad->vertex(1)(axis) >= radius*radius-1.e-10
- && quad->vertex(2).square()-quad->vertex(2)(axis)*quad->vertex(2)(axis) >= radius*radius-1.e-10
- && quad->vertex(3).square()-quad->vertex(3)(axis)*quad->vertex(3)(axis) >= radius*radius-1.e-10)
-
- {
- const double f = radius / std::sqrt(middle.square()-middle(axis)*middle(axis));
- for (unsigned int i=0; i<3; ++i)
- if (i!=axis)
- middle(i) *= f;
- }
- return middle;
+ // same algorithm as above
+ const unsigned int dim = 3;
+ const Point<dim> vector_from_axis = (middle-point_on_axis) -
+ ((middle-point_on_axis) * direction) * direction;
+ 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);
}
#endif
template <int dim>
void
CylinderBoundary<dim>::get_intermediate_points_between_points (
- const Point<dim> &v0, const Point<dim> &v1,
+ const Point<dim> &v0,
+ const Point<dim> &v1,
std::vector<Point<dim> > &points) const
{
const unsigned int n=points.size();
Assert(n>0, ExcInternalError());
+
// Do a simple linear interpolation
- // followed by projection. If this
- // is not sufficient, try to
- // understand the sophisticated
- // code in HyperBall later.
- Point<dim> ds = v1-v0;
- ds /= n+1;
-
- bool scale = (dim>=3
- && v0.square()-v0(axis)*v0(axis) >= radius*radius-1.e-10
- && v1.square()-v1(axis)*v1(axis) >= radius*radius-1.e-10);
-
+ // followed by projection, using
+ // the same algorithm as above
+ const Point<dim> ds = (v1-v0) / (n+1);
+
for (unsigned int i=0; i<n; ++i)
{
- if (i==0)
- points[i] = v0+ds;
- else
- points[i] = points[i-1]+ds;
+ const Point<dim> middle = v0 + (i+1)*ds;
- if (scale)
- {
- const double f = radius / std::sqrt(points[i].square()-points[i](axis)*points[i](axis));
- for (unsigned int d=0; d<dim; ++d)
- if (i!=axis)
- points[i](d) *= f;
- }
+ const Point<dim> vector_from_axis = (middle-point_on_axis) -
+ ((middle-point_on_axis) * direction) * direction;
+ if (vector_from_axis.norm() <= 1e-10 * middle.norm())
+ points[i] = middle;
+ else
+ points[i] = (vector_from_axis / vector_from_axis.norm() * radius +
+ ((middle-point_on_axis) * direction) * direction +
+ point_on_axis);
}
}
get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
{
- for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
{
- face_vertex_normals[vertex] = face->vertex(vertex);
- face_vertex_normals[vertex][axis] = 0.;
+ const Point<dim> vertex = face->vertex(v);
+
+ const Point<dim> vector_from_axis = (vertex-point_on_axis) -
+ ((vertex-point_on_axis) * direction) * direction;
+
+ face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm());
}
}