From: Wolfgang Bangerth Date: Sun, 10 May 1998 00:20:25 +0000 (+0000) Subject: Doc reformatting, update and add smoothing parameter to triangulation refinement. X-Git-Tag: v8.0.0~22982 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9953b022e839810917750670f621c9828a299c77;p=dealii.git Doc reformatting, update and add smoothing parameter to triangulation refinement. git-svn-id: https://svn.dealii.org/trunk@272 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index baf094f207..cd904d6a22 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -38,36 +38,36 @@ class ostream; /** - General template for information belonging to one level of a multilevel - hierarchy of a triangulation. This template is only declared to allow - specializations for different dimensions. - - @see TriangulationLevel<1> - @see TriangulationLevel<2> - */ + * General template for information belonging to one level of a multilevel + * hierarchy of a triangulation. This template is only declared to allow + * specializations for different dimensions. + * + * @see TriangulationLevel<1> + * @see TriangulationLevel<2> + */ template class TriangulationLevel; /** - Store all information which belongs to one level of the multilevel hierarchy. - - In #TriangulationLevel<0># all data is stored which is not - dependant on the dimension, e.g. a field to store the - refinement flag for the cells (what a cell actually is - is declared elsewhere), etc. Actually, it is only cell-based - data, like neighborship info or refinement flags. There is another - field, which may fit in here, namely the material data (for cells) - or the boundary indicators (for faces), but since we need for a line - or quad either boundary information or material data, we store them - with the lines and quads rather than with the common data. We may, - however lose some memory in three dimensions, when we need the - material data for cell, boundary data for the quads, but nothing - for the lines. Since we only store one byte per line, quad or hex, - this is a minor loss and we can live with that. - - @memo Information belonging to one level of the multilevel hierarchy. - */ + * Store all information which belongs to one level of the multilevel hierarchy. + * + * In #TriangulationLevel<0># all data is stored which is not + * dependant on the dimension, e.g. a field to store the + * refinement flag for the cells (what a cell actually is + * is declared elsewhere), etc. Actually, it is only cell-based + * data, like neighborship info or refinement flags. There is another + * field, which may fit in here, namely the material data (for cells) + * or the boundary indicators (for faces), but since we need for a line + * or quad either boundary information or material data, we store them + * with the lines and quads rather than with the common data. We may, + * however lose some memory in three dimensions, when we need the + * material data for cell, boundary data for the quads, but nothing + * for the lines. Since we only store one byte per line, quad or hex, + * this is a minor loss and we can live with that. + * + * @memo Information belonging to one level of the multilevel hierarchy. + */ class TriangulationLevel<0> { public: /** @@ -185,15 +185,15 @@ class TriangulationLevel<0> { /** - Store all information which belongs to one level of the multilevel hierarchy. - - In one dimension, this is a list of the lines associated with this level, - as well as a list with the indices of the children of these lines. - The #TriangulationsLevel# objects of higher dimensions are derived from - this one. - - @memo Information belonging to one level of the multilevel hierarchy. - */ + * Store all information which belongs to one level of the multilevel hierarchy. + * + * In one dimension, this is a list of the lines associated with this level, + * as well as a list with the indices of the children of these lines. + * The #TriangulationsLevel# objects of higher dimensions are derived from + * this one. + * + * @memo Information belonging to one level of the multilevel hierarchy. + */ class TriangulationLevel<1> : public TriangulationLevel<0> { private: @@ -322,17 +322,17 @@ class TriangulationLevel<1> : public TriangulationLevel<0> { /** - Store all information which belongs to one level of the multilevel hierarchy. - - In 2D this is a vector of the lines and one of the - quads on this levels, as well as a the two associated vectors holding - information about the children of these lines and quads. - - The vector of lines and their children is derived from - #TriangulationLevel<1>#. - - @memo Information belonging to one level of the multilevel hierarchy. - */ + * Store all information which belongs to one level of the multilevel hierarchy. + * + * In 2D this is a vector of the lines and one of the + * quads on this levels, as well as a the two associated vectors holding + * information about the children of these lines and quads. + * + * The vector of lines and their children is derived from + * #TriangulationLevel<1>#. + * + * @memo Information belonging to one level of the multilevel hierarchy. + */ class TriangulationLevel<2> : public TriangulationLevel<1> { /** @@ -430,38 +430,39 @@ class TriangulationLevel<2> : public TriangulationLevel<1> /** - This class implements some types which differ between the dimensions. - Declare it to have a template parameter, but do not actually declare - other types than those explicitely instantiated. - */ + * This class implements some types which differ between the dimensions. + * Declare it to have a template parameter, but do not actually declare + * other types than those explicitely instantiated. + */ template class TriaDimensionInfo; + /** - This class implements some types which differ between the dimensions. - - A #line_iterator# is typdef'd to an iterator operating on the - #lines# member variable of a #Triangulation<1># object. An - #active_line_iterator# only operates on the active lines. - #raw_line_iterator# objects operate on all lines, used or not. - - Since we are in one dimension, the following identities are declared: - \begin{verbatim} - typedef raw_line_iterator raw_cell_iterator; - typedef line_iterator cell_iterator; - typedef active_line_iterator active_cell_iterator; - \end{verbatim} - - To enable the declaration of #begin_quad# and the like in - #Triangulation<1>#, the #quad_iterator#s are declared as - #void *#. Thus these types exist, but are useless and will - certainly make any involuntary use visible. - - The same applies for the #face_iterator# types, since lines - have no substructures apart from vertices, which are handled in - a different way, however. - */ + * This class implements some types which differ between the dimensions. + * + * A #line_iterator# is typdef'd to an iterator operating on the + * #lines# member variable of a #Triangulation<1># object. An + * #active_line_iterator# only operates on the active lines. + * #raw_line_iterator# objects operate on all lines, used or not. + * + * Since we are in one dimension, the following identities are declared: + * \begin{verbatim} + * typedef raw_line_iterator raw_cell_iterator; + * typedef line_iterator cell_iterator; + * typedef active_line_iterator active_cell_iterator; + * \end{verbatim} + * + * To enable the declaration of #begin_quad# and the like in + * #Triangulation<1>#, the #quad_iterator#s are declared as + * #void *#. Thus these types exist, but are useless and will + * certainly make any involuntary use visible. + * + * The same applies for the #face_iterator# types, since lines + * have no substructures apart from vertices, which are handled in + * a different way, however. + */ class TriaDimensionInfo<1> { public: typedef TriaRawIterator<1,CellAccessor<1> > raw_line_iterator; @@ -483,31 +484,31 @@ class TriaDimensionInfo<1> { /** - This class implements some types which differ between the dimensions. - - A #line_iterator# is typdef'd to an iterator operating on the - #lines# member variable of a #Triangulation<1># object. An - #active_line_iterator# only operates on the active lines. - #raw_line_iterator# objects operate on all lines, used or not. - Using #active_line_iterator#s may not be particularly useful since it - only operates on unrefined lines. However, also refined lines may bound - unrefined cells if the neighboring cell is refined once more than the - present one. - - Similarly, #quad_iterator#, #raw_quad_iterator# and - #active_quad_iterator# are declared. - - Since we are in two dimension, the following identities are declared: - \begin{verbatim} - typedef raw_quad_iterator raw_cell_iterator; - typedef quad_iterator cell_iterator; - typedef active_quad_iterator active_cell_iterator; - - typedef raw_line_iterator raw_face_iterator; - typedef line_iterator face_iterator; - typedef active_line_iterator active_face_iterator; - \end{verbatim} - */ + * This class implements some types which differ between the dimensions. + * + * A #line_iterator# is typdef'd to an iterator operating on the + * #lines# member variable of a #Triangulation<1># object. An + * #active_line_iterator# only operates on the active lines. + * #raw_line_iterator# objects operate on all lines, used or not. + * Using #active_line_iterator#s may not be particularly useful since it + * only operates on unrefined lines. However, also refined lines may bound + * unrefined cells if the neighboring cell is refined once more than the + * present one. + * + * Similarly, #quad_iterator#, #raw_quad_iterator# and + * #active_quad_iterator# are declared. + * + * Since we are in two dimension, the following identities are declared: + * \begin{verbatim} + * typedef raw_quad_iterator raw_cell_iterator; + * typedef quad_iterator cell_iterator; + * typedef active_quad_iterator active_cell_iterator; + * + * typedef raw_line_iterator raw_face_iterator; + * typedef line_iterator face_iterator; + * typedef active_line_iterator active_face_iterator; + * \end{verbatim} + */ class TriaDimensionInfo<2> { public: typedef TriaRawIterator<2,LineAccessor<2> > raw_line_iterator; @@ -536,583 +537,698 @@ class TriaDimensionInfo<2> { /** - #Triangulation#s denote a hierarchy of levels of elements which together - form a region in #dim# spatial dimensions. - - This class is written to be as independent of the dimension as possible - (thus the complex construction of the #TriangulationLevel# classes) to - allow code-sharing, to allow reducing the need to mirror changes in the code - for one dimenion to the code for other dimensions. Nonetheless, some of - the functions are dependent of the dimension and there only exist - specialized versions for distinct dimensions. - - - \subsection{Structure and iterators} - - The actual data structure of a #Triangulation# object is rather complex - and quite inconvenient if one attempted to operate on it directly, since - data is spread over quite a lot of arrays and other places. However, - there are ways powerful enough to work on these data structures - without knowing their exact relations. This is done through the - concept of iterators (see the STL documentation and \Ref{TriaRawIterator}). - In order to make things as easy and dimension independent as possible, - use of class local typedefs is made, see below. - - In the base class #TriaDimensionInfo#, a #Cell# is typedef'd to be whatever - is reasonable for a cell in the respective dimension, i.e. a #Line# in - one dimension, a #Quad# in two dimensions, and so on. - - The #Triangulation# class provides iterator which enable looping over all - lines, cells, - etc without knowing the exact representation used to describe them. Their - names are typedef's in the #TriaDimensionInfo# base class (thus making them - local types to this class) and are as follows: - - #raw_line_iterator#: loop over all lines, used or not (declared for - all dimensions). - - #line_iterator#: loop over all used lines (declared for all dimensions). - - #active_line_iterator#: loop over all active lines (declared for all - dimensions). - - #raw_quad_iterator#: loop over all quads, used or not (declared only - for #dim>=2#). - - #quad_iterator#: loop over all quads (declared only for #dim#>=2). - - #active_quad_iterator#: loop over all active quads (declared only for - #dim#>=2). - - Additionaly, for #dim#==1, the following identities hold: - \begin{verbatim} - typedef raw_line_iterator raw_cell_iterator; - typedef line_iterator cell_iterator; - typedef active_line_iterator active_cell_iterator; - \end{verbatim} - while for #dim#==2 - \begin{verbatim} - typedef quad_line_iterator raw_cell_iterator; - typedef quad_iterator cell_iterator; - typedef active_quad_iterator active_cell_iterator; - - typedef raw_line_iterator raw_face_iterator; - typedef line_iterator face_iterator; - typedef active_line_iterator active_face_iterator; - \end{verbatim} - - By using the cell iterators, you can write code nearly independent of - the spatial dimension. The same applies for substructure iterators, - where a substructure is defined as a face of a cell. The face of a - cell is be a vertex in 1D and a line in 2D; however, vertices are - handled in a different way and therefore lines have no faces. - - The #Triangulation# class offers functions like #begin_active# which gives - you an iterator to the first active cell. There are quite a lot of functions - returning iterators. Take a look at the class doc to get an overview. - - Usage of these iterators works mostly like with the STL iterators. Some - examples taken from the #Triangulation# source code follow. - \begin{itemize} - \item {\it Counting the number of cells on a specific level} - \begin{verbatim} - template - int Triangulation::n_cells (const int level) const { - cell_iterator cell = begin (level), - endc = (level == levels.size()-1 ? - cell_iterator(end()) : - begin (level+1)); - int n=0; - for (; cell!=endc; ++cell) - ++n; - return n; - }; - \end{verbatim} - Another way which uses the STL #distance# function would be to write - \begin{verbatim} - template - int Triangulation::n_cells (const int level) const { - int n=0; - distance (begin(level), - (level == levels.size()-1 ? - cell_iterator(end()) : - begin (level+1)), - n); - return n; - }; - \end{verbatim} - Unfortunately, #g++# presently (version 2.7.2) fails to find the right - #distance# template instantiation; it seems we have to wait for future - #g++# versions :-( - - \item {\it Refining all cells of a triangulation} - \begin{verbatim} - template - void Triangulation::refine_global () { - active_cell_iterator cell = begin_active(), - endc = end(); - - for (; cell != endc; ++cell) - cell->set_refine_flag (); - execute_refinement (); - }; - \end{verbatim} - \end{itemize} - - - \subsection{Usage} - - Usage of a #Triangulation# is mainly done through the use of iterators. - An example probably shows best how to use it: - \begin{verbatim} - void main () { - Triangulation<2> tria; - - // read in a coarse grid file - - // we want to log the - // refinement history - ofstream history ("mesh.history"); - - // refine first cell - tria.begin_active()->set_refine_flag(); - tria.save_refine_flags (history); - tria.execute_refinement (); - - // refine first active cell - // on coarsest level - tria.begin_active()->set_refine_flag (); - tria.save_refine_flags (history); - tria.execute_refinement (); - - Triangulation<2>::active_cell_iterator cell; - for (int i=0; i<17; ++i) - { - // refine the presently - // second last cell 17 - // times - cell = tria.last_active(tria.n_levels()-1); - --cell; - cell->set_refine_flag (); - tria.save_refine_flags (history); - tria.execute_refinement (); - }; - // output the grid - ofstream out("grid.1"); - tria.print_gnuplot (out); - }; - \end{verbatim} - - - \subsection{Creating a triangulation} - - There are several possibilities to create a triangulation: - \begin{itemize} - \item Hypercube triangulations: a hypercube triangulation is a - domain which is the tensor product of an interval $[a,b]$ in - the given number of spatial dimensions. If you want to create such - a domain, which is a common test case for model problems, call - #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 - more information. The mentionned class uses the interface described - directly below to transfer the data into the triangulation. - - \item Explicitely creating a triangulation: you can create a triangulation - by providing a list of vertices and a list of cells. Each such cell - consists of a vector storing the indices of the vertices of this cell - in the vertex list. To see how this works, you can take a look at the - #DataIn::read_*# functions. The appropriate function to be - called is #Triangulation::create_triangulation (2)#. - - Creating the hierarchical information needed for this library from - 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 - should therefore not be too large, reading refined grids would be - inefficient. Apart from the performance aspect, refined grids do not - lend too well to multigrid algorithms, since solving on the coarsest - 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. - 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} - - The material id for each cell must be specified upon construction of - a triangulation. (There is a special section on material ids and - boundary indicators. See there for more information.) - The standard region functions (for hypercube, hyperball, - etc.) denote all cells the material id zero. You may change that afterwards, - but you should not use the material id 255. When reading a triangulation, - the material id must be specified in the input file (UCD format) or is - otherwise set to zero. When creating explicitely, the material id must - be given to the creation function. - - Regarding the boundary indicator for lines in two dimensions and quads - in three (subsummed by the word "faces"), all interior faces are denoted - the value 255. Trying to give an interior face another value results in - an error if in debug mode. Faces at the boundary of the domain are preset - with the boundary indicator zero, but you can give a list of faces with - different boundary indicators to the triangulation creation function. - The standard domain functions assume all faces to have boundary indicator - zero, which you may change manually afterwards. When reading from a file, - you have to give boundary indicators other than zero explicitely, e.g. in - UCD format by giving a list of lines with material id in the input file. - - Lines in two dimensions and quads in three dimensions inherit their - boundary indicator to their children upon refinement. You should therefore - make sure that if you have different boundary parts, the different parts - are separated by a vertex (in 2D) or a line (in 3D) such that each boundary - line or quad has a unique boundary indicator. - - Likewise, material data is inherited from mother to child cells. Place your - coarse level cells so, that the interface between cells is also the - interface between regions of different materials. - - - \subsection{Refinement of a triangulation} - - Refinement of a triangulation may be done through several ways. The most - low-level way is directly through iterators: let #i# be an iterator to - an active cell (i.e. the cell pointed to has no children), then the - function call #i->set_refine_flag()# marks the respective cell for - refinement. Marking non-active cells results in an error. - - After all the cells you wanted to mark for refinement, call the - #execute_refinement# function to actually perform the refinement. This - function itself first calls the #prepare_refinement# function to smooth - the resulting triangulation: since a face between to adjacent cells may - only be subdivided once (i.e. the levels of two adjacent cells may - differ by one at most; it is not possible to have a cell refined twice - while the neighboring one is not refined), some additional cells are - flagged for refinement to smooth the grid. This enlarges the number of - resulting cells but makes the grid more regular, thus leading to better - approximationonal properties and, above all, making the handling of data - structures and algorithms much much easier. - - Marking cells for refinement 'by hand' through iterators is one way to - produce a new grid, especially if you know what kind of grid you are - looking for, e.g. if you want to have a grid successively refined - towards the boundary or always at the center (see the example programs, - they do exactly these things). There are more advanced functions, - however, which are more suitable for automatic generation of hierarchical - grids in the context of a-posteriori error estimation and adaptive finite - elements. - - The central function to this is - #refine (const dVector &criterion, const double threshold)#: it takes a - vector of values, one per active cell, which denote the criterion according - to which the triangulation is to be refined. It marks all cells for which - the criterion is greater than the threshold being given as the second - argument. - - There are two variations of this function, which rely on #refine# by - computing the threshold from other information: - \begin{itemize} - \item #refine_fixed_number#: this function takes a vector as above and - a value between zero and one denoting the fraction of cells to be - refined. For this purpose, it sorts the criteria per cell and takes - the threshold to be the one belonging to the cell with the - #fraction times n_active_cells# highest criterion. For example, if - the fraction is $0.3$, the threshold is computed to a value such that - 30 per cent of cells have a criterion higher than the threshold and are - thus flagged for refinement. The flagging for refinement is done through - the central #refine# function. - - The sorting of criteria is not done actually, since we only need one - value, in the example above the criterion of the cell which is at - 30 per cent in the sorted list of cells. The order of cells with higher - and of those with lower criteria is irrelevant. Getting this value is - accomplished by the #nth_element# function of the #C++# standard - library, which takes only linear time in the number of elements, rather - than #N log N# for sorting all values. - - A typical value for the fraction of cells to be refined is 0.3. - However, for singular functions or error functionals, you may want to - chose a smaller value to avoid overrefinement in regions which do not - contribute much to the error. - - \item #refine_fixed_fraction#: this function computes the threshold such - that the number of cells getting flagged for refinement makes up for a - certain fraction of the total error. If this fraction is 50 per cent, - for example, the threshold is computed such that the cells with a - criterion greater than the threshold together account for half of the - total error. - - This strategy is more suited for singular functions and error - functionals, but may lead to very slow convergence of the grid - if only few cells are refined in each step. - - From the implementational point, this time we really need to - sort the array of criteria. However, it is not necessary to sort - the whole array, since for example if you chose the fraction at - 50 per cent of the total error, it is only necessary to sort at - most the 50 per cent of cells ranking topmost in the list of error - per cell. It is thus reasonable to use an algorithm like - #partial_sort# of the C++ standard library, which only sorts part - of the array and lets the rest unsorted. However, in many cases - much fewer than 50 per cent of the cells account for 50 per cent - of the error, so it may be possible to get away with sorting less - than 50 per cent of the cells. We therefore divide the whole lot - of 50 per cent of cells into, say, 5 parts, first sort for the - 10 per cent with highest error; look whether they together make up - for 50 per cent and if so thats ok, we can leave the rest unsorted; - if not, sort the next 10 per cent, and so on. The default is to - devide the maximum number of cells which may get refined (which - equals the fraction of the total error, as explained above) into - five parts, but this value may be given as a parameter to the - #refine_fixed_fraction# function. For highly singular error - functionals, it may be more efficient to chose a greater number - than five. Chosing a value which is too large should not lead to - a large performance drawback; chosing too small a value however - may lead to significantly higher computational costs for sorting - than necessary. - - Just like the other strategy described above, this function only - computes the threshold value and then passes over to #refine#. - - A typical value for the fraction of the total error is 0.5. - \end{itemize} - - For a more thorough discussion of advantages and disadvantages of the - different strategies for refinement, see the paper of R. Becker and - R. Rannacher titled "A Feed-Back Approach to Error Control in Finite - Element Methods: Basic Analysis and Examples". - - It is assumed that the criterion is a value in a certain norm over each - element, such that the square of the total error is the sum over the - squares of the criteria on the cells. - - - \subsection{Material and boundary information} - - Each line, quad, etc stores one byte of information denoting the material - a cell is made of (used in the computation of stiffness matrices) and to - which part of the boundary it belongs. Obviously, the material id is what - is needed for a cell, while for all structures with a dimension less than - the dimension of the domain (i.e. lines in 2D, lines and quads in 3D), the - boundary information is what is needed. Since either material or boundary - information is needed, but never both, only one field is used to store this - data, namely the #TriangulationLevel<1>::LinesData.material_id# and - #TriangulationLevel<2>::QuadsData.material_id# vectors. The data can be - read and written using line, quad and cell iterators. - - Material and boundary indicators are stored as one byte (an - #unsigned char#). They can therefore be in the range zero to 255, but - only zero to 254 is allowed. The value 255 is reserved to denote - interior lines (in 2D) and interior lines and quads (in 3D), which need - not have a boundary or material indicator. However, using this value, it - is possible to say whether a line in 2D is interior or not, which would - otherwise be impossible because the hierarchical structure of a - triangulation stores neighborship information and the like only with - cells. Finding out whether a line is an interior one would then only be - possible by looking at the cell it belongs to. There would be no way to - loop over all lines and for example do a contour integral, since there - would be no way to find out which of the lines we loop over actually are - on the contour. - - Since in one dimension, no substructures of lower dimension exist to - cells (of course apart from vertices, but these are handled - in another way than the structures and substructures with dimension one - or greater), there is no way to denote boundary indicators to boundary - vertices (the endpoints). This is not a big thing, however, since you - will normally not want to do a loop over all vertices, but rather work - on the cells, which in turn have a possibility to find out whether they - are at one of the endpoints. Only the handling of boundary values gets - a bit more complicated, but this seems to be the price to be paid for - the different handling of vertices from lines and quads. - - - \subsection{History of a triangulation} - - It is possible to reconstruct a grid from its refinement history, which - can be stored and loaded through the #save_refine_flags# and - #load_refine_flags# functions. Normally, the code will look like this: - \begin{verbatim} - // open output file - ofstream history("mesh.history"); - // do 10 refinement steps - for (int step=0; step<10; ++step) { - ...; - // flag cells according to some criterion - ...; - tria.save_refine_flags (history); - tria.execute_refinement (); - }; - \end{verbatim} - - If you want to re-create the grid from the stored information, you write: - \begin{verbatim} - // open input file - ifstream history("mesh.history"); - // do 10 refinement steps - for (int step=0; step<10; ++step) { - tria.load_refine_flags (history); - tria.execute_refinement (); - }; - \end{verbatim} - - You may write other information to the output file between different sets - of refinement information, as long as you read it upon re-creation of the - grid. - - - \subsection{User flags} - - A triangulation offers one bit per line, quad, etc for user data. - This field can be - accessed as all other data using iterators. Normally, this user flag is - used if an algorithm walks over all cells and needs information whether - another cell, e.g. a neighbor, has already been processed. It can also - be used to flag the lines subject to constraints in 2D, as for example the - functions in the #DoFHandler# classes do. - - There are two functions, #save_user_flags# and #load_user_flags# which - write and read these flags to and from a stream. Unlike - #save_refine_flags# and #load_refine_flags#, these two functions store - and read the flags of all used lines, quads, etc, not only of the - active ones (well, activity is a concept which really only applies to - cells, not for example to lines in 2D, so the abovementioned generalisation - to {\it all} lines, quads, etc seems plausible). - - If you want to store more specific user flags, you can use the functions - #save_user_flags_line# and #load_user_flags_line# and the generalizations - for quads, etc. - - - It is convention to clear the user flags using the - #Triangulation<>::clear_user_flags()# function before usage, since it is - often necessary to use the flags in more than one function consecutively and - is then error prone to dedicate one of these to clear the flags. - - It is recommended that a functions using the flags states so in its - documentation. - - - \subsection{Boundary approximation} - - You can specify a boundary function: if a new vertex is created on a - side or face at the boundary, this function is used to compute where - it will be placed. See \Ref{Boundary} for the details. Usage with - the #Triangulation# object is then like this (let #Ball# be a class - derived from #Boundary<2>#): - \begin{verbatim} - void main () { - Triangulation<2> tria; - // set the boundary function - Ball ball; - tria.set_boundary (&ball); - - // read some coarse grid - - - Triangulation<2>::active_cell_iterator cell, endc; - for (int i=0; i<8; ++i) - { - cell = tria.begin_active(); - endc = tria.end(); - - // refine all - // boundary cells - for (; cell!=endc; ++cell) - if (cell->at_boundary()) - cell->set_refine_flag(); - - tria.execute_refinement(); - }; - }; - \end{verbatim} - - You should take note of one caveat: if you have concave boundaries, you - must make sure that a new boundary vertex does not lie to much inside the - to be refined cell. The reason is that the center vertex is placed at the - point which is the arithmetic mean of the eight surrounding vertices. - Therefore if your new boundary vertex is too near the center of the old - quadrilateral or hexahedron, the distance to the midpoint vertex will become - too small, thus generating distorted cells. Remedy: you have to take care - of such situations when defining the coarse grid. - - - \subsection{Implementational conventions for two spatial dimensions} - - There is a convention about the direction of the bounding lines of quads in - 2D. The direction of a line is the direction of point 0 towards point 1. We - define, that allowed cells contain of lines of which the direction is - as follows: - - 2 - 3--->---2 - | | - 3^ ^1 - | | - 0--->---1 - 0 - The number of the vertices and lines is also indicated. This orientation of - lines has to be checked/generated upon construction of a grid and is - preserved upon refinement. - - Further we define, that child lines have the same direction as their parent, - i.e. that #subline(0).vertex(0)==line.vertex(0)# and - #subline(1).vertex(1)==line.vertex(1)#. This also implies, that the - first subline (#subline(0)#) is the one at vertex(0) of the old line. - - Similarly we define, that the four children of a quad are adjacent to the - vertex with the same number of the old quad. - - - \subsection{Warning} - - It seems impossible to preserve #const#ness of a triangulation through - iterator usage. Thus, if you declare pointers to a #const# triangulation - object, you should be well aware that you might involuntarily alter the - data stored in the triangulation. - - @memo Implementation of a multilevel triangulation of a domain - @see TriaRawIterator - @author Wolfgang Bangerth, 1998 - */ + * #Triangulation#s denote a hierarchy of levels of elements which together + * form a region in #dim# spatial dimensions. + * + * This class is written to be as independent of the dimension as possible + * (thus the complex construction of the #TriangulationLevel# classes) to + * allow code-sharing, to allow reducing the need to mirror changes in the code + * for one dimenion to the code for other dimensions. Nonetheless, some of + * the functions are dependent of the dimension and there only exist + * specialized versions for distinct dimensions. + * + * + * \subsection{Structure and iterators} + * + * The actual data structure of a #Triangulation# object is rather complex + * and quite inconvenient if one attempted to operate on it directly, since + * data is spread over quite a lot of arrays and other places. However, + * there are ways powerful enough to work on these data structures + * without knowing their exact relations. This is done through the + * concept of iterators (see the STL documentation and \Ref{TriaRawIterator}). + * In order to make things as easy and dimension independent as possible, + * use of class local typedefs is made, see below. + * + * In the base class #TriaDimensionInfo#, a #Cell# is typedef'd to be whatever + * is reasonable for a cell in the respective dimension, i.e. a #Line# in + * one dimension, a #Quad# in two dimensions, and so on. + * + * The #Triangulation# class provides iterator which enable looping over all + * lines, cells, + * etc without knowing the exact representation used to describe them. Their + * names are typedef's in the #TriaDimensionInfo# base class (thus making them + * local types to this class) and are as follows: + * + * #raw_line_iterator#: loop over all lines, used or not (declared for + * all dimensions). + * + * #line_iterator#: loop over all used lines (declared for all dimensions). + * + * #active_line_iterator#: loop over all active lines (declared for all + * dimensions). + * + * #raw_quad_iterator#: loop over all quads, used or not (declared only + * for #dim>=2#). + * + * #quad_iterator#: loop over all quads (declared only for #dim#>=2). + * + * #active_quad_iterator#: loop over all active quads (declared only for + * #dim#>=2). + * + * Additionaly, for #dim#==1, the following identities hold: + * \begin{verbatim} + * typedef raw_line_iterator raw_cell_iterator; + * typedef line_iterator cell_iterator; + * typedef active_line_iterator active_cell_iterator; + * \end{verbatim} + * while for #dim#==2 + * \begin{verbatim} + * typedef quad_line_iterator raw_cell_iterator; + * typedef quad_iterator cell_iterator; + * typedef active_quad_iterator active_cell_iterator; + * + * typedef raw_line_iterator raw_face_iterator; + * typedef line_iterator face_iterator; + * typedef active_line_iterator active_face_iterator; + * \end{verbatim} + * + * By using the cell iterators, you can write code nearly independent of + * the spatial dimension. The same applies for substructure iterators, + * where a substructure is defined as a face of a cell. The face of a + * cell is be a vertex in 1D and a line in 2D; however, vertices are + * handled in a different way and therefore lines have no faces. + * + * The #Triangulation# class offers functions like #begin_active# which gives + * you an iterator to the first active cell. There are quite a lot of functions + * returning iterators. Take a look at the class doc to get an overview. + * + * Usage of these iterators works mostly like with the STL iterators. Some + * examples taken from the #Triangulation# source code follow. + * \begin{itemize} + * \item {\it Counting the number of cells on a specific level} + * \begin{verbatim} + * template + * int Triangulation::n_cells (const int level) const { + * cell_iterator cell = begin (level), + * endc = (level == levels.size()-1 ? + * cell_iterator(end()) : + * begin (level+1)); + * int n=0; + * for (; cell!=endc; ++cell) + * ++n; + * return n; + * }; + * \end{verbatim} + * Another way which uses the STL #distance# function would be to write + * \begin{verbatim} + * template + * int Triangulation::n_cells (const int level) const { + * int n=0; + * distance (begin(level), + * (level == levels.size()-1 ? + * cell_iterator(end()) : + * begin (level+1)), + * n); + * return n; + * }; + * \end{verbatim} + * Unfortunately, #g++# presently (version 2.7.2) fails to find the right + * #distance# template instantiation; it seems we have to wait for future + * #g++# versions :-( + * + * \item {\it Refining all cells of a triangulation} + * \begin{verbatim} + * template + * void Triangulation::refine_global () { + * active_cell_iterator cell = begin_active(), + * endc = end(); + * + * for (; cell != endc; ++cell) + * cell->set_refine_flag (); + * execute_refinement (); + * }; + * \end{verbatim} + * \end{itemize} + * + * + * \subsection{Usage} + * + * Usage of a #Triangulation# is mainly done through the use of iterators. + * An example probably shows best how to use it: + * \begin{verbatim} + * void main () { + * Triangulation<2> tria; + * + * // read in a coarse grid file + * + * // we want to log the + * // refinement history + * ofstream history ("mesh.history"); + * + * // refine first cell + * tria.begin_active()->set_refine_flag(); + * tria.save_refine_flags (history); + * tria.execute_refinement (); + * + * // refine first active cell + * // on coarsest level + * tria.begin_active()->set_refine_flag (); + * tria.save_refine_flags (history); + * tria.execute_refinement (); + * + * Triangulation<2>::active_cell_iterator cell; + * for (int i=0; i<17; ++i) + * { + * // refine the presently + * // second last cell 17 + * // times + * cell = tria.last_active(tria.n_levels()-1); + * --cell; + * cell->set_refine_flag (); + * tria.save_refine_flags (history); + * tria.execute_refinement (); + * }; + * // output the grid + * ofstream out("grid.1"); + * tria.print_gnuplot (out); + * }; + * \end{verbatim} + * + * + * \subsection{Creating a triangulation} + * + * There are several possibilities to create a triangulation: + * \begin{itemize} + * \item Hypercube triangulations: a hypercube triangulation is a + * domain which is the tensor product of an interval $[a,b]$ in + * the given number of spatial dimensions. If you want to create such + * a domain, which is a common test case for model problems, call + * #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 + * more information. The mentioned class uses the interface described + * directly below to transfer the data into the triangulation. + * + * \item Explicitely creating a triangulation: you can create a triangulation + * by providing a list of vertices and a list of cells. Each such cell + * consists of a vector storing the indices of the vertices of this cell + * in the vertex list. To see how this works, you can take a look at the + * #DataIn::read_*# functions. The appropriate function to be + * called is #Triangulation::create_triangulation (2)#. + * + * Creating the hierarchical information needed for this library from + * 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 + * should therefore not be too large, reading refined grids would be + * inefficient. Apart from the performance aspect, refined grids do not + * lend too well to multigrid algorithms, since solving on the coarsest + * 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. + * 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} + * + * The material id for each cell must be specified upon construction of + * a triangulation. (There is a special section on material ids and + * boundary indicators. See there for more information.) + * The standard region functions (for hypercube, hyperball, + * etc.) denote all cells the material id zero. You may change that afterwards, + * but you should not use the material id 255. When reading a triangulation, + * the material id must be specified in the input file (UCD format) or is + * otherwise set to zero. When creating explicitely, the material id must + * be given to the creation function. + * + * Regarding the boundary indicator for lines in two dimensions and quads + * in three (subsummed by the word "faces"), all interior faces are denoted + * the value 255. Trying to give an interior face another value results in + * an error if in debug mode. Faces at the boundary of the domain are preset + * with the boundary indicator zero, but you can give a list of faces with + * different boundary indicators to the triangulation creation function. + * The standard domain functions assume all faces to have boundary indicator + * zero, which you may change manually afterwards. When reading from a file, + * you have to give boundary indicators other than zero explicitely, e.g. in + * UCD format by giving a list of lines with material id in the input file. + * + * Lines in two dimensions and quads in three dimensions inherit their + * boundary indicator to their children upon refinement. You should therefore + * make sure that if you have different boundary parts, the different parts + * are separated by a vertex (in 2D) or a line (in 3D) such that each boundary + * line or quad has a unique boundary indicator. + * + * Likewise, material data is inherited from mother to child cells. Place your + * coarse level cells so, that the interface between cells is also the + * interface between regions of different materials. + * + * + * \subsection{Refinement of a triangulation} + * + * Refinement of a triangulation may be done through several ways. The most + * low-level way is directly through iterators: let #i# be an iterator to + * an active cell (i.e. the cell pointed to has no children), then the + * function call #i->set_refine_flag()# marks the respective cell for + * refinement. Marking non-active cells results in an error. + * + * After all the cells you wanted to mark for refinement, call the + * #execute_refinement# function to actually perform the refinement. This + * function itself first calls the #prepare_refinement# function to regularise + * the resulting triangulation: since a face between two adjacent cells may + * only be subdivided once (i.e. the levels of two adjacent cells may + * differ by one at most; it is not possible to have a cell refined twice + * while the neighboring one is not refined), some additional cells are + * flagged for refinement to smooth the grid. This enlarges the number of + * resulting cells but makes the grid more regular, thus leading to better + * approximationonal properties and, above all, making the handling of data + * structures and algorithms much much easier. To be honest, this is mostly + * an algorithmic step than one needed by the finite element method. + * + * However, some degradation of approximation properties has been observed + * for grids which were refined more than once across a face, as descibed + * above, so there is also a practical justification for the above. + * It can also be shown, that such degradation occurs if the + * triangulation contains vertices which are member of cells with levels + * differing by more than one level. On such example is the following: + * \begin{verbatim} + * | | | | + * x-----x-----x--x--x-- + * | | | | | + * | | x--x--x + * | | | | | + * x-----x-----x--x--x-- + * | | | + * | | | + * | | | + * | x-----x-- + * | | | + * | | | + * | | | + * x-----------x-----x-- + * \end{verbatim} + * It seems that in two space dimensions, the maximu jump in levels between + * cells sharing a common vertex is two (as in the example above), this is + * not true if more than four cells meet at a vertex. It is not uncommon + * that a coarse (initial) mesh contains vertices at which six or even eight + * cells meet, when small features of the domain have to be resolved even on + * the coarsest mesh. In that case, the maximum difference in levels is + * three or four, respectively. The problem gets even worse in three space + * dimensions. + * + * Looking at an interpolation of the second derivative of the finite + * element solution (asuming bilinear finite elements), one sees that the + * numerical solution is almost totally wrong, compared with the true second + * derivative. Indeed, on regular meshes, there exist sharp estimations that + * the $H^2$-error is only $O(1)$, so we should not be suprised; however, the + * numerical solution may show a value for the second derivative which may + * be a factor of ten away from the true value. These problems are located + * on the small cell adjacent to the center vertex, where cells of + * non-subsequent levels meet, as well as on the upper and right neighbor + * of this cell (but with a less degree of deviation from the true value). + * + * Due to the approximational problems described above, the + * #Triangulation# constructor takes an argument specifying whether a + * smoothing step shall be performed on the grid each time #execute_refinement# + * is called. The default is that such a step not be done, since this results + * in additional cells being produced, which may not be necessary in all + * cases. If switched on, calling #execute_refinement# results in + * flagging additional cells for refinement to avoid + * vertices as the ones mentioned. The algorithms for both regularisation + * and smoothing of triangulations are described below in the section on + * technical issues. The reason why this parameter must be given to the + * constructor rather than to #execute_refinement# is that it would result + * in algorithmic problems if you called #execute_refinement# once without + * and once with smoothing, since then in some refinement steps would need + * to be refined twice. + * + * Marking cells for refinement 'by hand' through iterators is one way to + * produce a new grid, especially if you know what kind of grid you are + * looking for, e.g. if you want to have a grid successively refined + * towards the boundary or always at the center (see the example programs, + * they do exactly these things). There are more advanced functions, + * however, which are more suitable for automatic generation of hierarchical + * grids in the context of a-posteriori error estimation and adaptive finite + * elements. + * + * The central function to this is + * #refine (const dVector &criterion, const double threshold)#: it takes a + * vector of values, one per active cell, which denote the criterion according + * to which the triangulation is to be refined. It marks all cells for which + * the criterion is greater than the threshold being given as the second + * argument. + * + * There are two variations of this function, which rely on #refine# by + * computing the threshold from other information: + * \begin{itemize} + * \item #refine_fixed_number#: this function takes a vector as above and + * a value between zero and one denoting the fraction of cells to be + * refined. For this purpose, it sorts the criteria per cell and takes + * the threshold to be the one belonging to the cell with the + * #fraction times n_active_cells# highest criterion. For example, if + * the fraction is $0.3$, the threshold is computed to a value such that + * 30 per cent of cells have a criterion higher than the threshold and are + * thus flagged for refinement. The flagging for refinement is done through + * the central #refine# function. + * + * The sorting of criteria is not done actually, since we only need one + * value, in the example above the criterion of the cell which is at + * 30 per cent in the sorted list of cells. The order of cells with higher + * and of those with lower criteria is irrelevant. Getting this value is + * accomplished by the #nth_element# function of the #C++# standard + * library, which takes only linear time in the number of elements, rather + * than #N log N# for sorting all values. + * + * A typical value for the fraction of cells to be refined is 0.3. + * However, for singular functions or error functionals, you may want to + * chose a smaller value to avoid overrefinement in regions which do not + * contribute much to the error. + * + * \item #refine_fixed_fraction#: this function computes the threshold such + * that the number of cells getting flagged for refinement makes up for a + * certain fraction of the total error. If this fraction is 50 per cent, + * for example, the threshold is computed such that the cells with a + * criterion greater than the threshold together account for half of the + * total error. + * + * This strategy is more suited for singular functions and error + * functionals, but may lead to very slow convergence of the grid + * if only few cells are refined in each step. + * + * From the implementational point, this time we really need to + * sort the array of criteria. However, it is not necessary to sort + * the whole array, since for example if you chose the fraction at + * 50 per cent of the total error, it is only necessary to sort at + * most the 50 per cent of cells ranking topmost in the list of error + * per cell. It is thus reasonable to use an algorithm like + * #partial_sort# of the C++ standard library, which only sorts part + * of the array and lets the rest unsorted. However, in many cases + * much fewer than 50 per cent of the cells account for 50 per cent + * of the error, so it may be possible to get away with sorting less + * than 50 per cent of the cells. We therefore divide the whole lot + * of 50 per cent of cells into, say, 5 parts, first sort for the + * 10 per cent with highest error; look whether they together make up + * for 50 per cent and if so thats ok, we can leave the rest unsorted; + * if not, sort the next 10 per cent, and so on. The default is to + * devide the maximum number of cells which may get refined (which + * equals the fraction of the total error, as explained above) into + * five parts, but this value may be given as a parameter to the + * #refine_fixed_fraction# function. For highly singular error + * functionals, it may be more efficient to chose a greater number + * than five. Chosing a value which is too large should not lead to + * a large performance drawback; chosing too small a value however + * may lead to significantly higher computational costs for sorting + * than necessary. + * + * Just like the other strategy described above, this function only + * computes the threshold value and then passes over to #refine#. + * + * A typical value for the fraction of the total error is 0.5. + * \end{itemize} + * + * For a more thorough discussion of advantages and disadvantages of the + * different strategies for refinement, see the paper of R. Becker and + * R. Rannacher titled "A Feed-Back Approach to Error Control in Finite + * Element Methods: Basic Analysis and Examples". + * + * It is assumed that the criterion is a value in a certain norm over each + * element, such that the square of the total error is the sum over the + * squares of the criteria on the cells. + * + * + * \subsection{Material and boundary information} + * + * Each line, quad, etc stores one byte of information denoting the material + * a cell is made of (used in the computation of stiffness matrices) and to + * which part of the boundary it belongs. Obviously, the material id is what + * is needed for a cell, while for all structures with a dimension less than + * the dimension of the domain (i.e. lines in 2D, lines and quads in 3D), the + * boundary information is what is needed. Since either material or boundary + * information is needed, but never both, only one field is used to store this + * data, namely the #TriangulationLevel<1>::LinesData.material_id# and + * #TriangulationLevel<2>::QuadsData.material_id# vectors. The data can be + * read and written using line, quad and cell iterators. + * + * Material and boundary indicators are stored as one byte (an + * #unsigned char#). They can therefore be in the range zero to 255, but + * only zero to 254 is allowed. The value 255 is reserved to denote + * interior lines (in 2D) and interior lines and quads (in 3D), which need + * not have a boundary or material indicator. However, using this value, it + * is possible to say whether a line in 2D is interior or not, which would + * otherwise be impossible because the hierarchical structure of a + * triangulation stores neighborship information and the like only with + * cells. Finding out whether a line is an interior one would then only be + * possible by looking at the cell it belongs to. There would be no way to + * loop over all lines and for example do a contour integral, since there + * would be no way to find out which of the lines we loop over actually are + * on the contour. + * + * Since in one dimension, no substructures of lower dimension exist to + * cells (of course apart from vertices, but these are handled + * in another way than the structures and substructures with dimension one + * or greater), there is no way to denote boundary indicators to boundary + * vertices (the endpoints). This is not a big thing, however, since you + * will normally not want to do a loop over all vertices, but rather work + * on the cells, which in turn have a possibility to find out whether they + * are at one of the endpoints. Only the handling of boundary values gets + * a bit more complicated, but this seems to be the price to be paid for + * the different handling of vertices from lines and quads. + * + * + * \subsection{History of a triangulation} + * + * It is possible to reconstruct a grid from its refinement history, which + * can be stored and loaded through the #save_refine_flags# and + * #load_refine_flags# functions. Normally, the code will look like this: + * \begin{verbatim} + * // open output file + * ofstream history("mesh.history"); + * // do 10 refinement steps + * for (int step=0; step<10; ++step) { + * ...; + * // flag cells according to some criterion + * ...; + * tria.save_refine_flags (history); + * tria.execute_refinement (); + * }; + * \end{verbatim} + * + * If you want to re-create the grid from the stored information, you write: + * \begin{verbatim} + * // open input file + * ifstream history("mesh.history"); + * // do 10 refinement steps + * for (int step=0; step<10; ++step) { + * tria.load_refine_flags (history); + * tria.execute_refinement (); + * }; + * \end{verbatim} + * + * You may write other information to the output file between different sets + * of refinement information, as long as you read it upon re-creation of the + * grid. + * + * + * \subsection{User flags} + * + * A triangulation offers one bit per line, quad, etc for user data. + * This field can be + * accessed as all other data using iterators. Normally, this user flag is + * used if an algorithm walks over all cells and needs information whether + * another cell, e.g. a neighbor, has already been processed. It can also + * be used to flag the lines subject to constraints in 2D, as for example the + * functions in the #DoFHandler# classes do. + * + * There are two functions, #save_user_flags# and #load_user_flags# which + * write and read these flags to and from a stream. Unlike + * #save_refine_flags# and #load_refine_flags#, these two functions store + * and read the flags of all used lines, quads, etc, not only of the + * active ones (well, activity is a concept which really only applies to + * cells, not for example to lines in 2D, so the abovementioned generalisation + * to {\it all} lines, quads, etc seems plausible). + * + * If you want to store more specific user flags, you can use the functions + * #save_user_flags_line# and #load_user_flags_line# and the generalizations + * for quads, etc. + * + * + * It is convention to clear the user flags using the + * #Triangulation<>::clear_user_flags()# function before usage, since it is + * often necessary to use the flags in more than one function consecutively and + * is then error prone to dedicate one of these to clear the flags. + * + * It is recommended that a functions using the flags states so in its + * documentation. + * + * + * \subsection{Boundary approximation} + * + * You can specify a boundary function: if a new vertex is created on a + * side or face at the boundary, this function is used to compute where + * it will be placed. See \Ref{Boundary} for the details. Usage with + * the #Triangulation# object is then like this (let #Ball# be a class + * derived from #Boundary<2>#): + * \begin{verbatim} + * void main () { + * Triangulation<2> tria; + * // set the boundary function + * Ball ball; + * tria.set_boundary (&ball); + * + * // read some coarse grid + * + * + * Triangulation<2>::active_cell_iterator cell, endc; + * for (int i=0; i<8; ++i) + * { + * cell = tria.begin_active(); + * endc = tria.end(); + * + * // refine all + * // boundary cells + * for (; cell!=endc; ++cell) + * if (cell->at_boundary()) + * cell->set_refine_flag(); + * + * tria.execute_refinement(); + * }; + * }; + * \end{verbatim} + * + * You should take note of one caveat: if you have concave boundaries, you + * must make sure that a new boundary vertex does not lie to much inside the + * to be refined cell. The reason is that the center vertex is placed at the + * point which is the arithmetic mean of the eight surrounding vertices. + * Therefore if your new boundary vertex is too near the center of the old + * quadrilateral or hexahedron, the distance to the midpoint vertex will become + * too small, thus generating distorted cells. Remedy: you have to take care + * of such situations when defining the coarse grid. + * + * + * \subsection{Technical details} + * + * \subsubsection{Algorithms for mesh regularisation and smoothing upon refinement} + * + * We chose an inductive point of view: since upon creation of the + * triangulation all cells are on the same level, all regularity assumptions + * regarding the maximum difference in level of cells sharing a common face, + * edge or vertex hold. Since we use the regularisation and smoothing in + * each step of the mesh history, when coming to the point of refining it + * further the assumptions also hold. + * + * The regularisation and smoothing is done in the #prepare_refinement# + * function, which is called by #execute_refinement# at the very beginning. + * It decides which additional cells to flag for refinement by looking at the + * old grid and the refinement flags for each cell. + * + * \begin{itemize} + * \item {\it Regularisation:} The algorithm walks over all cells checking + * whether the present cell is flagged for refinement and a neighbor of the + * present cell is refined once less than the present one. If so, flag the + * neighbor for refinement. Because of the induction above, there may be no + * neighbor with level two less than the present one. + * + * The neighbor thus flagged for refinement may induce more cells which need + * to be refined. However, such cells which need additional refinement always + * are on one level lower than the present one, so we can get away with only + * one sweep over all cells if we do the loop in the reverse way, starting + * with those on the highest level. This way, we may flag additional cells + * on lower levels, but if these induce more refinement needed, this is + * performed later on when we visit them in out backward running loop. + * + * \item {\it Smoothing:} First a list is set up which stores for each vertex + * the highest level one of the adjacent cells belongs to. Now, since we did + * smoothing in the previous refinement steps also, each cell may only have + * vertices with levels at most one greater than the level of the present + * cell. + * + * However, if we store the level plus one for cells marked for refinement, + * we may end up with cells which have vertices of level two greater than + * the cells level. We need to refine this cell also, and need thus also + * update the levels of its vertices. This itself may lead to cells needing + * refinement, but these are on lower levels, as above, which is why we + * may do all kinds of additional flagging in one loop only. + * \end{itemize} + * + * Regularisation and smoothing are a bit complementary in that we check + * whether we need to set additional refinement flags when being on a cell + * flagged for refinement (regularisation) or on a cell not flagged for + * refinement. This makes readable programming easier. + * + * All the described algorithms apply only for more than one space dimension, + * since for one dimension no restrictions apply. It may be necessary to + * apply some smoothing for multigrid algorithms, but this has to be decided + * upon later. + * + * + * \subsubsection{Implementational conventions for two spatial dimensions} + * + * There is a convention about the direction of the bounding lines of quads in + * 2D. The direction of a line is the direction of point 0 towards point 1. We + * define, that allowed cells contain of lines of which the direction is + * as follows: + * \begin{verbatim} + * 2 + * 3--->---2 + * | | + * 3^ ^1 + * | | + * 0--->---1 + * 0 + * \end{verbatim} + * The number of the vertices and lines is also indicated. This orientation of + * lines has to be checked/generated upon construction of a grid and is + * preserved upon refinement. + * + * Further we define, that child lines have the same direction as their parent, + * i.e. that #subline(0).vertex(0)==line.vertex(0)# and + * #subline(1).vertex(1)==line.vertex(1)#. This also implies, that the + * first subline (#subline(0)#) is the one at vertex(0) of the old line. + * + * Similarly we define, that the four children of a quad are adjacent to the + * vertex with the same number of the old quad. + * + * + * \subsection{Warning} + * + * It seems impossible to preserve #const#ness of a triangulation through + * iterator usage. Thus, if you declare pointers to a #const# triangulation + * object, you should be well aware that you might involuntarily alter the + * data stored in the triangulation. + * + * @memo Implementation of a multilevel triangulation of a domain + * @see TriaRawIterator + * @author Wolfgang Bangerth, 1998 + */ template class Triangulation : public TriaDimensionInfo { public: @@ -1140,7 +1256,7 @@ class Triangulation : public TriaDimensionInfo { * the first level of the hierarchy. * Do not create any cells. */ - Triangulation (); + Triangulation (const bool smooth_grid = false); /** * Delete the object and all levels of @@ -2075,6 +2191,12 @@ class Triangulation : public TriaDimensionInfo { * (make sure that no two cells are * adjacent with a refinement level * differing with more than one), etc. + * It performs some mesh smoothing if + * the according flag was given to the + * constructor of this class. + * See the general + * doc of this class for more information. + * * This function is mostly dimension * independent. */ @@ -2107,7 +2229,14 @@ class Triangulation : public TriaDimensionInfo { * Pointer to a boundary object. */ const Boundary *boundary; - + + /** + * Do some smoothing in the process + * of refining the triangulation. See + * the general doc of this class for + * more information about this. + */ + bool smooth_grid; // Friendship includes local classes. friend class TriaAccessor;