]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move the grid refinement functions from the triangulation class to a class of its...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 28 Jul 2000 14:27:29 +0000 (14:27 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 28 Jul 2000 14:27:29 +0000 (14:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@3219 0785d39b-7218-0410-832d-ea1e28bc413d

13 files changed:
deal.II/deal.II/include/grid/grid_generator.h
deal.II/deal.II/include/grid/grid_refinement.h [new file with mode: 0644]
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/source/grid/grid_refinement.cc [new file with mode: 0644]
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/numerics/time_dependent.cc
deal.II/examples/step-6/step-6.cc
deal.II/examples/step-7/step-7.cc
deal.II/examples/step-8/step-8.cc
deal.II/examples/step-9/step-9.cc
tests/big-tests/error-estimation/error-estimation.cc
tests/big-tests/multigrid/multigrid.cc
tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc

index 25d2f04827a05580fa2c92afd10863f893ea28f0..f37d336a49f158b0fffe53d82d99203120b23da4 100644 (file)
@@ -29,12 +29,12 @@ template <int dim> class Triangulation;
  *       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
- *       @p{GridGenerator::hyper_cube (tria, a,b)}, which produces a
+ *       @ref{GridGenerator}@p{::hyper_cube (tria, 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 Generalized L-shape domain:
- *      using the @p{GridGenerator::hyper_L (tria, a,b)} function produces
+ *      using the @ref{GridGenerator}@p{::hyper_L (tria, a,b)} function produces
  *      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
@@ -44,7 +44,7 @@ template <int dim> class Triangulation;
  *    @item Hyper balls:
  *      You get the circle or ball (or generalized: hyperball) around origin
  *      @p{p} and with radius @p{r} by calling
- *      @p{GridGenerator::hyper_ball (tria, p, r)}. The circle is triangulated
+ *      @ref{GridGenerator}@p{::hyper_ball (tria, p, r)}. The circle is triangulated
  *      by five cells, the ball by seven cells. The diameter of the center cell is
  *      chosen so that the aspect ratio of the boundary cells after one refinement
  *      is minimized in some way. To create a hyperball in one dimension results in
@@ -57,7 +57,7 @@ template <int dim> class Triangulation;
  *    @item Hyper shell: A hyper shell is the region between two hyper
  *      sphere with the same origin. Therefore, it is a ring in two
  *      spatial dimensions. To triangulation it, call the function
- *      @pGridGenerator::hyper_shell (tria, origin, inner_radius, outer_radius, N)},
+ *      @ref{GridGenerator}@p{::hyper_shell (tria, origin, inner_radius, outer_radius, N)},
  *      where the center of the spheres as well as
  *      the inner and outer radius of the two spheres are given as
  *      shown.
@@ -72,7 +72,7 @@ template <int dim> class Triangulation;
  *      used in the radial direction.
  *
  *      You need to attach a boundary object to the triangulation. A
- *      suitable boundary class is provided as @p{HyperSphereBoundary}
+ *      suitable boundary class is provided as @ref{HyperSphereBoundary}
  *      in the library.
  *
  * @item Slit domain: The slit domain is a variant of the hyper cube
@@ -90,18 +90,23 @@ class GridGenerator
 {
   public:
                                     /**
-                                     * Initialize the given triangulation with a
-                                     * hypercube (line in 1D, square in 2D, etc)
-                                     * consisting of exactly one cell. The
-                                     * hypercube volume is the tensor product
-                                     * of the intervall $[left,right]$ in the
-                                     * present number of dimensions, where
-                                     * the limits are given as arguments. They
-                                     * default to zero and unity, then producing
-                                     * the unit hypercube.
+                                     * Initialize the given
+                                     * triangulation with a hypercube
+                                     * (line in 1D, square in 2D,
+                                     * etc) consisting of exactly one
+                                     * cell. The hypercube volume is
+                                     * the tensor product of the
+                                     * intervall $[left,right]$ in
+                                     * the present number of
+                                     * dimensions, where the limits
+                                     * are given as arguments. They
+                                     * default to zero and unity,
+                                     * then producing the unit
+                                     * hypercube.
                                      *
-                                     * The triangulation needs to be void
-                                     * upon calling this function.
+                                     * The triangulation needs to be
+                                     * void upon calling this
+                                     * function.
                                      */
     template <int dim>
     static void hyper_cube (Triangulation<dim> &tria,
@@ -136,15 +141,19 @@ class GridGenerator
                                     bool                colorize = false);
 
                                     /**
-                                     * Initialize the given triangulation with a
-                                     * hyperball, i.e. a circle or a ball.
-                                     * See the general documentation for a
-                                     * more concise description. The center of
-                                     * the hyperball default to the origin,
-                                     * the radius defaults to unity.
+                                     * Initialize the given
+                                     * triangulation with a
+                                     * hyperball, i.e. a circle or a
+                                     * ball.  See the general
+                                     * documentation for a more
+                                     * concise description. The
+                                     * center of the hyperball
+                                     * default to the origin, the
+                                     * radius defaults to unity.
                                      *
-                                     * The triangulation needs to be void
-                                     * upon calling this function.
+                                     * The triangulation needs to be
+                                     * void upon calling this
+                                     * function.
                                      */    
     template <int dim>
     static void hyper_ball (Triangulation<dim> &tria,
@@ -152,14 +161,18 @@ class GridGenerator
                            const double        radius = 1.);
 
                                     /**
-                                     * Initialize the given triangulation with a
-                                     * hyper-L consisting of exactly @p{2^dim-1}
-                                     * cells. See the general documentation for a
-                                     * description of the L-region. The limits
-                                     * default to minus unity and unity.
+                                     * Initialize the given
+                                     * triangulation with a hyper-L
+                                     * consisting of exactly
+                                     * @p{2^dim-1} cells. See the
+                                     * general documentation for a
+                                     * description of the
+                                     * L-region. The limits default
+                                     * to minus unity and unity.
                                      *
-                                     * The triangulation needs to be void
-                                     * upon calling this function.
+                                     * The triangulation needs to be
+                                     * void upon calling this
+                                     * function.
                                      */
     template <int dim>
     static void hyper_L (Triangulation<dim> &tria,
@@ -170,7 +183,8 @@ class GridGenerator
                                      * Initialize the given
                                      * Triangulation with a hypercube
                                      * with a slit. The slit goes
-                                     * from @p{(x=0,y=-1)} to @p{(0,0)} in 2d.
+                                     * from @p{(x=0,y=-1)} to
+                                     * @p{(0,0)} in 2d.
                                      *
                                      * The triangulation needs to be void
                                      * upon calling this function.
@@ -196,8 +210,9 @@ class GridGenerator
                                      * the resulting elements have
                                      * the least aspect ratio.
                                      *
-                                     * The triangulation needs to be void
-                                     * upon calling this function.
+                                     * The triangulation needs to be
+                                     * void upon calling this
+                                     * function.
                                      */
     template <int dim>
     static void hyper_shell (Triangulation<dim> &tria,
diff --git a/deal.II/deal.II/include/grid/grid_refinement.h b/deal.II/deal.II/include/grid/grid_refinement.h
new file mode 100644 (file)
index 0000000..b93e3bf
--- /dev/null
@@ -0,0 +1,266 @@
+//----------------------------  grid_refinement.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  grid_refinement.h  ---------------------------
+#ifndef __deal2__grid_refinement_h
+#define __deal2__grid_refinement_h
+
+
+// forward declarations
+template <int dim> class Triangulation;
+template <class T> class Vector;
+
+
+
+/**
+ *   This class provides several function that flag certain cells for
+ *   coarsening or refinement based on a vector of ``error''
+ *   indicators and some selection algorithm.  The central function is
+ *   @p{refine (const Vector<float> &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. Analogously,
+ *   @p{coarsen (const Vector<float> &criterion, const double threshold)}
+ *   flags those cells for
+ *   coarsening for which the criterion is less than the threshold.
+ *
+ *   There are two variations of these functions, which rely on @p{refine} and
+ *   @p{coarsen} by computing the thresholds from other information:
+ *   @begin{itemize}
+ *   @item @p{refine_and_coarsen_fixed_number}: this function takes a vector as
+ *     above and two values between zero and one denoting the fractions of cells to
+ *     be refined and coarsened. For this purpose, it sorts the criteria per cell
+ *     and takes the threshold to be the one belonging to the cell with the
+ *     @p{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 @p{refine} function. For coarsening, the same holds.
+ *
+ *     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 @p{nth_element} function of the @p{C++} standard
+ *     library, which takes only linear time in the number of elements, rather
+ *     than @p{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 singular error functionals, you may
+ *     want to chose a smaller value to avoid overrefinement in regions which
+ *     do not contribute much to the error.
+ *
+ *   @item @p{refine_and_coarsen_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. The definition of the fraction is a bit counterintuitive, since
+ *     the total error is the sum over all cells of the local contribution
+ *     squared. We define that the fraction $\alpha$ be such that those
+ *     elements with the greatest error are refined for which the condition
+ *     $\sum \eta_K^2 \le \alpha\eta^2$ holds. Note that $\alpha$ is not
+ *     squared. The sum runs over the mentioned
+ *     cells, $\eta_K$ are the local error indicators and $\eta$ is the global
+ *     indicator with $\eta^2 = \sum \eta_K^2$, with here the sum running over
+ *     all cells.
+ *
+ *     For the bottom fraction the same holds: the threshold for coarsening is
+ *     computed such that the cells with criterion less than the threshold
+ *     together make up for the fraction of the total error specified.
+ *
+ *     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 point of view of implementation, this time we really need to
+ *     sort the array of criteria.
+ *     Just like the other strategy described above, this function only
+ *     computes the threshold values and then passes over to @p{refine} and
+ *     @p{coarsen}.
+ *
+ *     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. The criteria shall be positive.
+ *
+ *   You can suppress coarsening or refining by giving zero as the fraction
+ *   for one of the operations.
+ *
+ * @author Wolfgang Bangerth, 1998, 2000
+ */
+class GridRefinement
+{
+  public:
+                                    /**
+                                     * Refine the triangulation
+                                     * according to the given
+                                     * criteria. The criterion is a
+                                     * @p{double} value for each cell
+                                     * which determines which cells
+                                     * are to be refined by
+                                     * comparison with the threshold:
+                                     * if the value for a cell is
+                                     * larger than the threshold, the
+                                     * cell is flagged for
+                                     * refinement. It is your duty to
+                                     * guarantee that the threshold
+                                     * value is in a resonable
+                                     * range. Please note that the
+                                     * @p{criteria} array may contain
+                                     * negative values (sometimes,
+                                     * error estimators are evaluated
+                                     * in a way which produces
+                                     * positive and negative values),
+                                     * but the comparison with
+                                     * @p{threshold} is done only on
+                                     * the absolute values of the
+                                     * criteria.
+                                     *
+                                     * The cells are only flagged for
+                                     * refinement, they are not
+                                     * actually refined. To do so,
+                                     * you have to call the
+                                     * @p{execute_coarsening_and_refinement}
+                                     * function.
+                                     *
+                                     * There are more sophisticated
+                                     * strategies for mesh
+                                     * refinement; refer to the
+                                     * following functions and to the
+                                     * general doc for this class for
+                                     * more information.
+                                     *
+                                     * Note that this function takes
+                                     * a vector of @p{float}s, rather
+                                     * than the usual @p{double}s,
+                                     * since accuracy is not so much
+                                     * needed here and saving memory
+                                     * may be a good goal when using
+                                     * many cells.
+                                     */
+    template <int dim, typename number>
+    static void refine (Triangulation<dim>   &tria,
+                       const Vector<number> &criteria,
+                       const double         threshold);
+
+                                    /**
+                                     * Analogue to the @p{refine}
+                                     * function: flag all cells for
+                                     * coarsening for which the
+                                     * absolute value of the
+                                     * criterion is less than the
+                                     * given threshold.
+                                     *
+                                     * Note that this function takes
+                                     * a vector of @p{float}s, rather
+                                     * than the usual @p{double}s,
+                                     * since accuracy is not so much
+                                     * needed here and saving memory
+                                     * may be a good goal when using
+                                     * many cells.
+                                     */
+    template <int dim, typename number>
+    static void coarsen (Triangulation<dim>   &tria,
+                        const Vector<number> &criteria,
+                        const double         threshold);
+    
+                                    /**
+                                     * Refine the triangulation by
+                                     * refining a certain fraction
+                                     * @p{top_fraction_of_cells} with
+                                     * the highest error. Likewise
+                                     * coarsen the fraction
+                                     * @p{bottom_fraction_of_cells}
+                                     * with the least error. To
+                                     * actually perform the
+                                     * refinement, call
+                                     * @p{execute_coarsening_and_refinement}.
+                                     *
+                                     * @p{fraction_of_cells} shall be
+                                     * a value between zero and one.
+                                     *
+                                     * Refer to the general doc of
+                                     * this class for more
+                                     * information.
+                                     *
+                                     * Note that this function takes
+                                     * a vector of @p{float}s, rather
+                                     * than the usual @p{double}s,
+                                     * since accuracy is not so much
+                                     * needed here and saving memory
+                                     * may be a good goal when using
+                                     * many cells.
+                                     */
+    template <int dim, typename number>
+    static void refine_and_coarsen_fixed_number (Triangulation<dim>   &tria,
+                                                const Vector<number> &criteria,
+                                                const double         top_fraction_of_cells,
+                                                const double         bottom_fraction_of_cells);
+
+                                    /**
+                                     * Refine the triangulation by
+                                     * flagging those cells which
+                                     * make up a certain
+                                     * @p{top_fraction} of the total
+                                     * error.  Likewise, coarsen all
+                                     * cells which make up only
+                                     * @p{bottom_fraction}.  To
+                                     * actually perform the
+                                     * refinement, call
+                                     * @p{execute_coarsening_and_refinement}.
+                                     *
+                                     * @p{*_fraction} shall be a
+                                     * values between zero and one.
+                                     *
+                                     * Refer to the general doc of
+                                     * this class for more
+                                     * information.
+                                     *
+                                     * Note that this function takes
+                                     * a vector of @p{float}s, rather
+                                     * than the usual @p{double}s,
+                                     * since accuracy is not so much
+                                     * needed here and saving memory
+                                     * may be a good goal when using
+                                     * many cells.
+                                     */
+    template<int dim, typename number>
+    static void refine_and_coarsen_fixed_fraction (Triangulation<dim>   &tria,
+                                                  const Vector<number> &criteria,
+                                                  const double         top_fraction,
+                                                  const double         bottom_fraction);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidVectorSize,
+                   int, int,
+                   << "The given vector has " << arg1
+                   << " elements, but " << arg2 << " were expected.");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidParameterValue);
+};
+
+
+
+#endif //__deal2__grid_refinement_h
index f4d5422639a0f66d2341904555e00efe2c88f176..d9df9075a65708bd66c0d72f57ef089587210e33 100644 (file)
@@ -34,12 +34,13 @@ template <int dim> class MGDoFHandler;
 
 
 /**
- *  Structure which is passed to the @p{Triangulation<dim>::create_triangulation}
+ *  Structure which is passed to the @ref{Triangulation}@p{<dim>::create_triangulation}
  *  function. It contains all data needed to construct a cell, namely the
  *  indices of the vertices and the material indicator.
  */
 template <int dim>
-struct CellData {
+struct CellData
+{
 #if !((__GNUC__==2) && (__GNUC_MINOR__==95))
     int           vertices[GeometryInfo<dim>::vertices_per_cell];
 #else
@@ -50,7 +51,7 @@ struct CellData {
 
 
 /**
- *  Structure to be passed to the @p{Triangulation<dim>::create_triangulation}
+ *  Structure to be passed to the @ref{Triangulation}@p{<dim>::create_triangulation}
  *  function to describe boundary information.
  *
  *  This structure is the same for all dimensions, since we use an input
@@ -343,7 +344,7 @@ struct TriaNumberCache<1>
  * Cache class used to store the number of used and active elements
  * (lines or quads etc) within the levels of a triangulation. This
  * specialization stores the numbers of quads. Due to the inheritance
- * from the base class @p{TriaNumberCache<1>}, the numbers of lines
+ * from the base class @ref{TriaNumberCache<1>}, the numbers of lines
  * are also within this class.
  *
  * In the old days, whenever one wanted to access one of these
@@ -396,7 +397,7 @@ struct TriaNumberCache<2> : public TriaNumberCache<1>
  * Cache class used to store the number of used and active elements
  * (lines or quads etc) within the levels of a triangulation. This
  * specialization stores the numbers of hexes. Due to the inheritance
- * from the base class @p{TriaNumberCache<2>}, the numbers of lines
+ * from the base class @ref{TriaNumberCache<2>}, the numbers of lines
  * and quads are also within this class.
  *
  * In the old days, whenever one wanted to access one of these
@@ -453,7 +454,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *  form a region in @p{dim} spatial dimensions.
  *
  *  This class is written to be as independent of the dimension as possible
- *  (thus the complex construction of the @p{TriangulationLevel} classes) to
+ *  (thus the complex construction of the @ref{TriangulationLevel} classes) to
  *  allow code-sharing, to allow reducing the need to mirror changes in the code
  *  for one dimension to the code for other dimensions. Nonetheless, some of
  *  the functions are dependent of the dimension and there only exist
@@ -471,31 +472,33 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *  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 @p{TriaDimensionInfo}, a @p{Cell} is typedef'd to be whatever
+ *  In the base class @ref{TriaDimensionInfo}, a @p{Cell} is typedef'd to be whatever
  *  is reasonable for a cell in the respective dimension, i.e. a @p{Line} in
  *  one dimension, a @p{Quad} in two dimensions, and so on.
  *
  *  The @p{Triangulation} class provides iterator which enable looping over all
  *  lines, cells,
  *  etc without knowing the exact representation used to describe them. Their
- *  names are typedefs in the @p{TriaDimensionInfo} base class (thus making them
+ *  names are typedefs in the @ref{TriaDimensionInfo} base class (thus making them
  *  local types to this class) and are as follows:
  *
- *  @p{raw_line_iterator}: loop over all lines, used or not (declared for
+ *  @begin{itemize}
+ *  @item @p{raw_line_iterator}: loop over all lines, used or not (declared for
  *  all dimensions).
  *  
- *  @p{line_iterator}: loop over all used lines (declared for all dimensions).
+ *  @item @p{line_iterator}: loop over all used lines (declared for all dimensions).
  *
- *  @p{active_line_iterator}: loop over all active lines (declared for all
+ *  @item @p{active_line_iterator}: loop over all active lines (declared for all
  *  dimensions).
  *
- *  @p{raw_quad_iterator}: loop over all quads, used or not (declared only
+ *  @item @p{raw_quad_iterator}: loop over all quads, used or not (declared only
  *  for @p{dim>=2}).
  *  
- *  @p{quad_iterator}: loop over all quads (declared only for @p{dim}>=2).
+ *  @item @p{quad_iterator}: loop over all quads (declared only for @p{dim}>=2).
  *
- *  @p{active_quad_iterator}: loop over all active quads (declared only for
+ *  @item @p{active_quad_iterator}: loop over all active quads (declared only for
  *  @p{dim}>=2).
+ *  @end{itemize}
  *
  *  Additionaly, for @p{dim}==1, the following identities hold:
  *  @begin{verbatim}
@@ -619,11 +622,11 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *    @item The most common domains, such as hypercubes (i.e. lines, squares,
  *       cubes, etc), hyper-balls (circles, balls, ...) and some other, more
  *       weird domains such as the L-shape region and higher dimensional
- *       generalizations and others, are provided by the @p{GridGenerator}
+ *       generalizations and others, are provided by the @ref{GridGenerator}
  *       class which takes a triangulation and fills it by a division
  *       of the required domain.
  *   
- *     @item Reading in a triangulation: By using an object of the @p{GridIn}
+ *     @item Reading in a triangulation: By using an object of the @ref{GridIn}
  *        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.
@@ -632,7 +635,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *        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
- *        @p{GridIn<dim>::read_*} functions. The appropriate function to be
+ *        @ref{GridIn}@p{<dim>::read_*} functions. The appropriate function to be
  *        called is @p{Triangulation<dim>::create_triangulation (2)}.
  *
  *        Creating the hierarchical information needed for this library from
@@ -777,86 +780,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *   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
- *   @p{refine (const Vector<float> &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. Analogously,
- *   @p{coarsen (const Vector<float> &criterion, const double threshold)} flags those
- *   cells for coarsening for which the criterion is less than the threshold.
- *
- *   There are two variations of these functions, which rely on @p{refine} and
- *   coarsen by computing the thresholds from other information:
- *   @begin{itemize}
- *   @item @p{refine_and_coarsen_fixed_number}: this function takes a vector as
- *     above and two values between zero and one denoting the fractions of cells to
- *     be refined and coarsened. For this purpose, it sorts the criteria per cell
- *     and takes the threshold to be the one belonging to the cell with the
- *     @p{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 @p{refine} function. For coarsening, the same holds.
- *
- *     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 @p{nth_element} function of the @p{C++} standard
- *     library, which takes only linear time in the number of elements, rather
- *     than @p{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 singular error functionals, you may
- *     want to chose a smaller value to avoid overrefinement in regions which
- *     do not contribute much to the error.
- *
- *   @item @p{refine_and_coarsen_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. The definition of the fraction is a bit counterintuitive, since
- *     the total error is the sum over all cells of the local contribution
- *     squared. We define that the fraction $\alpha$ be such that those
- *     elements with the greatest error are refined for which the condition
- *     $\sum \eta_K^2 \le \alpha\eta^2$ holds. Note that $\alpha$ is not
- *     squared. The sum runs over the mentioned
- *     cells, $\eta_K$ are the local error indicators and $\eta$ is the global
- *     indicator with $\eta^2 = \sum \eta_K^2$, with here the sum running over
- *     all cells.
- *
- *     For the bottom fraction the same holds: the threshold for coarsening is
- *     computed such that the cells with criterion less than the threshold
- *     together make up for the fraction of the total error specified.
- *
- *     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 point of view of implementation, this time we really need to
- *     sort the array of criteria.
- *     Just like the other strategy described above, this function only
- *     computes the threshold values and then passes over to @p{refine} and
- *     @p{coarsen}.
- *
- *     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. The criteria shall be positive.
- *
- *   You can suppress coarsening or refining by giving zero as the fraction
- *   for one of the operations.
+ *   elements. These functions can be found in the @ref{GridRefinement} class.
  *
  *
  *   @sect3{Smoothing of a triangulation}
@@ -993,7 +917,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *     Ensures patch level 1. As result the triangulation consists of
  *     patches, i.e. of cells that are refined once. It follows that
  *     if at least one of the children of a cell is or will be refined
- *     than all children need to be refined. If the @p{path_level_1} flag
+ *     than all children need to be refined. If the @p{patch_level_1} flag
  *     is set, than the flags @p{eliminate_unrefined_islands},
  *     @p{eliminate_refined_inner_islands} and
  *     @p{eliminate_refined_boundary_islands} will be ignored as they will
@@ -1111,7 +1035,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *   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 @p{DoFHandler} classes do.
+ *   functions in the @ref{DoFHandler} classes do.
  *
  *   There are two functions, @p{save_user_flags} and @p{load_user_flags} which
  *   write and read these flags to and from a stream. Unlike
@@ -1164,7 +1088,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *   placed. The boundary indicator of the face will be used to
  *   determine the proper component. See @ref{Boundary} for the
  *   details. Usage with the @p{Triangulation} object is then like this
- *   (let @p{Ball} be a class derived from @p{Boundary<2>}):
+ *   (let @p{Ball} be a class derived from @ref{Boundary}@p{<2>}):
  * 
  *   @begin{verbatim}
  *     void main () {
@@ -1515,7 +1439,7 @@ struct TriaNumberCache<3> : public TriaNumberCache<2>
  *   @item Face 4: children 3, 2, 6, 7;
  *   @item Face 5: children 0, 4, 7, 3.
  *   @end{itemize}
- *   You can get these numbers using the @p{GeometryInfo<3>::child_cell_on_face}
+ *   You can get these numbers using the @ref{GeometryInfo<3>}@p{::child_cell_on_face}
  *   function. Each child is adjacent to the vertex with the same number.
  *
  *
@@ -1571,8 +1495,10 @@ class Triangulation : public TriaDimensionInfo<dim>,
 {
   private:
                                     /**
-                                     * Default boundary object. This declaration is used
-                                     * for the default argument in @p{set_boundary}.
+                                     * Default boundary object. This
+                                     * declaration is used for the
+                                     * default argument in
+                                     * @p{set_boundary}.
                                      */
     static const StraightBoundary<dim>& straight_boundary;
 
@@ -1583,7 +1509,7 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                      * for mesh smoothing
                                      * algorithms. The meaning of
                                      * these flags is documented in
-                                     * the @p{Triangulation} class.
+                                     * the @ref{Triangulation} class.
                                      */
     enum MeshSmoothing
     {
@@ -1661,7 +1587,7 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                      *
                                      * Note that this operation is only allowed
                                      * if no subscriptions to this object exist
-                                     * any more, such as @p{DoFHandler} objects
+                                     * any more, such as @ref{DoFHandler} objects
                                      * using it.
                                      */
     void clear ();
@@ -1699,7 +1625,7 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                      * copied and MUST persist until
                                      * the triangulation is
                                      * destroyed. Otherwise, the
-                                     * @p{Subscriptor} class will issue
+                                     * @ref{Subscriptor} class will issue
                                      * @p{ExcObjectInUse}.  This is
                                      * also true for triangulations
                                      * generated from this one by
@@ -1724,62 +1650,75 @@ class Triangulation : public TriaDimensionInfo<dim>,
                       const Boundary<dim> &boundary_object = straight_boundary);
 
                                     /**
-                                     * Return a constant reference to a boundary
-                                     * object used for this triangulation.
-                                     * Number is the same as in @p{set_boundary}
+                                     * Return a constant reference to
+                                     * a boundary object used for
+                                     * this triangulation.  Number is
+                                     * the same as in
+                                     * @p{set_boundary}
                                      */
     const Boundary<dim> & get_boundary (unsigned int number) const;
     
                                     /**
-                                     *  Copy a triangulation. This operation is
-                                     *  not cheap, so you should be careful
-                                     *  with using this. We do not implement
-                                     *  this function as a copy constructor,
-                                     *  since it makes it easier to maintain
-                                     *  collections of triangulations if you
-                                     *  can assign them values later on.
+                                     *  Copy a triangulation. This
+                                     *  operation is not cheap, so
+                                     *  you should be careful with
+                                     *  using this. We do not
+                                     *  implement this function as a
+                                     *  copy constructor, since it
+                                     *  makes it easier to maintain
+                                     *  collections of triangulations
+                                     *  if you can assign them values
+                                     *  later on.
                                      *
-                                     *  Keep in mind that this function also
-                                     *  copies the pointer to the boundary
-                                     *  descriptor previously set by the
-                                     *  @p{set_boundary} function. You must
-                                     *  therefore also guarantee that the
-                                     *  boundary objects has a lifetime at
-                                     *  least as long as the copied
-                                     *  triangulation.
+                                     *  Keep in mind that this
+                                     *  function also copies the
+                                     *  pointer to the boundary
+                                     *  descriptor previously set by
+                                     *  the @p{set_boundary}
+                                     *  function. You must therefore
+                                     *  also guarantee that the
+                                     *  boundary objects has a
+                                     *  lifetime at least as long as
+                                     *  the copied triangulation.
                                      *
-                                     *  This triangulation must be empty
-                                     *  beforehand.
+                                     *  This triangulation must be
+                                     *  empty beforehand.
                                      *
                                      *  The function is made
-                                     *  @p{virtual} since some derived
-                                     *  classes might want to disable
-                                     *  or extend the functionality
-                                     *  of this function.
+                                     *  @p{virtual} since some
+                                     *  derived classes might want to
+                                     *  disable or extend the
+                                     *  functionality of this
+                                     *  function.
                                      */
     virtual void copy_triangulation (const Triangulation<dim> &old_tria);
 
                                     /**
-                                     * Create a triangulation from a list
-                                     * of vertices and a list of cells, each of
-                                     * the latter being a list of @p{1<<dim}
-                                     * vertex indices. The triangulation must
-                                     * be empty upon calling this function and
-                                     * the cell list should be useful (connected
-                                     * domain, etc.).
+                                     * Create a triangulation from a
+                                     * list of vertices and a list of
+                                     * cells, each of the latter
+                                     * being a list of @p{1<<dim}
+                                     * vertex indices. The
+                                     * triangulation must be empty
+                                     * upon calling this function and
+                                     * the cell list should be useful
+                                     * (connected domain, etc.).
                                      *
-                                     * Material data for the cells is given
-                                     * within the @p{cells} array, while boundary
+                                     * Material data for the cells is
+                                     * given within the @p{cells}
+                                     * array, while boundary
                                      * information is given in the
                                      * @p{subcelldata} field.
                                      *
-                                     * The numbering of vertices within the
-                                     * @p{cells} array is subject to some
-                                     * constraints; see the general class
+                                     * The numbering of vertices
+                                     * within the @p{cells} array is
+                                     * subject to some constraints;
+                                     * see the general class
                                      * documentation for this.
                                      *
-                                     * This function is made @p{virtual} to allow
-                                     * derived classes to set up some data
+                                     * This function is made
+                                     * @p{virtual} to allow derived
+                                     * classes to set up some data
                                      * structures as well.
                                      */
     virtual void create_triangulation (const vector<Point<dim> >    &vertices,
@@ -1787,18 +1726,21 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                       const SubCellData            &subcelldata);
 
                                     /**
-                                     * Distort the grid by randomly moving
-                                     * around all the vertices of the grid.
-                                     * The direction of moving is random,
-                                     * while the length of the shift vector
-                                     * has a value of @p{factor} times the
-                                     * minimal length of the active lines
-                                     * adjacent to this vertex. Note that
-                                     * @p{factor} should obviously be well
-                                     * below @p{0.5}.
+                                     * Distort the grid by randomly
+                                     * moving around all the vertices
+                                     * of the grid.  The direction of
+                                     * moving is random, while the
+                                     * length of the shift vector has
+                                     * a value of @p{factor} times
+                                     * the minimal length of the
+                                     * active lines adjacent to this
+                                     * vertex. Note that @p{factor}
+                                     * should obviously be well below
+                                     * @p{0.5}.
                                      *
-                                     * If @p{keep_boundary} is set to @p{true}
-                                     * (which is the default), then boundary
+                                     * If @p{keep_boundary} is set to
+                                     * @p{true} (which is the
+                                     * default), then boundary
                                      * vertices are not moved.
                                      */
     void distort_random (const double factor,
@@ -1810,140 +1752,40 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                      */
                                     /*@{*/
                                     /**
-                                     *  Flag all active cells for refinement.
-                                     *  This will refine
-                                     *  all cells of all levels which are not
-                                     *  already refined (i.e. only cells are
-                                     *  refined which do not yet have
-                                     *  children). The cells are only flagged,
-                                     *  not refined, thus you have the chance
-                                     *  to save the refinement flags.
+                                     *  Flag all active cells for
+                                     *  refinement.  This will refine
+                                     *  all cells of all levels which
+                                     *  are not already refined
+                                     *  (i.e. only cells are refined
+                                     *  which do not yet have
+                                     *  children). The cells are only
+                                     *  flagged, not refined, thus
+                                     *  you have the chance to save
+                                     *  the refinement flags.
                                      */
     void set_all_refine_flags ();
 
                                     /**
-                                     *  Refine all cells @p{times} times, by
-                                     *  alternatingly calling @p{refine_global()}
-                                     *  and @p{execute_coarsening_and_refinement()}.
-                                     *  This function actually starts the
-                                     *  refinement process, so you have no way
-                                     *  to store the refinement flags unless
-                                     *  you overload the
+                                     *  Refine all cells @p{times}
+                                     *  times, by alternatingly
+                                     *  calling
+                                     *  @p{set_all_refine_flags()}
+                                     *  and
+                                     *  @p{execute_coarsening_and_refinement()}.
+                                     *  This function actually starts
+                                     *  the refinement process, so
+                                     *  you have no way to store the
+                                     *  refinement flags unless you
+                                     *  overload the
                                      *  @p{execute_coarsening_and_refinement}
                                      *  function.
                                      */
     void refine_global (const unsigned int times);
 
                                     /**
-                                     * Refine the triangulation according to
-                                     * the given criteria. The criterion is a
-                                     * @p{double} value for each cell which
-                                     * determines which cells are to be refined
-                                     * by comparison with the threshold: if the
-                                     * value for a cell is larger than the
-                                     * threshold, the cell is flagged for
-                                     * refinement. It is your duty to guarantee
-                                     * that the threshold value is in a
-                                     * resonable range. Please note that the
-                                     * @p{criteria} array may contain negative
-                                     * values (sometimes, error estimators
-                                     * are evaluated in a way which produces
-                                     * positive and negative values), but the
-                                     * comparison with @p{threshold} is done only
-                                     * on the absolute values of the criteria.
-                                     *
-                                     * The cells are only flagged for
-                                     * refinement, they are not actually
-                                     * refined. To do so, you have to call the
-                                     * @p{execute_coarsening_and_refinement} function.
-                                     *
-                                     * There are more sophisticated strategies
-                                     * for mesh refinement; refer to the
-                                     * following functions and to the general
-                                     * doc for this class for more information.
-                                     *
-                                     * Note that this function takes a vector
-                                     * of @p{float}s, rather than the usual
-                                     * @p{double}s, since accuracy is not so
-                                     * much needed here and saving memory may
-                                     * be a good goal when using many cells.
-                                     */
-    template <typename number>
-    void refine (const Vector<number> &criteria,
-                const double         threshold);
-
-                                    /**
-                                     * Analogue to the @p{refine} function:
-                                     * flag all cells for coarsening for
-                                     * which the absolute value of the
-                                     * criterion is less than the
-                                     * given threshold.
-                                     *
-                                     * Note that this function takes a vector
-                                     * of @p{float}s, rather than the usual
-                                     * @p{double}s, since accuracy is not so
-                                     * much needed here and saving memory may
-                                     * be a good goal when using many cells.
-                                     */
-    template <typename number>
-    void coarsen (const Vector<number> &criteria,
-                 const double         threshold);
-    
-                                    /**
-                                     * Refine the triangulation by refining
-                                     * a certain fraction @p{top_fraction_of_cells}
-                                     * with the highest error. Likewise coarsen
-                                     * the fraction @p{bottom_fraction_of_cells}
-                                     * with the least error. To actually
-                                     * perform the refinement, call
-                                     * @p{execute_coarsening_and_refinement}.
-                                     *
-                                     * @p{fraction_of_cells} shall be a value
-                                     * between zero and one.
-                                     *
-                                     * Refer to the general doc of this class
-                                     * for more information.
-                                     *
-                                     * Note that this function takes a vector
-                                     * of @p{float}s, rather than the usual
-                                     * @p{double}s, since accuracy is not so
-                                     * much needed here and saving memory may
-                                     * be a good goal when using many cells.
-                                     */
-    template <typename number>
-    void refine_and_coarsen_fixed_number (const Vector<number> &criteria,
-                                         const double         top_fraction_of_cells,
-                                         const double         bottom_fraction_of_cells);
-
-                                    /**
-                                     * Refine the triangulation by flagging
-                                     * those cells which make up a certain
-                                     * @p{top_fraction} of the total error.
-                                     * Likewise, coarsen all cells which
-                                     * make up only @p{bottom_fraction}.
-                                     * To actually perform the refinement, call
-                                     * @p{execute_coarsening_and_refinement}.
-                                     *
-                                     * @p{*_fraction} shall be a values
-                                     * between zero and one.
-                                     *
-                                     * Refer to the general doc of this class
-                                     * for more information.
-                                     *
-                                     * Note that this function takes a vector
-                                     * of @p{float}s, rather than the usual
-                                     * @p{double}s, since accuracy is not so
-                                     * much needed here and saving memory may
-                                     * be a good goal when using many cells.
-                                     */
-    template<typename number>
-    void refine_and_coarsen_fixed_fraction (const Vector<number> &criteria,
-                                           const double         top_fraction,
-                                           const double         bottom_fraction);
-
-                                    /**
-                                     * Execute both refinement and coarsening
-                                     * of the triangulation.
+                                     * Execute both refinement and
+                                     * coarsening of the
+                                     * triangulation.
                                      *
                                      * The function resets all
                                      * refinement and coarsening
@@ -1953,60 +1795,72 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                       * See the general docs for more
                                       * information.
                                      *
-                                     * Note that this function is @p{virtual} to
-                                     * allow derived classes to insert hooks,
-                                     * such as saving refinement flags and the
-                                     * like.
+                                     * Note that this function is
+                                     * @p{virtual} to allow derived
+                                     * classes to insert hooks, such
+                                     * as saving refinement flags and
+                                     * the like.
                                      */
     virtual void execute_coarsening_and_refinement ();
     
                                     /**
-                                     * Do both preparation for refinement and
-                                     * coarsening as well as mesh smoothing.
+                                     * Do both preparation for
+                                     * refinement and coarsening as
+                                     * well as mesh smoothing.
                                      *
-                                     * Regarding the refinement process it fixes
-                                     * the closure of the refinement in @p{dim>=2}
-                                     * (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.
-                                     * The function returns whether additional
-                                     * cells have been flagged for refinement.
-                                     *  
-                                     * See the general doc of this class for
-                                     * more information on smoothing upon
+                                     * Regarding the refinement
+                                     * process it fixes the closure
+                                     * of the refinement in
+                                     * @p{dim>=2} (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.  The function
+                                     * returns whether additional
+                                     * cells have been flagged for
                                      * refinement.
+                                     *  
+                                     * See the general doc of this
+                                     * class for more information on
+                                     * smoothing upon refinement.
                                      *
-                                     * This part of the function is mostly
-                                     * dimension independent. However, for some
-                                     * dimension dependent things, it calls
+                                     * This part of the function is
+                                     * mostly dimension
+                                     * independent. However, for some
+                                     * dimension dependent things, it
+                                     * calls
                                      * @p{prepare_refinement_dim_dependent}.
                                      *
-                                     * Regarding the coarsening part, flagging
-                                     * and deflagging cells in preparation 
-                                     * of the actual coarsening step are
-                                     * done. This includes deleting coarsen
-                                      * flags from cells which may not be
-                                      * deleted (e.g. because one neighbor is
-                                      * more refined than the cell), doing
-                                      * some smoothing, etc.
+                                     * Regarding the coarsening part,
+                                     * flagging and deflagging cells
+                                     * in preparation of the actual
+                                     * coarsening step are done. This
+                                     * includes deleting coarsen 
+                                     * flags from cells which may not
+                                     * be deleted (e.g. because one
+                                     * neighbor is more refined
+                                     * than the cell), doing some
+                                     * smoothing, etc.
                                      *
-                                     * The effect is that only those cells
-                                     * are flagged for coarsening which
-                                     * will actually be coarsened. This
-                                     * includes the fact that all flagged
-                                     * cells belong to parent cells of which
-                                     * all children are flagged.
+                                     * The effect is that only those
+                                     * cells are flagged for
+                                     * coarsening which will actually
+                                     * be coarsened. This includes
+                                     * the fact that all flagged
+                                     * cells belong to parent cells
+                                     * of which all children are
+                                     * flagged.
                                      *
-                                     * The function returns whether some
-                                     * cells' flagging has been changed in
-                                     * the process.
+                                     * The function returns whether
+                                     * some cells' flagging has been
+                                     * changed in the process.
                                      *
-                                     * This function uses the user flags, so
-                                     * store them if you still need them
-                                     * afterwards.
+                                     * This function uses the user
+                                     * flags, so store them if you
+                                     * still need them afterwards.
                                      */
     bool prepare_coarsening_and_refinement ();
     
@@ -2017,9 +1871,10 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                      */
                                     /*@{*/
                                     /**
-                                     *  Save the addresses of the cells which
-                                     *  are flagged for refinement to @p{out}.
-                                     *  For usage, read the general
+                                     *  Save the addresses of the
+                                     *  cells which are flagged for
+                                     *  refinement to @p{out}.  For
+                                     *  usage, read the general
                                      *  documentation for this class.
                                      */
     void save_refine_flags (ostream &out) const;
@@ -3007,17 +2862,6 @@ class Triangulation : public TriaDimensionInfo<dim>,
                                     /**
                                      * Exception
                                      */
-    DeclException2 (ExcInvalidVectorSize,
-                   int, int,
-                   << "The given vector has " << arg1
-                   << " elements, but " << arg2 << " were expected.");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInvalidParameterValue);
-                                    /**
-                                     * Exception
-                                     */
     DeclException0 (ExcIO);
                                     /*@}*/
   protected:
diff --git a/deal.II/deal.II/source/grid/grid_refinement.cc b/deal.II/deal.II/source/grid/grid_refinement.cc
new file mode 100644 (file)
index 0000000..2ff5f96
--- /dev/null
@@ -0,0 +1,235 @@
+//----------------------------  grid_refinement.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  grid_refinement.cc  ---------------------------
+
+
+#include <lac/vector.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria.h>
+
+#include <numeric>
+#include <algorithm>
+
+
+
+
+template <int dim, typename number>
+void GridRefinement::refine (Triangulation<dim>   &tria,
+                            const Vector<number> &criteria,
+                            const double         threshold)
+{
+  Assert (criteria.size() == tria.n_active_cells(),
+         ExcInvalidVectorSize(criteria.size(), tria.n_active_cells()));
+
+                                  // nothing to do; especially we
+                                  // do not want to flag with zero
+                                  // error since then we may get
+                                  // into conflict with coarsening
+                                  // in some cases
+  if (threshold==0)
+    return;
+  
+  Triangulation<dim>::active_cell_iterator cell = tria.begin_active();
+  const unsigned int n_cells = criteria.size();
+  
+  for (unsigned int index=0; index<n_cells; ++cell, ++index)
+    if (fabs(criteria(index)) >= threshold)
+      cell->set_refine_flag();
+};
+
+
+
+template <int dim, typename number>
+void GridRefinement::coarsen (Triangulation<dim>   &tria,
+                             const Vector<number> &criteria,
+                             const double         threshold)
+{
+  Assert (criteria.size() == tria.n_active_cells(),
+         ExcInvalidVectorSize(criteria.size(), tria.n_active_cells()));
+
+  Triangulation<dim>::active_cell_iterator cell = tria.begin_active();
+  const unsigned int n_cells = criteria.size();
+  
+  for (unsigned int index=0; index<n_cells; ++cell, ++index)
+    if (fabs(criteria(index)) <= threshold)
+      cell->set_coarsen_flag();
+};
+
+
+
+template <int dim, typename number>
+void
+GridRefinement::refine_and_coarsen_fixed_number (Triangulation<dim>   &tria,
+                                                const Vector<number> &criteria,
+                                                const double         top_fraction,
+                                                const double         bottom_fraction)
+{
+                                  // correct number of cells is
+                                  // checked in @p{refine}
+  Assert ((top_fraction>=0) && (top_fraction<=1), ExcInvalidParameterValue());
+  Assert ((bottom_fraction>=0) && (bottom_fraction<=1), ExcInvalidParameterValue());
+  Assert (top_fraction+bottom_fraction <= 1, ExcInvalidParameterValue());
+
+  const int refine_cells=static_cast<int>(top_fraction*criteria.size());
+  const int coarsen_cells=static_cast<int>(bottom_fraction*criteria.size());
+
+  if (refine_cells || coarsen_cells)
+    {
+      Vector<number> tmp(criteria);
+      if (refine_cells)
+       {
+         nth_element (tmp.begin(), tmp.begin()+refine_cells,
+                      tmp.end(),
+                      greater<double>());
+         refine (tria, criteria, *(tmp.begin() + refine_cells));
+       };
+
+      if (coarsen_cells)
+       {
+         nth_element (tmp.begin(), tmp.begin()+tmp.size()-coarsen_cells,
+                      tmp.end(),
+                      greater<double>());
+         coarsen (tria, criteria, *(tmp.begin() + tmp.size() - coarsen_cells));
+       };
+    };
+};
+
+
+
+template <int dim, typename number>
+void
+GridRefinement::refine_and_coarsen_fixed_fraction (Triangulation<dim>   &tria,
+                                                  const Vector<number> &criteria,
+                                                  const double         top_fraction,
+                                                  const double         bottom_fraction)
+{
+                                  // correct number of cells is
+                                  // checked in @p{refine}
+  Assert ((top_fraction>=0) && (top_fraction<=1), ExcInvalidParameterValue());
+  Assert ((bottom_fraction>=0) && (bottom_fraction<=1), ExcInvalidParameterValue());
+  Assert (top_fraction+bottom_fraction <= 1, ExcInvalidParameterValue());
+
+                                  // let tmp be the cellwise square of the
+                                  // error, which is what we have to sum
+                                  // up and compare with
+                                  // @p{fraction_of_error*total_error}.
+  Vector<number> tmp(criteria);
+  const double total_error = tmp.l1_norm();
+
+  Vector<float> partial_sums(criteria.size());
+  
+                                  // sort the largest criteria to the
+                                  // beginning of the vector
+  sort (tmp.begin(), tmp.end(), greater<double>());
+  partial_sum (tmp.begin(), tmp.end(), partial_sums.begin());
+
+                                  // compute thresholds
+  const Vector<float>::const_iterator
+    q = lower_bound (partial_sums.begin(), partial_sums.end(),
+                    top_fraction*total_error),
+    p = upper_bound (partial_sums.begin(), partial_sums.end(),
+                    total_error*(1-bottom_fraction));
+  
+  double bottom_threshold = tmp(p != partial_sums.end() ?
+                               p-partial_sums.begin() :
+                               criteria.size()-1),
+        top_threshold    = tmp(q-partial_sums.begin());
+
+                                  // in some rare cases it may happen that
+                                  // both thresholds are the same (e.g. if
+                                  // there are many cells with the same
+                                  // error indicator). That would mean that
+                                  // all cells will be flagged for
+                                  // refinement or coarsening, but some will
+                                  // be flagged for both, namely those for
+                                  // which the indicator equals the
+                                  // thresholds. This is forbidden, however.
+                                  //
+                                  // In some rare cases with very few cells
+                                  // we also could get integer round off
+                                  // errors and get problems with
+                                  // the top and bottom fractions.
+                                  //
+                                  // In these case we arbitrarily reduce the
+                                  // bottom threshold by one permille below
+                                  // the top threshold
+                                  //
+                                  // Finally, in some cases
+                                  // (especially involving symmetric
+                                  // solutions) there are many cells
+                                  // with the same error indicator
+                                  // values. if there are many with
+                                  // indicator equal to the top
+                                  // threshold, no refinement will
+                                  // take place below; to avoid this
+                                  // case, we also lower the top
+                                  // threshold if it equals the
+                                  // largest indicator and the
+                                  // top_fraction!=1
+  if ((top_threshold == *max_element(criteria.begin(), criteria.end())) &&
+      (top_fraction != 1))
+    top_threshold *= 0.999;
+  
+  if (bottom_threshold>=top_threshold)
+    bottom_threshold = 0.999*top_threshold;
+  
+                                  // actually flag cells
+  if (top_threshold < *max_element(criteria.begin(), criteria.end()))
+    refine (tria, criteria, top_threshold);
+  if (bottom_threshold > *min_element(criteria.begin(), criteria.end()))
+    coarsen (tria, criteria, bottom_threshold);
+};
+
+
+
+
+
+// explicit instantiations
+template void GridRefinement
+::refine (Triangulation<deal_II_dimension> &, const Vector<float> &, const double);
+
+template void GridRefinement
+::refine (Triangulation<deal_II_dimension> &, const Vector<double> &, const double);
+
+template void GridRefinement
+::coarsen (Triangulation<deal_II_dimension> &, const Vector<float> &, const double);
+
+template void GridRefinement
+::coarsen (Triangulation<deal_II_dimension> &, const Vector<double> &, const double);
+
+
+template void GridRefinement
+::refine_and_coarsen_fixed_number (Triangulation<deal_II_dimension> &,
+                                  const Vector<double> &,
+                                  const double         top_fraction,
+                                  const double         bottom_fraction);
+
+template void GridRefinement
+::refine_and_coarsen_fixed_number (Triangulation<deal_II_dimension> &,
+                                  const Vector<float> &criteria,
+                                  const double         top_fraction,
+                                  const double         bottom_fraction);
+
+template void GridRefinement
+::refine_and_coarsen_fixed_fraction (Triangulation<deal_II_dimension> &,
+                                    const Vector<double> &criteria,
+                                    const double         top_fraction,
+                                    const double         bottom_fraction);
+
+template void GridRefinement
+::refine_and_coarsen_fixed_fraction (Triangulation<deal_II_dimension> &,
+                                    const Vector<float> &criteria,
+                                    const double         top_fraction,
+                                    const double         bottom_fraction);
+
index 8bdc92e9603927b6ece1e1c7f7b61f5b688bcd0e..09b176a26da57d7a51323dc05af756c439f2e7cf 100644 (file)
@@ -3254,178 +3254,6 @@ unsigned int Triangulation<dim>::max_adjacent_cells () const {
 };
 
 
-template <int dim>
-template <typename number>
-void Triangulation<dim>::refine (const Vector<number> &criteria,
-                                const double         threshold)
-{
-  Assert (criteria.size() == n_active_cells(),
-         ExcInvalidVectorSize(criteria.size(), n_active_cells()));
-
-                                  // nothing to do; especially we
-                                  // do not want to flag with zero
-                                  // error since then we may get
-                                  // into conflict with coarsening
-                                  // in some cases
-  if (threshold==0)
-    return;
-  
-  active_cell_iterator cell = begin_active();
-  const unsigned int n_cells = criteria.size();
-  
-  for (unsigned int index=0; index<n_cells; ++cell, ++index)
-    if (fabs(criteria(index)) >= threshold)
-      cell->set_refine_flag();
-};
-
-
-template <int dim>
-template <typename number>
-void Triangulation<dim>::coarsen (const Vector<number> &criteria,
-                                 const double         threshold)
-{
-  Assert (criteria.size() == n_active_cells(),
-         ExcInvalidVectorSize(criteria.size(), n_active_cells()));
-
-  active_cell_iterator cell = begin_active();
-  const unsigned int n_cells = criteria.size();
-  
-  for (unsigned int index=0; index<n_cells; ++cell, ++index)
-    if (fabs(criteria(index)) <= threshold)
-      cell->set_coarsen_flag();
-};
-
-
-template <int dim>
-template <typename number>
-void
-Triangulation<dim>::refine_and_coarsen_fixed_number (const Vector<number> &criteria,
-                                                    const double         top_fraction,
-                                                    const double         bottom_fraction)
-{
-                                  // correct number of cells is
-                                  // checked in @p{refine}
-  Assert ((top_fraction>=0) && (top_fraction<=1), ExcInvalidParameterValue());
-  Assert ((bottom_fraction>=0) && (bottom_fraction<=1), ExcInvalidParameterValue());
-  Assert (top_fraction+bottom_fraction <= 1, ExcInvalidParameterValue());
-
-  const int refine_cells=static_cast<int>(top_fraction*criteria.size());
-  const int coarsen_cells=static_cast<int>(bottom_fraction*criteria.size());
-
-  if (refine_cells || coarsen_cells)
-    {
-      Vector<number> tmp(criteria);
-      if (refine_cells)
-       {
-         nth_element (tmp.begin(), tmp.begin()+refine_cells,
-                      tmp.end(),
-                      greater<double>());
-         refine (criteria, *(tmp.begin() + refine_cells));
-       }
-
-      if (coarsen_cells)
-       {
-         nth_element (tmp.begin(), tmp.begin()+tmp.size()-coarsen_cells,
-                      tmp.end(),
-                      greater<double>());
-         coarsen (criteria, *(tmp.begin() + tmp.size() - coarsen_cells));
-       }
-    }
-};
-
-
-static
-inline
-double sqr(double a) {
-  return a*a;
-};
-
-
-template <int dim>
-template <typename number>
-void
-Triangulation<dim>::refine_and_coarsen_fixed_fraction (const Vector<number> &criteria,
-                                                      const double         top_fraction,
-                                                      const double         bottom_fraction) {
-                                  // correct number of cells is
-                                  // checked in @p{refine}
-  Assert ((top_fraction>=0) && (top_fraction<=1), ExcInvalidParameterValue());
-  Assert ((bottom_fraction>=0) && (bottom_fraction<=1), ExcInvalidParameterValue());
-  Assert (top_fraction+bottom_fraction <= 1, ExcInvalidParameterValue());
-
-                                  // let tmp be the cellwise square of the
-                                  // error, which is what we have to sum
-                                  // up and compare with
-                                  // @p{fraction_of_error*total_error}.
-  Vector<number> tmp(criteria);
-  const double total_error = tmp.l1_norm();
-
-  Vector<float> partial_sums(criteria.size());
-  
-                                  // sort the largest criteria to the
-                                  // beginning of the vector
-  sort (tmp.begin(), tmp.end(), greater<double>());
-  partial_sum (tmp.begin(), tmp.end(), partial_sums.begin());
-
-                                  // compute thresholds
-  const Vector<float>::const_iterator
-    q = lower_bound (partial_sums.begin(), partial_sums.end(),
-                    top_fraction*total_error),
-    p = upper_bound (partial_sums.begin(), partial_sums.end(),
-                    total_error*(1-bottom_fraction));
-  
-  double bottom_threshold = tmp(p != partial_sums.end() ?
-                               p-partial_sums.begin() :
-                               criteria.size()-1),
-        top_threshold    = tmp(q-partial_sums.begin());
-
-                                  // in some rare cases it may happen that
-                                  // both thresholds are the same (e.g. if
-                                  // there are many cells with the same
-                                  // error indicator). That would mean that
-                                  // all cells will be flagged for
-                                  // refinement or coarsening, but some will
-                                  // be flagged for both, namely those for
-                                  // which the indicator equals the
-                                  // thresholds. This is forbidden, however.
-                                  //
-                                  // In some rare cases with very few cells
-                                  // we also could get integer round off
-                                  // errors and get problems with
-                                  // the top and bottom fractions.
-                                  //
-                                  // In these case we arbitrarily reduce the
-                                  // bottom threshold by one permille below
-                                  // the top threshold
-                                  //
-                                  // Finally, in some cases
-                                  // (especially involving symmetric
-                                  // solutions) there are many cells
-                                  // with the same error indicator
-                                  // values. if there are many with
-                                  // indicator equal to the top
-                                  // threshold, no refinement will
-                                  // take place below; to avoid this
-                                  // case, we also lower the top
-                                  // threshold if it equals the
-                                  // largest indicator and the
-                                  // top_fraction!=1
-  if ((top_threshold == *max_element(criteria.begin(), criteria.end())) &&
-      (top_fraction != 1))
-    top_threshold *= 0.999;
-  
-  if (bottom_threshold>=top_threshold)
-    bottom_threshold = 0.999*top_threshold;
-  
-                                  // actually flag cells
-  if (top_threshold < *max_element(criteria.begin(), criteria.end()))
-    refine (criteria, top_threshold);
-  if (bottom_threshold > *min_element(criteria.begin(), criteria.end()))
-    coarsen (criteria, bottom_threshold);
-
-  prepare_coarsening_and_refinement ();
-};
-
 
 template <int dim>
 void Triangulation<dim>::execute_coarsening_and_refinement () {
@@ -7387,36 +7215,3 @@ void Triangulation<dim>::update_number_cache ()
 // explicit instantiations
 template class Triangulation<deal_II_dimension>;
 
-template void Triangulation<deal_II_dimension>
-::refine (const Vector<float> &, const double);
-
-template void Triangulation<deal_II_dimension>
-::refine (const Vector<double> &, const double);
-
-template void Triangulation<deal_II_dimension>
-::coarsen (const Vector<float> &, const double);
-
-template void Triangulation<deal_II_dimension>
-::coarsen (const Vector<double> &, const double);
-
-
-template void Triangulation<deal_II_dimension>
-::refine_and_coarsen_fixed_number (const Vector<double> &,
-                                  const double         top_fraction,
-                                  const double         bottom_fraction);
-
-template void Triangulation<deal_II_dimension>
-::refine_and_coarsen_fixed_number (const Vector<float> &criteria,
-                                  const double         top_fraction,
-                                  const double         bottom_fraction);
-
-template void Triangulation<deal_II_dimension>
-::refine_and_coarsen_fixed_fraction (const Vector<double> &criteria,
-                                    const double         top_fraction,
-                                    const double         bottom_fraction);
-
-template void Triangulation<deal_II_dimension>
-::refine_and_coarsen_fixed_fraction (const Vector<float> &criteria,
-                                    const double         top_fraction,
-                                    const double         bottom_fraction);
-
index b82cbe520752ae8bf75f8143d963896bf86a77a7..b05e6b3781a0d3094bc5dcf5490109b128a257a6 100644 (file)
@@ -16,6 +16,7 @@
 #include <grid/tria.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
+#include <grid/grid_refinement.h>
 #include <lac/vector.h>
 
 #ifdef DEAL_II_USE_MT
@@ -711,8 +712,8 @@ void TimeStepBase_Tria<dim>::refine_grid (const RefinementData refinement_data)
 
 
                                   // actually flag cells the first time
-  tria->refine (criteria, refinement_threshold);
-  tria->coarsen (criteria, coarsening_threshold);
+  GridRefinement::refine (*tria, criteria, refinement_threshold);
+  GridRefinement::coarsen (*tria, criteria, coarsening_threshold);
 
                                   // store this number for the following
                                   // since its computation is rather
@@ -996,8 +997,8 @@ void TimeStepBase_Tria<dim>::refine_grid (const RefinementData refinement_data)
 
 
                                         // flag cells finally
-       tria->refine (criteria, refinement_threshold);
-       tria->coarsen (criteria, coarsening_threshold);
+       GridRefinement::refine (*tria, criteria, refinement_threshold);
+       GridRefinement::coarsen (*tria, criteria, coarsening_threshold);
       };
   
                                   // if step number is greater than
index 26f74bccfbcdb131e49cddb2d1b00bf0257a060a..42e56590b047978aecfc674dc6cc94d5d5567e8a 100644 (file)
                                 // instead of ``grid_in.h'':
 #include <grid/grid_out.h>
 
+                                // In order to refine our grids
+                                // locally, we need a function from
+                                // the library that decides which
+                                // cells to flag for refinement or
+                                // coarsening based on the error
+                                // indicators we have computed. This
+                                // function is defined here:
+#include <grid/grid_refinement.h>
+
                                 // When using locally refined grids,
                                 // we will get so-called ``hanging
                                 // nodes''. However, the standard
@@ -760,8 +769,20 @@ void LaplaceProblem<dim>::refine_grid ()
                                   // over-refinement may have taken
                                   // place. Thus a small, non-zero
                                   // value is appropriate here.
-  triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
-                                                0.3, 0.03);
+                                  //
+                                  // The following function now takes
+                                  // these refinement indicators and
+                                  // flags some cells of the
+                                  // triangulation for refinement or
+                                  // coarsening using the method
+                                  // described above. It is from a
+                                  // class that implements
+                                  // several different algorithms to
+                                  // refine a triangulation based on
+                                  // cellwise error indicators.
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                  estimated_error_per_cell,
+                                                  0.3, 0.03);
 
                                   // After the previous function has
                                   // exited, some cells are flagged
index 4df187039733dc8f318f44737f20193f04ca0c36..e2e35aec01f51823aa7dced0db636b96efe05fd2 100644 (file)
@@ -16,6 +16,7 @@
 #include <lac/precondition.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary_lib.h>
@@ -1086,8 +1087,9 @@ void LaplaceProblem<dim>::refine_grid ()
                                            solution,
                                            estimated_error_per_cell);
 
-       triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
-                                                      0.3, 0.03);
+       GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                        estimated_error_per_cell,
+                                                        0.3, 0.03);
        
        triangulation.execute_coarsening_and_refinement ();
 
index b03ead2eb615ce8a89fee6f12f6447926d281e71..264a93a8bb3712bdd0b8b16f328fe94d707f9256 100644 (file)
@@ -16,6 +16,7 @@
 #include <grid/tria.h>
 #include <dofs/dof_handler.h>
 #include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary_lib.h>
@@ -800,8 +801,9 @@ void ElasticProblem<dim>::refine_grid ()
                                      solution,
                                      estimated_error_per_cell);
 
-  triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
-                                                0.3, 0.03);
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                  estimated_error_per_cell,
+                                                  0.3, 0.03);
 
   triangulation.execute_coarsening_and_refinement ();
 };
@@ -983,7 +985,7 @@ void ElasticProblem<dim>::run ()
                                           // cell as well, and the
                                           // call to
                                           // ``refine_and_coarsen_fixed_number''
-                                          // of the ``triangulation''
+                                          // on the ``triangulation''
                                           // object will not flag any
                                           // cells for refinement
                                           // (why should it if the
@@ -1011,7 +1013,27 @@ void ElasticProblem<dim>::run ()
                                           // it needs to be able to
                                           // see the right hand
                                           // side. Thus, we refine
-                                          // twice globally.
+                                          // twice globally. (Note
+                                          // that the
+                                          // ``refine_global''
+                                          // function is not part of
+                                          // the ``GridRefinement''
+                                          // class in which
+                                          // ``refine_and_coarsen_fixed_number''
+                                          // is declared, for
+                                          // example. The reason is
+                                          // first that it is not an
+                                          // algorithm that computed
+                                          // refinement flags from
+                                          // indicators, but more
+                                          // importantly that it
+                                          // actually performs the
+                                          // refinement, in contrast
+                                          // to the functions in
+                                          // ``GridRefinement'' that
+                                          // only flag cells without
+                                          // actually refining the
+                                          // grid.)
          triangulation.refine_global (2);
        }
       else
index dc169920c8cdbb613396c81921dbba712e7b8570..7e4eae0a6ce56e096a85219bf6e33f6e8d3bda1c 100644 (file)
@@ -16,6 +16,7 @@
 #include <lac/precondition.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary_lib.h>
@@ -1329,8 +1330,9 @@ void AdvectionProblem<dim>::refine_grid ()
                                solution,
                                estimated_error_per_cell);
 
-  triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
-                                                0.5, 0.03);
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                  estimated_error_per_cell,
+                                                  0.5, 0.03);
 
   triangulation.execute_coarsening_and_refinement ();
 };
index 6d181e222ea536ab163d0f4ad2f971bf9ee348a9..8c8f107312b733fc0408809bc316cd50d26c6f81 100644 (file)
@@ -11,6 +11,7 @@
 #include <grid/tria_accessor.h>
 #include <grid/tria_boundary_lib.h>
 #include <grid/tria_iterator.h>
+#include <grid/grid_refinement.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_constraints.h>
 #include <dofs/dof_handler.h>
@@ -623,14 +624,18 @@ void PoissonProblem<dim>::run (ParameterHandler &prm) {
                tria->set_all_refine_flags ();
                break;
          case true_error:
-               tria->refine_and_coarsen_fixed_number (h1_error_per_cell,
-                                                      prm.get_double("Refinement fraction"),
-                                                      prm.get_double("Coarsening fraction"));
+               GridRefinement::
+                 refine_and_coarsen_fixed_number (*tria,
+                                                  h1_error_per_cell,
+                                                  prm.get_double("Refinement fraction"),
+                                                  prm.get_double("Coarsening fraction"));
                break;
          case error_estimator:
-               tria->refine_and_coarsen_fixed_number (estimated_error_per_cell,
-                                                      prm.get_double("Refinement fraction"),
-                                                      prm.get_double("Coarsening fraction"));
+               GridRefinement::
+                 refine_and_coarsen_fixed_number (*tria,
+                                                  estimated_error_per_cell,
+                                                  prm.get_double("Refinement fraction"),
+                                                  prm.get_double("Coarsening fraction"));
                break;
        };
 
index d28e210a3e5ee3fdaa02e32ab759237bd5d92e1f..477465cab8e21172d462cd34518913dc5d5a72cb 100644 (file)
@@ -13,6 +13,7 @@
 #include <lac/vector_memory.h>
 #include <lac/precondition.h>
 #include <grid/tria.h>
+#include <grid/grid_refinement.h>
 #include <dofs/dof_handler.h>
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
@@ -404,8 +405,9 @@ void LaplaceProblem<dim>::refine_grid ()
                                      solution,
                                      estimated_error_per_cell);
 
-  triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
-                                                0.3, 0.03);
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                  estimated_error_per_cell,
+                                                  0.3, 0.03);
   triangulation.execute_coarsening_and_refinement ();
 };
 
index 1ab5f1b2a9bd6f369559853d6834fa66a7040185..aaf36230c8c67f3806211f547fd8ef13d766cc3b 100644 (file)
@@ -12,6 +12,7 @@
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary.h>
 #include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
 #include <numerics/data_out.h>
 #include <base/function.h>
 #include <fe/fe_lib.lagrange.h>
@@ -235,7 +236,7 @@ void NonlinearProblem<dim>::run () {
                   KellyErrorEstimator<dim>::FunctionMap(),
                   solution,
                   error_indicator);
-      tria->refine_and_coarsen_fixed_number (error_indicator, 0.3, 0);
+      GridRefinement::refine_and_coarsen_fixed_number (*tria, error_indicator, 0.3, 0);
       tria->execute_coarsening_and_refinement ();
     };
   

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.