* 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
* @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
* @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.
* 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
{
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,
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,
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,
* 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.
* 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,
--- /dev/null
+//---------------------------- 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
/**
- * 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
/**
- * 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
* 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
* 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
* 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
* 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}
* @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.
* 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
* 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}
* 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
* 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
* 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 () {
* @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.
*
*
{
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;
* for mesh smoothing
* algorithms. The meaning of
* these flags is documented in
- * the @p{Triangulation} class.
+ * the @ref{Triangulation} class.
*/
enum MeshSmoothing
{
*
* 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 ();
* 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
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,
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,
*/
/*@{*/
/**
- * 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
* 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 ();
*/
/*@{*/
/**
- * 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;
/**
* Exception
*/
- DeclException2 (ExcInvalidVectorSize,
- int, int,
- << "The given vector has " << arg1
- << " elements, but " << arg2 << " were expected.");
- /**
- * Exception
- */
- DeclException0 (ExcInvalidParameterValue);
- /**
- * Exception
- */
DeclException0 (ExcIO);
/*@}*/
protected:
--- /dev/null
+//---------------------------- 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);
+
};
-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 () {
// 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);
-
#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
// 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
// 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
// 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
// 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
#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>
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 ();
#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>
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 ();
};
// 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
// 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
#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>
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 ();
};
#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>
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;
};
#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>
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 ();
};
#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>
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 ();
};
projects / dealii.git / commitdiff