From e6a6145cc09abf02548d022085ee4cef6086c2db Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 18 Feb 1999 14:36:53 +0000 Subject: [PATCH] Replace the old scheme to describe the boundary, which was really archaic, by a new and much more flexible one. git-svn-id: https://svn.dealii.org/trunk@836 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/tria_boundary.h | 175 ++++++++++--------- deal.II/deal.II/source/fe/fe.cc | 1 + deal.II/deal.II/source/grid/tria.cc | 36 +--- deal.II/deal.II/source/grid/tria_boundary.cc | 71 +++++++- 4 files changed, 168 insertions(+), 115 deletions(-) diff --git a/deal.II/deal.II/include/grid/tria_boundary.h b/deal.II/deal.II/include/grid/tria_boundary.h index 78a3696427..ec61dc0a1c 100644 --- a/deal.II/deal.II/include/grid/tria_boundary.h +++ b/deal.II/deal.II/include/grid/tria_boundary.h @@ -10,52 +10,8 @@ #include +template class Triangulation; -/** - * Workaround for a bug in egcs snapshot 1998/08/03. - */ -template -struct BoundaryHelper; - -/** - * Workaround for a bug in egcs snapshot 1998/08/03. - */ -template <> -struct BoundaryHelper<1> { - /** - * Declare a data type for the derived - * classes. - * - * actually, this does not make much - * sense, but declaring a zero-sized - * array is forbidden nowadays. - */ - typedef const Point<1> *PointArray[1]; -}; - -/** - * Workaround for a bug in egcs snapshot 1998/08/03. - */ -template <> -struct BoundaryHelper<2> { - /** - * Declare a data type for the derived - * classes. - */ - typedef const Point<2> *PointArray[GeometryInfo<2>::vertices_per_face]; -}; - -/** - * Workaround for a bug in egcs snapshot 1998/08/03. - */ -template <> -struct BoundaryHelper<3> { - /** - * Declare a data type for the derived - * classes. - */ - typedef const Point<3> *PointArray[GeometryInfo<3>::vertices_per_face]; -}; /** @@ -65,56 +21,97 @@ struct BoundaryHelper<3> { * following code (here in two dimensions): * \begin{verbatim} * ... - * const Point<2> *neighbors[2] = {&neighbor1, &neighbor2}; - * Point<2> new_vertex = boundary.in_between (neighbors); + * Point<2> new_vertex = boundary.get_new_point_on_line (line); * ... * \end{verbatim} - * #neighbor1# and #neighbor2# are the two vertices bounding the old - * line on the boundary, which is to be subdivided. #boundary# is an - * object of type #Boundary#. + * #line# denotes the line at the boundary that shall be refined + * and for which we seek the common point of the two child lines. * * In 3D, a new vertex may be placed on the middle of a line or on - * the middle of a side. In the both cases, an array with four points - * has to be passed to #in_between#; in the latter case the two end - * points of the line have to be given consecutively twice, as - * elements 0 and 1, and 2 and 3, respectively. + * the middle of a side. Respectively, the library calls + * \begin{verbatim} + * ... + * Point<3> new_line_vertices[4] + * = { boundary.get_new_point_on_line (face->line(0)), + * boundary.get_new_point_on_line (face->line(1)), + * boundary.get_new_point_on_line (face->line(2)), + * boundary.get_new_point_on_line (face->line(3)) }; + * ... + * \end{verbatim} + * to get the four midpoints of the lines bounding the quad at the + * boundary, and after that + * \begin{verbatim} + * ... + * Point<3> new_quad_vertex = boundary.get_new_point_on_quad (face); + * ... + * \end{verbatim} + * to get the midpoint of the face. It is guaranteed that this order + * (first lines, then faces) holds, so you can use information from + * the children of the four lines of a face, since these already exist + * at the time the midpoint of the face is to be computed. * + * Since iterators are passed to the functions, you may use information + * about boundary indicators and the like, as well as all other information + * provided by these objects. + * * There are specialisations, #StraightBoundary#, which places * the new point right into the middle of the given points, and * #HyperBallBoundary# creating a hyperball with given radius * around a given center point. * - * @author Wolfgang Bangerth, 1998 + * @author Wolfgang Bangerth, 1999 */ template class Boundary : public Subscriptor { public: - /** - * Typedef an array of the needed number - * of old points to compute the new - * middle point of a face. This does not - * make much sense in 1D, so we set the - * array size to a dummy value. - */ - typedef typename BoundaryHelper::PointArray PointArray; -// this is the way it should be, but egcs throws an internal compiler error -// on this, so we invented the above workaround -// typedef const Point* PointArray[((dim==1) ? -// 1 : -// GeometryInfo::vertices_per_face)]; - /** * Destructor. Does nothing here, but * needs to be declared to make it * virtual. */ virtual ~Boundary (); - + + /** + * Return the point which shall become + * the new middle vertex of the two + * children of a regular line. In 2D, + * this line is a line at the boundary, + * while in 3d, it is bounding a face + * at the boundary (the lines therefore + * is also on the boundary). + */ + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const = 0; + + /** + * Return the point which shall become + * the common point of the four children + * of a quad at the boundary in three + * or more spatial dimensions. This + * function therefore is only useful in + * at least three dimensions and should + * not be called for lower dimensions. + * + * This function is called after the + * four lines bounding the given #quad# + * are refined, so you may want to use + * the information provided by + * #quad->line(i)->child(j)#, #i=0...3#, + * #j=0,1#. + * + * Because in 2D, this function is not + * needed, it is not made pure virtual, + * to avoid the need to overload it. + * The default implementation throws + * an error in any case, however. + */ + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; + /** - * This function calculates the position - * of the new vertex. + * Exception. */ - virtual Point in_between (const PointArray &neighbors) const = 0; + DeclException0 (ExcPureVirtualFunctionCalled); }; @@ -135,10 +132,20 @@ template class StraightBoundary : public Boundary { public: /** - * This function calculates the position - * of the new vertex. + * Refer to the general documentation of + * this class and the documentation of the + * base class. */ - virtual Point in_between (const PointArray &neighbors) const; + 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; }; @@ -167,10 +174,20 @@ class HyperBallBoundary : public StraightBoundary { center(p), radius(radius) {}; /** - * This function calculates the position - * of the new vertex. + * 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 in_between (const PointArray &neighbors) const; + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; private: diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index a7e0e8761e..bf65e97833 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -4,6 +4,7 @@ #include #include +#include #include #include #include diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 4b6e4b3d2d..252765c8fe 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -4494,13 +4494,9 @@ void Triangulation<2>::execute_refinement () { // vertex? Point<2> new_point; if (cell->at_boundary(nb)) - { - // boundary vertex - const Point<2> *neighbors[2] = - {&vertices[new_vertices[2*nb]], - &vertices[new_vertices[(2*nb+2)%8]]}; - new_point = boundary->in_between (neighbors); - } else { + // boundary vertex + new_point = boundary->get_new_point_on_line (cell->line(nb)); + else { // vertex between two // normal cells new_point = vertices[new_vertices[2*nb]]; @@ -4945,17 +4941,9 @@ void Triangulation<3>::execute_refinement () { ExcTooFewVerticesAllocated()); vertices_used[next_unused_vertex] = true; - if (line->at_boundary()) - { - const Boundary::PointArray surrounding_points - = { &line->vertex(0), - &line->vertex(0), - &line->vertex(1), - &line->vertex(1) - }; - vertices[next_unused_vertex] - = boundary->in_between (surrounding_points); - } + if (line->at_boundary()) + vertices[next_unused_vertex] + = boundary->get_new_point_on_line (line); else vertices[next_unused_vertex] = (line->vertex(0) + line->vertex(1)) / 2; @@ -5036,16 +5024,8 @@ void Triangulation<3>::execute_refinement () { vertices_used[next_unused_vertex] = true; if (quad->at_boundary()) - { - const Boundary::PointArray surrounding_points - = { &quad->vertex(0), - &quad->vertex(1), - &quad->vertex(2), - &quad->vertex(3) - }; - vertices[next_unused_vertex] - = boundary->in_between (surrounding_points); - } + vertices[next_unused_vertex] + = boundary->get_new_point_on_quad (quad); else vertices[next_unused_vertex] = (quad->vertex(0) + quad->vertex(1) + diff --git a/deal.II/deal.II/source/grid/tria_boundary.cc b/deal.II/deal.II/source/grid/tria_boundary.cc index 7d21d33eb5..0bc85aa88b 100644 --- a/deal.II/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/deal.II/source/grid/tria_boundary.cc @@ -3,6 +3,9 @@ #include +#include +#include +#include #include @@ -12,22 +15,74 @@ Boundary::~Boundary () {}; + + +template +Point +Boundary::get_new_point_on_quad (const typename Triangulation::quad_iterator &) const +{ + Assert (false, ExcPureVirtualFunctionCalled()); + return Point(); +}; + + + +template +Point +StraightBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + return (line->vertex(0) + line->vertex(1)) / 2; +}; + + + +#if deal_II_dimension < 3 + +template +Point +StraightBoundary::get_new_point_on_quad (const typename Triangulation::quad_iterator &) const +{ + Assert (false, ExcPureVirtualFunctionCalled()); + return Point(); +}; + + +#else + + +template +Point +StraightBoundary::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + return (quad->vertex(0) + quad->vertex(1) + + quad->vertex(1) + quad->vertex(2)) / 2; +}; + +#endif + + + template -Point StraightBoundary::in_between (const PointArray &neighbors) const +Point +HyperBallBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const { - Point p; - for (int i=0; i<(1<<(dim-1)); ++i) - p += *neighbors[i]; - p /= (1<<(dim-1))*1.0; - return p; + 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::in_between (const PointArray &neighbors) const +Point +HyperBallBoundary::get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { - Point middle = StraightBoundary::in_between(neighbors); + Point middle = StraightBoundary::get_new_point_on_quad (quad); middle -= center; // project to boundary -- 2.39.5