From 5b3ea890e99c5d3b6d673026021ae6e9fc259f04 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 16 Mar 1998 08:34:45 +0000 Subject: [PATCH] Add grid input functionality and new boundary objects. git-svn-id: https://svn.dealii.org/trunk@67 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/tria.h | 107 ++++++++++++++++++- deal.II/deal.II/include/grid/tria_boundary.h | 59 +++++++++- 2 files changed, 159 insertions(+), 7 deletions(-) diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 684787ce82..72dc910872 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -8,7 +8,7 @@ #include #include #include -#include + //#ifdef __GNUC__ @@ -19,6 +19,7 @@ //forward declaration needed +template class Boundary; template class TriaAccessor; template class LineAccessor; @@ -346,7 +347,7 @@ class TriangulationLevel<2> : public TriangulationLevel<1> /** * Assert that enough space is allocated - * to accomodate #new_quads# new lines. + * to accomodate #new_quads# new quads. */ void reserve_space (const unsigned int new_quads); @@ -638,6 +639,26 @@ class TriaDimensionInfo<2> { #Triangulation::create_hypercube (a,b)#, which produces a hypercube domain triangulated with exactly one element. You can get tensor product meshes by successive refinement of this cell. + + \item Other standard regions: you can get the generalized L-shape domain + using the #Triangulation::create_L_region (a,b)# function, which + is the hypercube with the interval $[a,b]$ without the hypercube + made out of the interval $[(a+b)/2,b]$. Let, for example, be $a=-1$ + and $b=1$, then the hpyer-L in two dimensions is the region + $[-1,1]^2 - [0,1]^2$. To create a hyper-L in one dimension results in + an error. + + You get the circle or ball (or generalized: hyperball) around origin + #p# and with radius #r# by calling + #Triangulation::create_hyper_ball (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 + an error. + + Do not forget to attach a suitable + boundary approximation object if you want the triangulation to be refined + at the outer boundaries. \item Reading in a triangulation: By using an object of the \Ref{#DataIn#} class, you can read in fairly general triangulations. See there for @@ -652,7 +673,7 @@ class TriaDimensionInfo<2> { called is #Triangulation::create_triangulation (2)#. Creating the hierarchical information needed for this library from - cells storing only vertex information can be a quite complex task. + cells storing only vertex information can be quite a complex task. For example in 2d, we have to create lines between vertices (but only once, though there are two cells which link these two vertices) and we have to create neighborship information. Grids being read in @@ -662,7 +683,25 @@ class TriaDimensionInfo<2> { level is expensive. It is wiser in any case to read in a grid as coarse as possible and then do the needed refinement steps. - It is your duty to guarantee that cells have the correct orientation. + It is your duty to guarantee that cells have the correct orientation. + To guarantee this, in the input vector keeping the cell list, the + vertex indices for each cell have to be in a defined order. In one + dimension, the first vertex index must refer to that vertex with the + lower coordinate value. In two dimensions, the four vertices must be + given in an order representing a counterclockwise sense. This + condition is not easy to verify and no full attempt to do so is made. + If you violate this condition, you may end up with matrix entries + having the wrong sign (clockwise vertex numbering, which results in + a negative area element) of with wrong matrix elements (twisted + quadrilaterals, i.e. two vertices interchanged; this results in + a wrong area element). + + There are more subtle conditions which must be imposed upon the + vertex numbering within cells. See the documentation for the + \Ref{DataIn} class for more details on this. They do not only + hold for the data read from an UCD or any other input file, but + also for the data passed to the + #Triangulation::create_triangulation (2)# function. \end{itemize} @@ -852,6 +891,11 @@ class Triangulation : public TriaDimensionInfo { * be empty upon calling this function and * the cell list should be useful (connected * domain, etc.). + * + * The numbering of vertices within the + * #cells# array is subject to some + * constraints; see the general class + * documentation for this. */ void create_triangulation (const vector > &vertices, const vector > &cells); @@ -873,8 +917,34 @@ class Triangulation : public TriaDimensionInfo { void create_hypercube (const double left = 0., const double right= 1.); + /** + * Initialize the triangulation with a + * hyper-L consisting of exactly #2^dim-1# + * cells. See the general documentation for a + * description of the L-region. The limits + * default to minus unity and unity. + * + * The triangulation needs to be void + * upon calling this function. + */ + void create_hyper_L (const double left = -1., + const double right= 1.); - + /** + * Initialize the triangulation with a + * hyperball, i.e. a circle or a ball. + * See the general documentation for a + * more concise description. The center of + * the hyperball default to the origin, + * the radius defaults to unity. + * + * The triangulation needs to be void + * upon calling this function. + */ + void create_hyper_ball (const Point center = Point(), + const double radius = 1.); + + /** * @name Mesh refinement */ @@ -1435,6 +1505,33 @@ class Triangulation : public TriaDimensionInfo { /** * Exception */ + DeclException1 (ExcGridHasInvalidCell, + int, + << "Something went wrong when making cell " << arg1 + << ". Read the docs and the source code " + << "for more information."); + /** + * Exception + */ + DeclException0 (ExcGridHasInvalidVertices); + /** + * Exception + */ + DeclException1 (ExcInternalErrorOnCell, + int, + << "Something went wrong upon construction of cell " + << arg1); + /** + * Exception + */ + DeclException3 (ExcInvalidVertexIndex, + int, int, int, + << "Error while creating cell " << arg1 + << ": the vertex index " << arg2 << " must be between 0 and " + << arg3 << "."); + /** + * Exception + */ DeclException0 (ExcInternalError); //@} protected: diff --git a/deal.II/deal.II/include/grid/tria_boundary.h b/deal.II/deal.II/include/grid/tria_boundary.h index 4bf0245e90..4df6ee48e2 100644 --- a/deal.II/deal.II/include/grid/tria_boundary.h +++ b/deal.II/deal.II/include/grid/tria_boundary.h @@ -27,8 +27,10 @@ points of the line have to be given consecutively twice, as elements 0 and 1, and 2 and 3, respectively. - There is a specialisation, #StraightBoundary#, which places - the new point right into the middle of the given points. + 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. */ template class Boundary { @@ -74,6 +76,59 @@ 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 have to be given upon construction of + an object of this type. + + 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, const double radius) : + center(p), radius(radius) {}; + + /** + * This function calculates the position + * of the new vertex. + */ + virtual Point in_between (const PointArray &neighbors) const { + Point middle = StraightBoundary::in_between(neighbors); + + middle -= center; + // project to boundary + middle *= radius / sqrt(middle.square()); + + middle += center; + return middle; + }; + + + 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 */ #endif -- 2.39.5