#include <grid/tria_accessor.h>
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
+#include <grid/tria_boundary_lib.h>
#include <grid/dof_constraints.h>
#include <grid/grid_generator.h>
#include <base/function.h>
#include <grid/tria_accessor.h>
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
+#include <grid/tria_boundary_lib.h>
#include <grid/dof_constraints.h>
#include <grid/grid_generator.h>
#include <base/function.h>
#include "poisson.h"
#include <lac/vector.h>
#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
#include <basic/data_out.h>
#include <base/forward-declarations.h>
#include <basic/forward-declarations.h>
+#include <base/exceptions.h>
+
+
+
/**
* This class offers triangulations of some standard domains such as hypercubes,
* \item Hyper balls:
* You get the circle or ball (or generalized: hyperball) around origin
* #p# and with radius #r# by calling
- * #Triangulation<dim>::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
* 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)
* from the middle of the top
* boundary to the middle of the
* area.
+ *
+ * The triangulation needs to be void
+ * upon calling this function.
*/
template <int dim>
static void hyper_cube_slit (Triangulation<dim> &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 <int dim>
+ static void hyper_shell (Triangulation<dim> &tria,
+ const Point<dim> ¢er,
+ const double inner_radius,
+ const double outer_radius,
+ const unsigned int n_cells = 0);
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInvalidRadii);
};
-/**
- * Specialisation of \Ref{Boundary}<dim>, 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<dim>::in_between(neighbors)# function.
- */
-template <int dim>
-class HyperBallBoundary : public StraightBoundary<dim> {
- public:
- /**
- * Constructor
- */
- HyperBallBoundary (const Point<dim> p=Point<dim>(), 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<dim>
- get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
-
- /**
- * Refer to the general documentation of
- * this class and the documentation of the
- * base class.
- */
- virtual Point<dim>
- get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
-
-
- private:
- /**
- * Center point of the hyperball.
- */
- const Point<dim> center;
-
- /**
- * Radius of the hyperball.
- */
- const double radius;
-};
-
-
-
/*---------------------------- boundary-function.h ---------------------------*/
/* end of #ifndef __tria_boundary_H */
--- /dev/null
+/*---------------------------- tria_boundary_lib.h ---------------------------*/
+/* $Id$ */
+#ifndef __tria_boundary_lib_H
+#define __tria_boundary_lib_H
+/*---------------------------- tria_boundary_lib.h ---------------------------*/
+
+
+#include <grid/tria_boundary.h>
+
+
+
+
+/**
+ * Specialisation of \Ref{Boundary}<dim>, 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<dim>::in_between(neighbors)# function.
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class HyperBallBoundary : public StraightBoundary<dim> {
+ public:
+ /**
+ * Constructor
+ */
+ HyperBallBoundary (const Point<dim> p=Point<dim>(), 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<dim>
+ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+
+ /**
+ * Refer to the general documentation of
+ * this class and the documentation of the
+ * base class.
+ */
+ virtual Point<dim>
+ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+
+
+ private:
+ /**
+ * Center point of the hyperball.
+ */
+ const Point<dim> 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 <int dim>
+class HyperShellBoundary : public StraightBoundary<dim>
+{
+ public:
+ /**
+ * Constructor. The center of the
+ * spheres defaults to the
+ * origin.
+ */
+ HyperShellBoundary (const Point<dim> ¢er = Point<dim>());
+
+ /**
+ * Construct a new point on a line.
+ */
+ virtual Point<dim>
+ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+
+ /**
+ * Construct a new point on a quad.
+ */
+ virtual Point<dim>
+ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+
+ private:
+ /**
+ * Store the center of the spheres.
+ */
+ const Point<dim> center;
+};
+
+
+
+
+/*---------------------------- tria_boundary_lib.h ---------------------------*/
+/* end of #ifndef __tria_boundary_lib_H */
+#endif
+/*---------------------------- tria_boundary_lib.h ---------------------------*/
-
#if deal_II_dimension == 1
template <>
Assert (false, ExcInternalError());
};
+
+
+template <>
+void GridGenerator::hyper_shell<> (Triangulation<1> &,
+ const Point<1> &,
+ const double,
+ const double,
+ const unsigned int) {
+ Assert (false, ExcInternalError());
+};
+
#endif
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<unsigned int>
+ (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<Point<2> > vertices(2*N);
+ for (unsigned int i=0; i<N; ++i)
+ {
+ vertices[i] = Point<2>( 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<CellData<2> > cells (N, CellData<2>());
+
+ for (unsigned int i=0; i<N; ++i)
+ {
+ cells[i].vertices[0] = i;
+ cells[i].vertices[1] = (i+1)%N;
+ cells[i].vertices[2] = N+((i+1)%N);
+ cells[i].vertices[3] = N+i;
+
+ cells[i].material_id = 0;
+ };
+
+ tria.create_triangulation (vertices, cells, SubCellData());
+};
+
+
#endif
void GridGenerator::hyper_cube_slit<> (Triangulation<3> &,
const double,
const double) {
- Assert (false, ExcInternalError());
+ Assert (false, ExcNotImplemented());
};
+template <>
+void GridGenerator::hyper_shell<> (Triangulation<3> &,
+ const Point<3> &,
+ const double ,
+ const double ,
+ const unsigned int )
+{
+ Assert (false, ExcNotImplemented());
+};
+
+
#endif
-template <int dim>
-Point<dim>
-HyperBallBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-{
- Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
-
- middle -= center;
- // project to boundary
- middle *= radius / sqrt(middle.square());
-
- middle += center;
- return middle;
-};
-
-
-
-template <int dim>
-Point<dim>
-HyperBallBoundary<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-{
- Point<dim> middle = StraightBoundary<dim>::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<deal_II_dimension>;
template class StraightBoundary<deal_II_dimension>;
-template class HyperBallBoundary<deal_II_dimension>;
+
--- /dev/null
+/* $Id$ */
+/* Copyright W. Bangerth, University of Heidelberg, 1998 */
+
+#include <grid/tria_boundary_lib.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <cmath>
+
+
+
+
+template <int dim>
+Point<dim>
+HyperBallBoundary<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+ Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+
+ middle -= center;
+ // project to boundary
+ middle *= radius / sqrt(middle.square());
+
+ middle += center;
+ return middle;
+};
+
+
+
+template <int dim>
+Point<dim>
+HyperBallBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+{
+ Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad (quad);
+
+ middle -= center;
+ // project to boundary
+ middle *= radius / sqrt(middle.square());
+
+ middle += center;
+ return middle;
+};
+
+
+
+
+
+template <int dim>
+HyperShellBoundary<dim>::HyperShellBoundary (const Point<dim> ¢er) :
+ center (center)
+{};
+
+
+template <int dim>
+Point<dim>
+HyperShellBoundary<dim>::
+get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+{
+ const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+ // compute the position of the points relative to the origin
+ const Point<dim> 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 <int dim>
+Point<dim>
+HyperShellBoundary<dim>::
+get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+{
+ const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_quad (quad);
+ // compute the position of the points relative to the origin
+ const Point<dim> 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<deal_II_dimension>;
+template class HyperShellBoundary<deal_II_dimension>;
#include <grid/tria_accessor.h>
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
+#include <grid/tria_boundary_lib.h>
#include <grid/dof_constraints.h>
#include <grid/grid_generator.h>
#include <base/function.h>
#include <grid/tria_accessor.h>
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
+#include <grid/tria_boundary_lib.h>
#include <grid/dof_constraints.h>
#include <grid/grid_generator.h>
#include <base/function.h>
#include "poisson.h"
#include <lac/vector.h>
#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
#include <basic/data_out.h>