]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Handle arbitrarily oriented and located cylinders.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 8 Mar 2007 19:42:35 +0000 (19:42 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 8 Mar 2007 19:42:35 +0000 (19:42 +0000)
git-svn-id: https://svn.dealii.org/trunk@14554 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/tria_boundary_lib.h
deal.II/deal.II/source/grid/tria_boundary_lib.cc
deal.II/doc/news/changes.html

index 862e63dbcbba479ccc0fb67d01eab332e8a84b5f..6ad89a2349a81c8de7e19f8dbd3147ba734f3989 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -22,10 +22,13 @@ DEAL_II_NAMESPACE_OPEN
 /**
  * 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
@@ -37,7 +40,7 @@ DEAL_II_NAMESPACE_OPEN
  *
  * @ingroup boundary
  *
- * @author Guido Kanschat, 2001
+ * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007
  */
 template <int dim>
 class CylinderBoundary : public StraightBoundary<dim>
@@ -52,8 +55,24 @@ 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
@@ -132,16 +151,14 @@ class CylinderBoundary : public StraightBoundary<dim>
     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:
 
@@ -158,7 +175,14 @@ class CylinderBoundary : public StraightBoundary<dim>
                                      * 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);
 };
 
 
index fc58f0772ab1c10169704d5b20a50e17add94351..a1bc4609a3d8d67ee1ce4defb4b1e99a24865450 100644 (file)
@@ -2,7 +2,7 @@
 //    $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);
 }
 
 
@@ -57,21 +95,18 @@ Point<3>
 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
@@ -104,37 +139,30 @@ CylinderBoundary<dim>::get_intermediate_points_on_line (
 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);
     }
 }
 
@@ -205,10 +233,14 @@ CylinderBoundary<dim>::
 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());
     }
 }
 
index deb41f8b5570eab31624ef8ea6cca7348b4d7e1f..ef97f3433d762b21a5c3659c3b8da3f17de87fd3 100644 (file)
@@ -1023,6 +1023,14 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+   <li> <p>Improved: The
+       <code>CylinderBoundary</code> class can now describe the
+       boundaries of arbitrarily oriented cylinders, not only does
+       oriented parallel to the axes and going through the origin.
+       <br>
+       (WB 2007/03/08)
+       </p>
+
    <li> <p>Improved: The
        <code>GridRefinement::refine_and_coarsen_fixed_number</code> and
        <code>GridRefinement::refine_and_coarsen_fixed_fraction</code> functions

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.