From ecbcf05ba98f885ed5d95c45565c13f68ef6ac0d Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 26 Oct 1999 16:59:44 +0000 Subject: [PATCH] Move the HyperBallBoundary from tria_boundary to tria_boundary_lib, and add a way to produce ring-shaped domains. git-svn-id: https://svn.dealii.org/trunk@1796 0785d39b-7218-0410-832d-ea1e28bc413d --- .../Attic/examples/convergence/convergence.cc | 2 +- .../error-estimation/error-estimation.cc | 2 +- .../deal.II/Attic/examples/poisson/problem.cc | 1 + deal.II/deal.II/include/grid/grid_generator.h | 74 +++++++++++- deal.II/deal.II/include/grid/tria_boundary.h | 53 -------- .../deal.II/include/grid/tria_boundary_lib.h | 113 ++++++++++++++++++ deal.II/deal.II/source/grid/grid_generator.cc | 86 ++++++++++++- deal.II/deal.II/source/grid/tria_boundary.cc | 34 +----- .../deal.II/source/grid/tria_boundary_lib.cc | 102 ++++++++++++++++ tests/big-tests/convergence/convergence.cc | 2 +- .../error-estimation/error-estimation.cc | 2 +- tests/big-tests/poisson/problem.cc | 1 + 12 files changed, 376 insertions(+), 96 deletions(-) create mode 100644 deal.II/deal.II/include/grid/tria_boundary_lib.h create mode 100644 deal.II/deal.II/source/grid/tria_boundary_lib.cc diff --git a/deal.II/deal.II/Attic/examples/convergence/convergence.cc b/deal.II/deal.II/Attic/examples/convergence/convergence.cc index afe9ac1f4b..9cdba31161 100644 --- a/deal.II/deal.II/Attic/examples/convergence/convergence.cc +++ b/deal.II/deal.II/Attic/examples/convergence/convergence.cc @@ -7,7 +7,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc b/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc index aebd9cf999..ee14debf93 100644 --- a/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc +++ b/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc @@ -7,7 +7,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/deal.II/deal.II/Attic/examples/poisson/problem.cc b/deal.II/deal.II/Attic/examples/poisson/problem.cc index 4c48dabf0b..b3a6b81192 100644 --- a/deal.II/deal.II/Attic/examples/poisson/problem.cc +++ b/deal.II/deal.II/Attic/examples/poisson/problem.cc @@ -6,6 +6,7 @@ #include "poisson.h" #include #include +#include #include diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index 088f487bd3..f425ccc041 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -7,6 +7,10 @@ #include #include +#include + + + /** * This class offers triangulations of some standard domains such as hypercubes, @@ -32,7 +36,7 @@ * \item Hyper balls: * You get the circle or ball (or generalized: hyperball) around origin * #p# and with radius #r# by calling - * #Triangulation::hyper_ball (p, r)#. The circle is triangulated + * #GridGenerator::hyper_ball (tria, p, r)#. The circle is triangulated * by five cells, the ball by seven cells. The diameter of the center cell is * chosen so that the aspect ratio of the boundary cells after one refinement * is minimized in some way. To create a hyperball in one dimension results in @@ -41,14 +45,42 @@ * Do not forget to also attach a suitable boundary approximation object * to the triangulation object you passed to this function if you later want * the triangulation to be refined at the outer boundaries. + * + * \item Hyper shell: A hyper shell is the region between two hyper + * sphere with the same origin. Therefore, it is a ring in two + * spatial dimensions. To triangulation it, call the function + * #GridGenerator::hyper_shell (tria, origin, inner_radius, + * outer_radius, N)#, where the center of the spheres as well as + * the inner and outer radius of the two spheres are given as + * shown. + * + * The parameter #N# denotes how many cells are to be used for + * this coarse triangulation. It defaults to zero, which tells + * the function to chose the number itself; this, then, is done + * such that the aspect ration of the resulting cells is as small + * as possible. However, it should be mentioned that this + * function does not work very well if the inner radius is much + * smaller than the outer radius since only one layer of cells is + * used in the radial direction. + * + * You need to attach a boundary object to the triangulation. A + * suitable boundary class is provided as #HyperSphereBoundary# + * in the library. + * + * \item Slit domain: The slit domain is a variant of the hyper cube + * domain. In two spatial dimensions, it is a square into which a slit + * is sawed; if the initial square is though to be composed of four + * smaller squares, then two of them are not connected even though + * they are neighboring each other. Analogously, into the cube in + * three spatial dimensions, a half-plane is sawed, disconnecting four + * of the eight child-cubes from one of their neighbors. * \end{itemize} * - * @author Wolfgang Bangerth, 1998, 1999 + * @author Wolfgang Bangerth, 1998, 1999. Slit domain by Stefan Nauber, 1999 */ -class GridGenerator +class GridGenerator { public: - /** * Initialize the given triangulation with a * hypercube (line in 1D, square in 2D, etc) @@ -106,11 +138,45 @@ class GridGenerator * from the middle of the top * boundary to the middle of the * area. + * + * The triangulation needs to be void + * upon calling this function. */ template static void hyper_cube_slit (Triangulation &tria, const double left = 0., const double right= 1.); + + /** + * Produce a hyper-shell, + * i.e. the space between two + * circles in two space + * dimensions and the region + * between two spheres in 3d, + * with given inner and outer + * radius and a given number of + * elements for this initial + * triangulation. If the number + * of initial cells is zero (as + * is the default), then it is + * computed adaptively such that + * the resulting elements have + * the least aspect ratio. + * + * The triangulation needs to be void + * upon calling this function. + */ + template + static void hyper_shell (Triangulation &tria, + const Point ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells = 0); + + /** + * Exception + */ + DeclException0 (ExcInvalidRadii); }; diff --git a/deal.II/deal.II/include/grid/tria_boundary.h b/deal.II/deal.II/include/grid/tria_boundary.h index bc2a07c4a6..81ce9da8a0 100644 --- a/deal.II/deal.II/include/grid/tria_boundary.h +++ b/deal.II/deal.II/include/grid/tria_boundary.h @@ -158,59 +158,6 @@ class StraightBoundary : public Boundary { -/** - * Specialisation of \Ref{Boundary}, which places the new point on - * the boundary of a ball in arbitrary dimension. It works by projecting - * the point in the middle of the old points onto the ball. The middle is - * defined as the arithmetic mean of the points. - * - * The center of the ball and its radius may be given upon construction of - * an object of this type. They default to the origin and a radius of 1.0. - * - * This class is derived from #StraightBoundary# rather than from - * #Boundary#, which would seem natural, since this way we can use the - * #StraightBoundary::in_between(neighbors)# function. - */ -template -class HyperBallBoundary : public StraightBoundary { - public: - /** - * Constructor - */ - HyperBallBoundary (const Point p=Point(), const double radius=1.0) : - center(p), radius(radius) {}; - - /** - * Refer to the general documentation of - * this class and the documentation of the - * base class. - */ - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - /** - * Refer to the general documentation of - * this class and the documentation of the - * base class. - */ - virtual Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; - - - private: - /** - * Center point of the hyperball. - */ - const Point center; - - /** - * Radius of the hyperball. - */ - const double radius; -}; - - - /*---------------------------- boundary-function.h ---------------------------*/ /* end of #ifndef __tria_boundary_H */ diff --git a/deal.II/deal.II/include/grid/tria_boundary_lib.h b/deal.II/deal.II/include/grid/tria_boundary_lib.h new file mode 100644 index 0000000000..5f4839d907 --- /dev/null +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -0,0 +1,113 @@ +/*---------------------------- tria_boundary_lib.h ---------------------------*/ +/* $Id$ */ +#ifndef __tria_boundary_lib_H +#define __tria_boundary_lib_H +/*---------------------------- tria_boundary_lib.h ---------------------------*/ + + +#include + + + + +/** + * Specialisation of \Ref{Boundary}, which places the new point on + * the boundary of a ball in arbitrary dimension. It works by projecting + * the point in the middle of the old points onto the ball. The middle is + * defined as the arithmetic mean of the points. + * + * The center of the ball and its radius may be given upon construction of + * an object of this type. They default to the origin and a radius of 1.0. + * + * This class is derived from #StraightBoundary# rather than from + * #Boundary#, which would seem natural, since this way we can use the + * #StraightBoundary::in_between(neighbors)# function. + * + * @author Wolfgang Bangerth, 1998 + */ +template +class HyperBallBoundary : public StraightBoundary { + public: + /** + * Constructor + */ + HyperBallBoundary (const Point p=Point(), const double radius=1.0) : + center(p), radius(radius) {}; + + /** + * Refer to the general documentation of + * this class and the documentation of the + * base class. + */ + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + + /** + * Refer to the general documentation of + * this class and the documentation of the + * base class. + */ + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; + + + private: + /** + * Center point of the hyperball. + */ + const Point center; + + /** + * Radius of the hyperball. + */ + const double radius; +}; + + + + +/** + * Class describing the boundaries of a hyper shell. Only the center + * of the two spheres needs to be given, the radii of inner and outer + * sphere are computed automatically upon calling one of the virtual + * functions. + * + * @author Wolfgang Bangerth, 1999 + */ +template +class HyperShellBoundary : public StraightBoundary +{ + public: + /** + * Constructor. The center of the + * spheres defaults to the + * origin. + */ + HyperShellBoundary (const Point ¢er = Point()); + + /** + * Construct a new point on a line. + */ + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + + /** + * Construct a new point on a quad. + */ + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; + + private: + /** + * Store the center of the spheres. + */ + const Point center; +}; + + + + +/*---------------------------- tria_boundary_lib.h ---------------------------*/ +/* end of #ifndef __tria_boundary_lib_H */ +#endif +/*---------------------------- tria_boundary_lib.h ---------------------------*/ diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index 6b52fed0f8..bc023fc627 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -8,7 +8,6 @@ - #if deal_II_dimension == 1 template <> @@ -53,6 +52,17 @@ void GridGenerator::hyper_ball<> (Triangulation<1> &, Assert (false, ExcInternalError()); }; + + +template <> +void GridGenerator::hyper_shell<> (Triangulation<1> &, + const Point<1> &, + const double, + const double, + const unsigned int) { + Assert (false, ExcInternalError()); +}; + #endif @@ -182,6 +192,67 @@ void GridGenerator::hyper_ball<> (Triangulation<2> &tria, SubCellData()); // no boundary information }; + + +template <> +void GridGenerator::hyper_shell<> (Triangulation<2> &tria, + const Point<2> ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells) +{ + Assert ((inner_radius > 0) && (inner_radius < outer_radius), + ExcInvalidRadii ()); + + const double pi = 3.1415926536; + // determine the number of cells + // for the grid. if not provided by + // the user determine it such that + // the length of each cell on the + // median (in the middle between + // the two circles) is equal to its + // radial extent (which is the + // difference between the two + // radii) + const unsigned int N = (n_cells == 0 ? + static_cast + (ceil((2*pi* (outer_radius + inner_radius)/2) / + (outer_radius - inner_radius))) : + n_cells); + + // set up N vertices on the + // outer and N vertices on + // the inner circle. the + // first N ones are on the + // outer one, and all are + // numbered counter-clockwise + vector > vertices(2*N); + for (unsigned int i=0; i( cos(2*pi * i/N), + sin(2*pi * i/N)) * outer_radius; + vertices[i+N] = vertices[i] * (inner_radius/outer_radius); + + vertices[i] += center; + vertices[i+N] += center; + }; + + vector > cells (N, CellData<2>()); + + for (unsigned int i=0; i void GridGenerator::hyper_cube_slit<> (Triangulation<3> &, const double, const double) { - Assert (false, ExcInternalError()); + Assert (false, ExcNotImplemented()); }; @@ -345,6 +416,17 @@ void GridGenerator::hyper_ball<> (Triangulation<3> &tria, +template <> +void GridGenerator::hyper_shell<> (Triangulation<3> &, + const Point<3> &, + const double , + const double , + const unsigned int ) +{ + Assert (false, ExcNotImplemented()); +}; + + #endif diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc index 7a076d862e..81280b632c 100644 --- a/deal.II/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/deal.II/source/grid/tria_boundary.cc @@ -66,40 +66,8 @@ StraightBoundary::get_new_point_on_quad (const typename Triangulation: -template -Point -HyperBallBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - Point middle = StraightBoundary::get_new_point_on_line (line); - - middle -= center; - // project to boundary - middle *= radius / sqrt(middle.square()); - - middle += center; - return middle; -}; - - - -template -Point -HyperBallBoundary::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const -{ - Point middle = StraightBoundary::get_new_point_on_quad (quad); - - middle -= center; - // project to boundary - middle *= radius / sqrt(middle.square()); - - middle += center; - return middle; -}; - - - // explicit instantiations template class Boundary; template class StraightBoundary; -template class HyperBallBoundary; + diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc new file mode 100644 index 0000000000..6266ffaeab --- /dev/null +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -0,0 +1,102 @@ +/* $Id$ */ +/* Copyright W. Bangerth, University of Heidelberg, 1998 */ + +#include +#include +#include +#include +#include + + + + +template +Point +HyperBallBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + Point middle = StraightBoundary::get_new_point_on_line (line); + + middle -= center; + // project to boundary + middle *= radius / sqrt(middle.square()); + + middle += center; + return middle; +}; + + + +template +Point +HyperBallBoundary:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + Point middle = StraightBoundary::get_new_point_on_quad (quad); + + middle -= center; + // project to boundary + middle *= radius / sqrt(middle.square()); + + middle += center; + return middle; +}; + + + + + +template +HyperShellBoundary::HyperShellBoundary (const Point ¢er) : + center (center) +{}; + + +template +Point +HyperShellBoundary:: +get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + const Point middle = StraightBoundary::get_new_point_on_line (line); + // compute the position of the points relative to the origin + const Point middle_relative = middle - center, + vertex_relative = line->vertex(0) - center; + + // take vertex(0) to gauge the + // radius corresponding to the line + // under consideration + const double radius = sqrt(vertex_relative.square()); + + // scale and shift back to the + // original coordinate system + return (middle_relative * (radius / sqrt(middle_relative.square()))) + center; +}; + + + +template +Point +HyperShellBoundary:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + const Point middle = StraightBoundary::get_new_point_on_quad (quad); + // compute the position of the points relative to the origin + const Point middle_relative = middle - center, + vertex_relative = quad->vertex(0) - center; + + // take vertex(0) to gauge the + // radius corresponding to the line + // under consideration + const double radius = sqrt(vertex_relative.square()); + + // scale and shift back to the + // original coordinate system + return (middle_relative * (radius / sqrt(middle_relative.square()))) + center; +}; + + + + + +// explicit instantiations +template class HyperBallBoundary; +template class HyperShellBoundary; diff --git a/tests/big-tests/convergence/convergence.cc b/tests/big-tests/convergence/convergence.cc index afe9ac1f4b..9cdba31161 100644 --- a/tests/big-tests/convergence/convergence.cc +++ b/tests/big-tests/convergence/convergence.cc @@ -7,7 +7,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/tests/big-tests/error-estimation/error-estimation.cc b/tests/big-tests/error-estimation/error-estimation.cc index aebd9cf999..ee14debf93 100644 --- a/tests/big-tests/error-estimation/error-estimation.cc +++ b/tests/big-tests/error-estimation/error-estimation.cc @@ -7,7 +7,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/tests/big-tests/poisson/problem.cc b/tests/big-tests/poisson/problem.cc index 4c48dabf0b..b3a6b81192 100644 --- a/tests/big-tests/poisson/problem.cc +++ b/tests/big-tests/poisson/problem.cc @@ -6,6 +6,7 @@ #include "poisson.h" #include #include +#include #include -- 2.39.5