From: bangerth Date: Thu, 8 Mar 2007 19:42:35 +0000 (+0000) Subject: Handle arbitrarily oriented and located cylinders. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=307c49642109c73f45433010685d328e78185f1a;p=dealii-svn.git Handle arbitrarily oriented and located cylinders. git-svn-id: https://svn.dealii.org/trunk@14554 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria_boundary_lib.h b/deal.II/deal.II/include/grid/tria_boundary_lib.h index 862e63dbcb..6ad89a2349 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -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 x-, - * y- or z-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). + * y- or z-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 class CylinderBoundary : public StraightBoundary @@ -52,8 +55,24 @@ class CylinderBoundary : public StraightBoundary * 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 direction, + const Point point_on_axis); /** * Refer to the general documentation of @@ -132,16 +151,14 @@ class CylinderBoundary : public StraightBoundary const double radius; /** - * Denotes the axis of the - * circular - * tube. axis=0, - * axis=1 and - * axis=2 denote the - * x-axis, - * y-axis and - * z-axis, respectively. + * The direction vector of the axis. + */ + const Point direction; + + /** + * An arbitrary point on the axis. */ - const unsigned int axis; + const Point point_on_axis; private: @@ -158,7 +175,14 @@ class CylinderBoundary : public StraightBoundary * base class. */ void get_intermediate_points_between_points (const Point &p0, const Point &p1, - std::vector > &points) const; + std::vector > &points) const; + + /** + * Given a number for the axis, + * return a vector that denotes + * this direction. + */ + static Point get_axis_vector (const unsigned int axis); }; diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc index fc58f0772a..a1bc4609a3 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -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,31 +22,69 @@ DEAL_II_NAMESPACE_OPEN + + + template CylinderBoundary::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()) {} +template +CylinderBoundary::CylinderBoundary (const double radius, + const Point direction, + const Point point_on_axis) + : + radius(radius), + direction (direction / direction.norm()), + point_on_axis (point_on_axis) +{} + template Point -CylinderBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +CylinderBoundary::get_axis_vector (const unsigned int axis) { - Point middle = StraightBoundary::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 axis_vector; + axis_vector[axis] = 1; + return axis_vector; +} + + + +template +Point +CylinderBoundary:: +get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + // compute a proposed new point + const Point middle = StraightBoundary::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 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 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::get_intermediate_points_on_line ( template void CylinderBoundary::get_intermediate_points_between_points ( - const Point &v0, const Point &v1, + const Point &v0, + const Point &v1, std::vector > &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 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 ds = (v1-v0) / (n+1); + for (unsigned int i=0; i 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 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:: get_normals_at_vertices (const typename Triangulation::face_iterator &face, typename Boundary::FaceVertexNormals &face_vertex_normals) const { - for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + for (unsigned int v=0; v::vertices_per_face; ++v) { - face_vertex_normals[vertex] = face->vertex(vertex); - face_vertex_normals[vertex][axis] = 0.; + const Point vertex = face->vertex(v); + + const Point vector_from_axis = (vertex-point_on_axis) - + ((vertex-point_on_axis) * direction) * direction; + + face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm()); } } diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index deb41f8b55..ef97f3433d 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -1023,6 +1023,14 @@ inconvenience this causes.

deal.II

    +
  1. Improved: The + CylinderBoundary class can now describe the + boundaries of arbitrarily oriented cylinders, not only does + oriented parallel to the axes and going through the origin. +
    + (WB 2007/03/08) +

    +
  2. Improved: The GridRefinement::refine_and_coarsen_fixed_number and GridRefinement::refine_and_coarsen_fixed_fraction functions