From 0b8a86e4c529904f963459aecec136d7d29a01b7 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 8 Aug 2000 11:42:59 +0000 Subject: [PATCH] Implement HalfHyperShellBoundary and the corresponding grid generation function. git-svn-id: https://svn.dealii.org/trunk@3234 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_generator.h | 53 ++++++++++- .../deal.II/include/grid/tria_boundary_lib.h | 55 ++++++++++- deal.II/deal.II/source/grid/grid_generator.cc | 84 ++++++++++++++++- .../deal.II/source/grid/tria_boundary_lib.cc | 92 +++++++++++++++++++ 4 files changed, 278 insertions(+), 6 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index f37d336a49..e46d51dbaa 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -41,7 +41,7 @@ template class Triangulation; * $[-1,1]^2 - [0,1]^2$. To create a hyper-L in one dimension results in * an error. The function is also implemented for three space dimensions. * - * @item Hyper balls: + * @item Hyper ball: * You get the circle or ball (or generalized: hyperball) around origin * @p{p} and with radius @p{r} by calling * @ref{GridGenerator}@p{::hyper_ball (tria, p, r)}. The circle is triangulated @@ -75,7 +75,13 @@ template class Triangulation; * suitable boundary class is provided as @ref{HyperSphereBoundary} * in the library. * - * @item Slit domain: The slit domain is a variant of the hyper cube + * @item Half shells: This is a variant of the preceeding domain, + * where only one half of the shell is taken, defined by the first + * coordinate having non-negative values. This class is useful in + * computations with rotational symmetry, as then a 3d shell can be + * modeled by a 2d half shell. + * + * @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 @@ -221,6 +227,49 @@ class GridGenerator const double outer_radius, const unsigned int n_cells = 0); + + /** + * Produce a half 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. However, + * opposed to the previous + * function, it does not produce + * a whole shell, but only one + * half of it, namely that part + * for which the first component + * is restricted to non-negative + * values. The purpose of this + * class is to enable + * computations for solutions + * which have rotational + * symmetry, in which case the + * half shell in 2d represents a + * shell in 3d. + * + * 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 half_hyper_shell (Triangulation &tria, + const Point ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells = 0); + /** * Exception */ 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 70a8695cd6..501272620b 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -33,7 +33,8 @@ * @author Wolfgang Bangerth, 1998 */ template -class HyperBallBoundary : public StraightBoundary { +class HyperBallBoundary : public StraightBoundary +{ public: /** * Constructor @@ -80,16 +81,20 @@ class HyperBallBoundary : public StraightBoundary { }; + /** * Variant of @p{HyperBallBoundary} which denotes a half hyper ball * where the first coordinate is restricted to the range $x>=0$ (or * $x>=center(0)$). In two dimensions, this equals the right half - * circle, in three space dimensions it is a half ball. + * circle, in three space dimensions it is a half ball. This class + * might be useful for computations with rotational symmetry, where + * one dimension is the radius from the axis of rotation. * * @author Wolfgang Bangerth, 1999 */ template -class HalfHyperBallBoundary : public HyperBallBoundary { +class HalfHyperBallBoundary : public HyperBallBoundary +{ public: /** * Constructor @@ -115,6 +120,7 @@ class HalfHyperBallBoundary : public HyperBallBoundary { }; + /** * 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 @@ -154,4 +160,47 @@ class HyperShellBoundary : public StraightBoundary }; + +/** + * Variant of @p{HyperShellBoundary} which denotes a half hyper shell + * where the first coordinate is restricted to the range $x>=0$ (or + * $x>=center(0)$). In two dimensions, this equals the right half arc, + * in three space dimensions it is a half shell. This class might be + * useful for computations with rotational symmetry, where one + * dimension is the radius from the axis of rotation. + * + * @author Wolfgang Bangerth, 2000 + */ +template +class HalfHyperShellBoundary : public HyperShellBoundary +{ + public: + /** + * Constructor. The center of the + * spheres defaults to the + * origin. + */ + HalfHyperShellBoundary (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; +}; + + + #endif diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index b44f8eadb0..8356e2f4a2 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -35,6 +35,7 @@ void GridGenerator::hyper_cube<> (Triangulation<1> &tria, }; + template <> void GridGenerator::hyper_cube_slit<> (Triangulation<1> &, const double, @@ -43,6 +44,7 @@ void GridGenerator::hyper_cube_slit<> (Triangulation<1> &, }; + template <> void GridGenerator::hyper_L<> (Triangulation<1> &, const double, @@ -51,6 +53,7 @@ void GridGenerator::hyper_L<> (Triangulation<1> &, }; + template <> void GridGenerator::hyper_ball<> (Triangulation<1> &, const Point<1> &, @@ -59,6 +62,7 @@ void GridGenerator::hyper_ball<> (Triangulation<1> &, }; + template <> void GridGenerator::hyper_shell<> (Triangulation<1> &, const Point<1> &, @@ -71,6 +75,7 @@ void GridGenerator::hyper_shell<> (Triangulation<1> &, #endif + #if deal_II_dimension == 2 template <> @@ -136,7 +141,9 @@ void GridGenerator::enclosed_hyper_cube (Triangulation<2> &tria, cells, SubCellData()); // no boundary information } - + + + template <> void GridGenerator::hyper_cube_slit<> (Triangulation<2> &tria, const double left, @@ -169,6 +176,7 @@ void GridGenerator::hyper_cube_slit<> (Triangulation<2> &tria, }; + template <> void GridGenerator::hyper_L<> (Triangulation<2> &tria, const double a, @@ -201,6 +209,7 @@ void GridGenerator::hyper_L<> (Triangulation<2> &tria, }; + template <> void GridGenerator::hyper_ball<> (Triangulation<2> &tria, const Point<2> &p, @@ -238,6 +247,7 @@ void GridGenerator::hyper_ball<> (Triangulation<2> &tria, }; + template <> void GridGenerator::hyper_shell<> (Triangulation<2> &tria, const Point<2> ¢er, @@ -297,6 +307,76 @@ void GridGenerator::hyper_shell<> (Triangulation<2> &tria, }; + +template <> +void GridGenerator::half_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.14159265359; + // 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((pi* (outer_radius + inner_radius)/2) / + (outer_radius - inner_radius))) : + n_cells); + + // set up N+1 vertices on the + // outer and N+1 vertices on + // the inner circle. the + // first N+1 ones are on the + // outer one, and all are + // numbered counter-clockwise + vector > vertices(2*(N+1)); + for (unsigned int i=0; i<=N; ++i) + { + // enforce that the x-coordinates + // of the first and last point of + // each half-circle are exactly + // zero (contrary to what we may + // compute using the imprecise + // value of pi) + vertices[i] = Point<2>( ( (i==0) || (i==N) ? + 0 : + cos(pi * i/N - pi/2) ), + sin(pi * i/N - pi/2)) * outer_radius; + vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius); + + vertices[i] += center; + vertices[i+N+1] += center; + }; + + + vector > cells (N, CellData<2>()); + + for (unsigned int i=0; i (Triangulation<3> &tria, }; + template <> void GridGenerator::hyper_cube_slit<> (Triangulation<3> &, const double, @@ -335,6 +416,7 @@ void GridGenerator::hyper_cube_slit<> (Triangulation<3> &, }; + template<> void GridGenerator::enclosed_hyper_cube (Triangulation<3> &tria, const double l, 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 945d4c6fab..67eff897c5 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -26,6 +26,7 @@ HyperBallBoundary::HyperBallBoundary (const Point p, {}; + template Point HyperBallBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const @@ -41,6 +42,7 @@ HyperBallBoundary::get_new_point_on_line (const typename Triangulation }; + #if deal_II_dimension == 1 template <> @@ -55,6 +57,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const #endif + template Point HyperBallBoundary:: @@ -71,6 +74,7 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c }; + template Point HyperBallBoundary::get_center () const @@ -79,6 +83,7 @@ HyperBallBoundary::get_center () const }; + template double HyperBallBoundary::get_radius () const @@ -87,6 +92,9 @@ HyperBallBoundary::get_radius () const }; +/* ---------------------------------------------------------------------- */ + + template HalfHyperBallBoundary::HalfHyperBallBoundary (const Point center, const double radius) : @@ -94,6 +102,7 @@ HalfHyperBallBoundary::HalfHyperBallBoundary (const Point center, {}; + template Point HalfHyperBallBoundary:: @@ -107,6 +116,7 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c }; + #if deal_II_dimension == 1 template <> @@ -121,6 +131,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const #endif + template Point HalfHyperBallBoundary:: @@ -134,12 +145,18 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c }; + +/* ---------------------------------------------------------------------- */ + + + template HyperShellBoundary::HyperShellBoundary (const Point ¢er) : center (center) {}; + template Point HyperShellBoundary:: @@ -161,6 +178,7 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c }; + #if deal_II_dimension == 1 template <> @@ -175,6 +193,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const #endif + template Point HyperShellBoundary:: @@ -196,7 +215,80 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c }; +/* ---------------------------------------------------------------------- */ + + + + +template +HalfHyperShellBoundary::HalfHyperShellBoundary (const Point ¢er) : + HyperShellBoundary (center) +{}; + + + +template +Point +HalfHyperShellBoundary:: +get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + // first check whether the two end + // points of the line are on the + // axis of symmetry. if so, then + // return the mid point + if ((line->vertex(0)(0) == center(0)) && + (line->vertex(1)(0) == center(0))) + return (line->vertex(0) + line->vertex(1))/2; + + + // otherwise we are on the outer or + // inner part of the shell. proceed + // as in the base class + return HyperShellBoundary::get_new_point_on_line (line); +}; + + + +#if deal_II_dimension == 1 + +template <> +Point<1> +HalfHyperShellBoundary<1>:: +get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const +{ + Assert (false, ExcInternalError()); + return Point<1>(); +}; + +#endif + + + +template +Point +HalfHyperShellBoundary:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + // same thing as for the new point + // on the line + if ((quad->vertex(0)(0) == center(0)) && + (quad->vertex(1)(0) == center(0)) && + (quad->vertex(2)(0) == center(0)) && + (quad->vertex(3)(0) == center(0))) + return (quad->vertex(0) + quad->vertex(1) + + quad->vertex(2) + quad->vertex(3) )/4; + + + // otherwise we are on the outer or + // inner part of the shell. proceed + // as in the base class + return HyperShellBoundary::get_new_point_on_quad (quad); +}; + + + // explicit instantiations template class HyperBallBoundary; template class HalfHyperBallBoundary; template class HyperShellBoundary; +template class HalfHyperShellBoundary; -- 2.39.5