From a5f3b79f089478d25aaa1b1bd2150028fe467e6a Mon Sep 17 00:00:00 2001 From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Fri, 15 Oct 2004 07:00:24 +0000 Subject: [PATCH] Extension of CylinderBoundary class to represent circular tubes also along the y- and z-axis. git-svn-id: https://svn.dealii.org/trunk@9710 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/grid/tria_boundary_lib.h | 32 ++++++++++--- .../deal.II/source/grid/tria_boundary_lib.cc | 47 ++++++++++--------- deal.II/doc/news/c-5.0.html | 9 ++++ 3 files changed, 61 insertions(+), 27 deletions(-) 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 8b2c7c9216..4d0676c3ad 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -20,10 +20,11 @@ /** * Boundary object for the hull of a cylinder. In three dimensions, - * points are projected on a circular tube along the @p x-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. + * 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). * * This class was developed to be used in conjunction with the * @p cylinder function of GridGenerator. It should be used for @@ -40,9 +41,16 @@ class CylinderBoundary : public StraightBoundary<dim> { public: /** - * Constructor + * Constructor. Per default + * 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); + CylinderBoundary (const double radius = 1.0, + const unsigned int axis = 0); /** * Refer to the general documentation of @@ -120,6 +128,18 @@ 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. + */ + const unsigned int axis; + private: /** 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 b61bb505d6..97b6a9bdf2 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) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -21,8 +21,10 @@ template <int dim> -CylinderBoundary<dim>::CylinderBoundary (const double radius) : - radius(radius) +CylinderBoundary<dim>::CylinderBoundary (const double radius, + const unsigned int axis) : + radius(radius), + axis(axis) {} @@ -34,12 +36,13 @@ CylinderBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>: Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line); // project to boundary if (dim>=3 - && line->vertex(0).square()-line->vertex(0)(0)*line->vertex(0)(0) >= radius*radius-1.e-12 - && line->vertex(1).square()-line->vertex(1)(0)*line->vertex(1)(0) >= radius*radius-1.e-12) + && 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(0)*middle(0)); - for (unsigned int i=1;i<dim;++i) - middle(i) *= f; + 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; } @@ -55,15 +58,16 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad); // project to boundary - if (quad->vertex(0).square()-quad->vertex(0)(0)*quad->vertex(0)(0) >= radius*radius-1.e-12 - && quad->vertex(1).square()-quad->vertex(1)(0)*quad->vertex(1)(0) >= radius*radius-1.e-12 - && quad->vertex(2).square()-quad->vertex(2)(0)*quad->vertex(2)(0) >= radius*radius-1.e-12 - && quad->vertex(3).square()-quad->vertex(3)(0)*quad->vertex(3)(0) >= radius*radius-1.e-12) + 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(0)*middle(0)); - for (unsigned int i=1;i<3;++i) - middle(i) *= f; + 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; } @@ -112,8 +116,8 @@ CylinderBoundary<dim>::get_intermediate_points_between_points ( ds /= n+1; bool scale = (dim>=3 - && v0.square()-v0(0)*v0(0) >= radius*radius-1.e-12 - && v1.square()-v1(0)*v1(0) >= radius*radius-1.e-12); + && v0.square()-v0(axis)*v0(axis) >= radius*radius-1.e-10 + && v1.square()-v1(axis)*v1(axis) >= radius*radius-1.e-10); for (unsigned int i=0; i<n; ++i) { @@ -124,9 +128,10 @@ CylinderBoundary<dim>::get_intermediate_points_between_points ( if (scale) { - const double f = radius / std::sqrt(points[i].square()-points[i](0)*points[i](0)); - for (unsigned int d=1;d<dim;++d) - points[i](d) *= f; + 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; } } } @@ -196,7 +201,7 @@ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face, for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex) { face_vertex_normals[vertex] = face->vertex(vertex); - face_vertex_normals[vertex][0] = 0.; + face_vertex_normals[vertex][axis] = 0.; } } diff --git a/deal.II/doc/news/c-5.0.html b/deal.II/doc/news/c-5.0.html index f78d446a8c..38f023982f 100644 --- a/deal.II/doc/news/c-5.0.html +++ b/deal.II/doc/news/c-5.0.html @@ -359,6 +359,15 @@ inconvenience this causes. <h3>deal.II</h3> <ol> + <li> <p> + Improved: The <code class="class">CylinderBoundary</code> + represented the hull of a circular tubes along the x-axis. It + is now extended to allow for circular tubes also along the y- + and z-axis. + <br> + (RH 2004/10/15) + </p> + <li> <p> Fixed: The <code class="class">ConstraintMatrix</code> class had some algorithms that were linear in the number of constraints. Since these -- 2.39.5