From 473193dfd70136b01892590de670eaa1a6a75a8c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 28 Jul 2000 14:27:29 +0000 Subject: [PATCH] Move the grid refinement functions from the triangulation class to a class of its own. git-svn-id: https://svn.dealii.org/trunk@3219 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_generator.h | 83 +-- .../deal.II/include/grid/grid_refinement.h | 266 +++++++++ deal.II/deal.II/include/grid/tria.h | 514 ++++++------------ .../deal.II/source/grid/grid_refinement.cc | 235 ++++++++ deal.II/deal.II/source/grid/tria.cc | 205 ------- .../deal.II/source/numerics/time_dependent.cc | 9 +- deal.II/examples/step-6/step-6.cc | 25 +- deal.II/examples/step-7/step-7.cc | 6 +- deal.II/examples/step-8/step-8.cc | 30 +- deal.II/examples/step-9/step-9.cc | 6 +- .../error-estimation/error-estimation.cc | 17 +- tests/big-tests/multigrid/multigrid.cc | 6 +- .../fixed-point-iteration/nonlinear.cc | 3 +- 13 files changed, 808 insertions(+), 597 deletions(-) create mode 100644 deal.II/deal.II/include/grid/grid_refinement.h create mode 100644 deal.II/deal.II/source/grid/grid_refinement.cc diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index 25d2f04827..f37d336a49 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -29,12 +29,12 @@ template 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 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 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 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 static void hyper_cube (Triangulation &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 static void hyper_ball (Triangulation &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 static void hyper_L (Triangulation &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 static void hyper_shell (Triangulation &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 index 0000000000..b93e3bfb7e --- /dev/null +++ b/deal.II/deal.II/include/grid/grid_refinement.h @@ -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 class Triangulation; +template 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 &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 &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 + static void refine (Triangulation &tria, + const Vector &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 + static void coarsen (Triangulation &tria, + const Vector &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 + static void refine_and_coarsen_fixed_number (Triangulation &tria, + const Vector &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 + static void refine_and_coarsen_fixed_fraction (Triangulation &tria, + const Vector &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 diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index f4d5422639..d9df9075a6 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -34,12 +34,13 @@ template class MGDoFHandler; /** - * Structure which is passed to the @p{Triangulation::create_triangulation} + * Structure which is passed to the @ref{Triangulation}@p{::create_triangulation} * function. It contains all data needed to construct a cell, namely the * indices of the vertices and the material indicator. */ template -struct CellData { +struct CellData +{ #if !((__GNUC__==2) && (__GNUC_MINOR__==95)) int vertices[GeometryInfo::vertices_per_cell]; #else @@ -50,7 +51,7 @@ struct CellData { /** - * Structure to be passed to the @p{Triangulation::create_triangulation} + * Structure to be passed to the @ref{Triangulation}@p{::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::read_*} functions. The appropriate function to be + * @ref{GridIn}@p{::read_*} functions. The appropriate function to be * called is @p{Triangulation::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 &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 &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, { 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& straight_boundary; @@ -1583,7 +1509,7 @@ class Triangulation : public TriaDimensionInfo, * 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, * * 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, * 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, const Boundary &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 & 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 &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< > &vertices, @@ -1787,18 +1726,21 @@ class Triangulation : public TriaDimensionInfo, 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, */ /*@{*/ /** - * 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 - void refine (const Vector &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 - void coarsen (const Vector &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 - void refine_and_coarsen_fixed_number (const Vector &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 - void refine_and_coarsen_fixed_fraction (const Vector &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, * 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, */ /*@{*/ /** - * 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, /** * 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 index 0000000000..2ff5f96d67 --- /dev/null +++ b/deal.II/deal.II/source/grid/grid_refinement.cc @@ -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 +#include +#include +#include +#include + +#include +#include + + + + +template +void GridRefinement::refine (Triangulation &tria, + const Vector &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::active_cell_iterator cell = tria.begin_active(); + const unsigned int n_cells = criteria.size(); + + for (unsigned int index=0; index= threshold) + cell->set_refine_flag(); +}; + + + +template +void GridRefinement::coarsen (Triangulation &tria, + const Vector &criteria, + const double threshold) +{ + Assert (criteria.size() == tria.n_active_cells(), + ExcInvalidVectorSize(criteria.size(), tria.n_active_cells())); + + Triangulation::active_cell_iterator cell = tria.begin_active(); + const unsigned int n_cells = criteria.size(); + + for (unsigned int index=0; indexset_coarsen_flag(); +}; + + + +template +void +GridRefinement::refine_and_coarsen_fixed_number (Triangulation &tria, + const Vector &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(top_fraction*criteria.size()); + const int coarsen_cells=static_cast(bottom_fraction*criteria.size()); + + if (refine_cells || coarsen_cells) + { + Vector tmp(criteria); + if (refine_cells) + { + nth_element (tmp.begin(), tmp.begin()+refine_cells, + tmp.end(), + greater()); + refine (tria, criteria, *(tmp.begin() + refine_cells)); + }; + + if (coarsen_cells) + { + nth_element (tmp.begin(), tmp.begin()+tmp.size()-coarsen_cells, + tmp.end(), + greater()); + coarsen (tria, criteria, *(tmp.begin() + tmp.size() - coarsen_cells)); + }; + }; +}; + + + +template +void +GridRefinement::refine_and_coarsen_fixed_fraction (Triangulation &tria, + const Vector &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 tmp(criteria); + const double total_error = tmp.l1_norm(); + + Vector partial_sums(criteria.size()); + + // sort the largest criteria to the + // beginning of the vector + sort (tmp.begin(), tmp.end(), greater()); + partial_sum (tmp.begin(), tmp.end(), partial_sums.begin()); + + // compute thresholds + const Vector::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 &, const Vector &, const double); + +template void GridRefinement +::refine (Triangulation &, const Vector &, const double); + +template void GridRefinement +::coarsen (Triangulation &, const Vector &, const double); + +template void GridRefinement +::coarsen (Triangulation &, const Vector &, const double); + + +template void GridRefinement +::refine_and_coarsen_fixed_number (Triangulation &, + const Vector &, + const double top_fraction, + const double bottom_fraction); + +template void GridRefinement +::refine_and_coarsen_fixed_number (Triangulation &, + const Vector &criteria, + const double top_fraction, + const double bottom_fraction); + +template void GridRefinement +::refine_and_coarsen_fixed_fraction (Triangulation &, + const Vector &criteria, + const double top_fraction, + const double bottom_fraction); + +template void GridRefinement +::refine_and_coarsen_fixed_fraction (Triangulation &, + const Vector &criteria, + const double top_fraction, + const double bottom_fraction); + diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 8bdc92e960..09b176a26d 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -3254,178 +3254,6 @@ unsigned int Triangulation::max_adjacent_cells () const { }; -template -template -void Triangulation::refine (const Vector &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= threshold) - cell->set_refine_flag(); -}; - - -template -template -void Triangulation::coarsen (const Vector &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; indexset_coarsen_flag(); -}; - - -template -template -void -Triangulation::refine_and_coarsen_fixed_number (const Vector &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(top_fraction*criteria.size()); - const int coarsen_cells=static_cast(bottom_fraction*criteria.size()); - - if (refine_cells || coarsen_cells) - { - Vector tmp(criteria); - if (refine_cells) - { - nth_element (tmp.begin(), tmp.begin()+refine_cells, - tmp.end(), - greater()); - refine (criteria, *(tmp.begin() + refine_cells)); - } - - if (coarsen_cells) - { - nth_element (tmp.begin(), tmp.begin()+tmp.size()-coarsen_cells, - tmp.end(), - greater()); - coarsen (criteria, *(tmp.begin() + tmp.size() - coarsen_cells)); - } - } -}; - - -static -inline -double sqr(double a) { - return a*a; -}; - - -template -template -void -Triangulation::refine_and_coarsen_fixed_fraction (const Vector &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 tmp(criteria); - const double total_error = tmp.l1_norm(); - - Vector partial_sums(criteria.size()); - - // sort the largest criteria to the - // beginning of the vector - sort (tmp.begin(), tmp.end(), greater()); - partial_sum (tmp.begin(), tmp.end(), partial_sums.begin()); - - // compute thresholds - const Vector::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 void Triangulation::execute_coarsening_and_refinement () { @@ -7387,36 +7215,3 @@ void Triangulation::update_number_cache () // explicit instantiations template class Triangulation; -template void Triangulation -::refine (const Vector &, const double); - -template void Triangulation -::refine (const Vector &, const double); - -template void Triangulation -::coarsen (const Vector &, const double); - -template void Triangulation -::coarsen (const Vector &, const double); - - -template void Triangulation -::refine_and_coarsen_fixed_number (const Vector &, - const double top_fraction, - const double bottom_fraction); - -template void Triangulation -::refine_and_coarsen_fixed_number (const Vector &criteria, - const double top_fraction, - const double bottom_fraction); - -template void Triangulation -::refine_and_coarsen_fixed_fraction (const Vector &criteria, - const double top_fraction, - const double bottom_fraction); - -template void Triangulation -::refine_and_coarsen_fixed_fraction (const Vector &criteria, - const double top_fraction, - const double bottom_fraction); - diff --git a/deal.II/deal.II/source/numerics/time_dependent.cc b/deal.II/deal.II/source/numerics/time_dependent.cc index b82cbe5207..b05e6b3781 100644 --- a/deal.II/deal.II/source/numerics/time_dependent.cc +++ b/deal.II/deal.II/source/numerics/time_dependent.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #ifdef DEAL_II_USE_MT @@ -711,8 +712,8 @@ void TimeStepBase_Tria::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::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 diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc index 26f74bccfb..42e56590b0 100644 --- a/deal.II/examples/step-6/step-6.cc +++ b/deal.II/examples/step-6/step-6.cc @@ -45,6 +45,15 @@ // instead of ``grid_in.h'': #include + // 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 + // When using locally refined grids, // we will get so-called ``hanging // nodes''. However, the standard @@ -760,8 +769,20 @@ void LaplaceProblem::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 diff --git a/deal.II/examples/step-7/step-7.cc b/deal.II/examples/step-7/step-7.cc index 4df1870397..e2e35aec01 100644 --- a/deal.II/examples/step-7/step-7.cc +++ b/deal.II/examples/step-7/step-7.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -1086,8 +1087,9 @@ void LaplaceProblem::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 (); diff --git a/deal.II/examples/step-8/step-8.cc b/deal.II/examples/step-8/step-8.cc index b03ead2eb6..264a93a8bb 100644 --- a/deal.II/examples/step-8/step-8.cc +++ b/deal.II/examples/step-8/step-8.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -800,8 +801,9 @@ void ElasticProblem::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::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::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 diff --git a/deal.II/examples/step-9/step-9.cc b/deal.II/examples/step-9/step-9.cc index dc169920c8..7e4eae0a6c 100644 --- a/deal.II/examples/step-9/step-9.cc +++ b/deal.II/examples/step-9/step-9.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -1329,8 +1330,9 @@ void AdvectionProblem::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 (); }; diff --git a/tests/big-tests/error-estimation/error-estimation.cc b/tests/big-tests/error-estimation/error-estimation.cc index 6d181e222e..8c8f107312 100644 --- a/tests/big-tests/error-estimation/error-estimation.cc +++ b/tests/big-tests/error-estimation/error-estimation.cc @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -623,14 +624,18 @@ void PoissonProblem::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; }; diff --git a/tests/big-tests/multigrid/multigrid.cc b/tests/big-tests/multigrid/multigrid.cc index d28e210a3e..477465cab8 100644 --- a/tests/big-tests/multigrid/multigrid.cc +++ b/tests/big-tests/multigrid/multigrid.cc @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -404,8 +405,9 @@ void LaplaceProblem::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 (); }; diff --git a/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc b/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc index 1ab5f1b2a9..aaf36230c8 100644 --- a/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc +++ b/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -235,7 +236,7 @@ void NonlinearProblem::run () { KellyErrorEstimator::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 (); }; -- 2.39.5