]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Doc reformatting, update and add smoothing parameter to triangulation refinement.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 10 May 1998 00:20:25 +0000 (00:20 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 10 May 1998 00:20:25 +0000 (00:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@272 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/tria.h

index baf094f207815c6f8e3598b341b698067ababd3a..cd904d6a2295ac5eee1d55e31d72d579c147affb 100644 (file)
@@ -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 <int dim>
 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 <int dim>
 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 dim>
-       int Triangulation<dim>::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 dim>
-        int Triangulation<dim>::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 <int dim>
-        void Triangulation<dim>::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<dim>::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<dim>::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<dim>::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<dim>::read_*# functions. The appropriate function to be
-         called is #Triangulation<dim>::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<dim>::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 dim>
+ *     int Triangulation<dim>::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 dim>
+ *      int Triangulation<dim>::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 <int dim>
+ *      void Triangulation<dim>::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<dim>::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<dim>::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<dim>::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<dim>::read_*# functions. The appropriate function to be
+ *        called is #Triangulation<dim>::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<dim>::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 <int dim>
 class Triangulation : public TriaDimensionInfo<dim> {
   public:
@@ -1140,7 +1256,7 @@ class Triangulation : public TriaDimensionInfo<dim> {
                                      *  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<dim> {
                                      *  (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<dim> {
                                      *  Pointer to a boundary object.
                                      */
     const Boundary<dim>             *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<dim>;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.