KellyErrorEstimator< dim, spacedim > Class Template Reference
[Numerical algorithms]

List of all members.

Classes

class  ExcIncompatibleNumberOfElements
class  ExcInvalidBoundaryFunction
class  ExcInvalidBoundaryIndicator
class  ExcInvalidCoefficient
class  ExcInvalidComponentMask
class  ExcInvalidSolutionVector
class  ExcNoSolutions
struct  FaceIntegrals
struct  PerThreadData

Static Public Member Functions

template<typename InputVector , class DH >
static void estimate (const Mapping< dim, spacedim > &mapping, const DH &dof, const Quadrature< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const InputVector &solution, Vector< float > &error, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)
template<typename InputVector , class DH >
static void estimate (const DH &dof, const Quadrature< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const InputVector &solution, Vector< float > &error, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)
template<typename InputVector , class DH >
static void estimate (const Mapping< dim, spacedim > &mapping, const DH &dof, const Quadrature< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const std::vector< const InputVector * > &solutions, std::vector< Vector< float > * > &errors, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)
template<typename InputVector , class DH >
static void estimate (const DH &dof, const Quadrature< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const std::vector< const InputVector * > &solutions, std::vector< Vector< float > * > &errors, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)
template<typename InputVector , class DH >
static void estimate (const Mapping< dim, spacedim > &mapping, const DH &dof, const hp::QCollection< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const InputVector &solution, Vector< float > &error, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)
template<typename InputVector , class DH >
static void estimate (const DH &dof, const hp::QCollection< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const InputVector &solution, Vector< float > &error, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)
template<typename InputVector , class DH >
static void estimate (const Mapping< dim, spacedim > &mapping, const DH &dof, const hp::QCollection< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const std::vector< const InputVector * > &solutions, std::vector< Vector< float > * > &errors, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)
template<typename InputVector , class DH >
static void estimate (const DH &dof, const hp::QCollection< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const std::vector< const InputVector * > &solutions, std::vector< Vector< float > * > &errors, const std::vector< bool > &component_mask=std::vector< bool >(), const Function< dim > *coefficients=0, const unsigned int n_threads=multithread_info.n_default_threads, const unsigned int subdomain_id=numbers::invalid_unsigned_int, const unsigned int material_id=numbers::invalid_unsigned_int)

Static Private Member Functions

template<typename InputVector , class DH >
static void estimate_some (const hp::MappingCollection< dim, spacedim > &mapping, const DH &dof_handler, const hp::QCollection< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const std::vector< const InputVector * > &solutions, const std::pair< const std::vector< bool > *, const Function< dim > * > &component_mask_and_coefficients, const std::pair< unsigned int, unsigned int > this_thread, typename FaceIntegrals< DH >::type &face_integrals, PerThreadData &per_thread_data)
template<typename InputVector , class DH >
static void integrate_over_regular_face (const DH &dof_handler, const hp::QCollection< dim-1 > &quadrature, const typename FunctionMap< dim >::type &neumann_bc, const std::vector< const InputVector * > &solutions, const std::vector< bool > &component_mask, const Function< dim > *coefficients, typename FaceIntegrals< DH >::type &face_integrals, PerThreadData &per_thread_data, const typename DH::active_cell_iterator &cell, const unsigned int face_no, hp::FEFaceValues< dim > &fe_face_values_cell, hp::FEFaceValues< dim > &fe_face_values_neighbor)
template<typename InputVector , class DH >
static void integrate_over_irregular_face (const DH &dof_handler, const hp::QCollection< dim-1 > &quadrature, const std::vector< const InputVector * > &solutions, const std::vector< bool > &component_mask, const Function< dim > *coefficients, typename FaceIntegrals< DH >::type &face_integrals, PerThreadData &per_thread_data, const typename DH::active_cell_iterator &cell, const unsigned int face_no, hp::FEFaceValues< dim > &fe_face_values, hp::FESubfaceValues< dim > &fe_subface_values)

Friends

class PerThreadData

Detailed Description

template<int dim, int spacedim = dim>
class KellyErrorEstimator< dim, spacedim >

Implementation of the error estimator by Kelly, Gago, Zienkiewicz and Babuska. This error estimator tries to approximate the error per cell by integration of the jump of the gradient of the solution along the faces of each cell. It can be understood as a gradient recovery estimator; see the survey of Ainsworth for a complete discussion.

It seem as if this error estimator should only be valid for linear trial spaces, and there are indications that for higher order trial spaces the integrals computed here show superconvergence properties, i.e. they tend to zero faster than the error itself, thus ruling out the values as error indicators.

The error estimator really only estimates the error for the generalized Poisson equation $-\nabla\cdot a(x) \nabla u = f$ with either Dirichlet boundary conditions or generalized Neumann boundary conditions involving the conormal derivative $a\frac{du}{dn} = g$.

The error estimator returns a vector of estimated errors per cell which can be used to feed the Triangulation<dim>::refine_* functions. This vector contains elements of data type float, rather than double, since accuracy is not so important here, and since this can save rather a lot of memory, when using many cells.

Implementation

In principle, the implementation of the error estimation is simple: let

\[ \eta_K^2 = \frac h{24} \int_{\partial K} \left[a \frac{\partial u_h}{\partial n}\right]^2 do \]

be the error estimator for cell $K$. $[\cdot]$ denotes the jump of the argument at the face. In the paper of Ainsworth, $h$ is divided by $24$, but this factor is a bit esoteric, stemming from interpolation estimates and stability constants which may hold for the Poisson problem, but may not hold for more general situations. In the implementation, this factor is considered, but may lead to wrong results. You may scale the vector appropriately afterwards.

To perform the integration, use is made of the FEFaceValues and FESubfaceValues classes. The integration is performed by looping over all cells and integrating over faces that are not yet treated. This way we avoid integration on faces twice, once for each time we visit one of the adjacent cells. In a second loop over all cells, we sum up the contributions of the faces (which are the integrated square of the jumps) of each cell and take the square root.

The integration is done using a quadrature formula on the face. For linear trial functions (FEQ1), the QGauss2 or even the QMidpoint rule will suffice. For higher order elements, it is necessary to utilize higher order quadrature formulae as well.

We store the contribution of each face in a map, as provided by the C++ standard library, with the iterator pointing to that face being the key into the map. In fact, we do not store the indicator per face, but only the integral listed above. When looping the second time over all cells, we have to sum up the contributions of the faces, multiply them with $\frac h{24}$ and take the square root. By doing the multiplication with $h$ in the second loop, we avoid problems to decide with which $h$ to multiply, that of the cell on the one or that of the cell on the other side of the face.

$h$ is taken to be the greatest length of the diagonals of the cell. For more or less uniform cells without deformed angles, this coincides with the diameter of the cell.

Vector-valued functions

If the finite element field for which the error is to be estimated is vector-valued, i.e. the finite element has more than one component, then all the above can be applied to all or only some components at the same time. The main function of this class takes a list of flags denoting the components for which components the error estimator is to be applied; by default, it is a list of only trues, meaning that all variables shall be treated.

In case the different components of a field have different physical meaning (for example velocity and pressure in the Stokes equations), it would be senseless to use the same coefficient for all components. In that case, you can pass a function with as many components as there are components in the finite element field and each component of the error estimator will then be weighted by the respective component in this coefficient function. In the other case, when all components have the same meaning (for example the displacements in Lame's equations of elasticity), you can specify a scalar coefficient which will then be used for all components.

Boundary values

If the face is at the boundary, i.e. there is no neighboring cell to which the jump in the gradiend could be computed, there are two possibilities:

The object describing the boundary conditions is obtained from the triangulation.

Thanks go to Franz-Theo Suttmeier for clarifications about boundary conditions.

Handling of hanging nodes

The integration along faces with hanging nodes is quite tricky, since one of the elements has to be shifted one level up or down. See the documentation for the FESubFaceValues class for more information about technical issues regarding this topic.

In praxi, since we integrate over each face only once, we do this when we are on the coarser one of the two cells adjacent to a subface (a subface is defined to be the child of a face; seen from the coarse cell, it is a subface, while seen from the refined cell it is one of its faces). The reason is that finding neighborship information is a bit easier then, but that's all practical reasoning, nothing fundamental.

Since we integrate from the coarse side of the face, we have the mother face readily at hand and store the result of the integration over that mother face (being the sum of the integrals along the subfaces) in the abovementioned map of integrals as well. This consumes some memory more than needed, but makes the summing up of the face contributions to the cells easier, since then we have the information from all faces of all cells at hand and need not think about explicitly determining whether a face was refined or not. The same applies for boundary faces, see above.

Multiple solution vectors

In some cases, for example in time-dependent problems, one would like to compute the error estimates for several solution vectors on the same grid at once, with the same coefficients, boundary condition object, etc, e.g. for the solutions on several successive time steps. One could then call the functions of this class several times for each solution. However, the largest factor in the computation of the error estimates (in terms of computing time) is initialization of FEFaceValues and FESubFaceValues objects, and iterating through all faces and subfaces. If the solution vectors live on the same grid, this effort can be reduced significantly by treating all solution vectors at the same time, initializing the FEFaceValues objects only once per cell and for all solution vectors at once, and also only looping through the triangulation only once. For this reason, besides the estimate function in this class that takes a single input vector and returns a single output vector, there is also a function that accepts several in- and output vectors at the same time.

Author:
Wolfgang Bangerth, 1998, 1999, 2000, 2004, 2006; parallelization by Thomas Richter, 2000

Member Function Documentation

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const Mapping< dim, spacedim > &  mapping,
const DH &  dof,
const Quadrature< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const InputVector &  solution,
Vector< float > &  error,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Implementation of the error estimator described above. You may give a coefficient, but there is a default value which denotes the constant coefficient with value one. The coefficient function may either be a scalar one, in which case it is used for all components of the finite element, or a vector-valued one with as many components as there are in the finite element; in the latter case, each component is weighted by the respective component in the coefficient.

You might give a list of components you want to evaluate, in case the finite element used by the DoFHandler object is vector-valued. You then have to set those entries to true in the bit-vector component_mask for which the respective component is to be used in the error estimator. The default is to use all components, which is done by either providing a bit-vector with all-set entries, or an empty bit-vector.

The estimator supports multithreading and splits the cells to multithread_info.n_default_threads (default) threads. The number of threads to be used in multithreaded mode can be set with the last parameter of the error estimator. Multithreading is only implemented in two and three dimensions.

The subdomain_id parameter indicates whether we shall compute indicators for all cells (in case its value is the default, numbers::invalid_unsigned_int), or only for the cells belonging to a certain subdomain with the given indicator. The latter case is used for parallel computations where all processor nodes have the global grid stored, and could well compute all the indicators for all cells themselves, but where it is more efficient to have each process compute only indicators for the cells it owns, and have them exchange the resulting information afterwards. This is in particular true for the case where meshes are very large and computing indicators for every cells is too expensive, while computing indicators for only local cells is acceptable. Note that if you only ask for the indicators of a certain subdomain to be computed, you must nevertheless make sure that this function has access to the correct node values of all degrees of freedom. This is since the function needs to access DoF values on neighboring cells as well, even if they belong to a different subdomain.

The material_id parameter has a similar meaning: if not set to its default value, it means that indicators will only be computed for cells with this particular material id.

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const DH &  dof,
const Quadrature< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const InputVector &  solution,
Vector< float > &  error,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Calls the estimate function, see above, with mapping=MappingQ1<dim>().

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const Mapping< dim, spacedim > &  mapping,
const DH &  dof,
const Quadrature< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const std::vector< const InputVector * > &  solutions,
std::vector< Vector< float > * > &  errors,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Same function as above, but accepts more than one solution vector and returns one error vector for each solution vector. For the reason of existence of this function, see the general documentation of this class.

Since we do not want to force the user of this function to copy around their solution vectors, the vector of solution vectors takes pointers to the solutions, rather than being a vector of vectors. This makes it simpler to have the solution vectors somewhere in memory, rather than to have them collected somewhere special. (Note that it is not possible to construct of vector of references, so we had to use a vector of pointers.)

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const DH &  dof,
const Quadrature< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const std::vector< const InputVector * > &  solutions,
std::vector< Vector< float > * > &  errors,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Calls the estimate function, see above, with mapping=MappingQ1<dim>().

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const Mapping< dim, spacedim > &  mapping,
const DH &  dof,
const hp::QCollection< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const InputVector &  solution,
Vector< float > &  error,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Equivalent to the set of functions above, except that this one takes a quadrature collection for hp finite element dof handlers.

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const DH &  dof,
const hp::QCollection< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const InputVector &  solution,
Vector< float > &  error,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Equivalent to the set of functions above, except that this one takes a quadrature collection for hp finite element dof handlers.

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const Mapping< dim, spacedim > &  mapping,
const DH &  dof,
const hp::QCollection< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const std::vector< const InputVector * > &  solutions,
std::vector< Vector< float > * > &  errors,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Equivalent to the set of functions above, except that this one takes a quadrature collection for hp finite element dof handlers.

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate ( const DH &  dof,
const hp::QCollection< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const std::vector< const InputVector * > &  solutions,
std::vector< Vector< float > * > &  errors,
const std::vector< bool > &  component_mask = std::vector< bool >(),
const Function< dim > *  coefficients = 0,
const unsigned int  n_threads = multithread_info.n_default_threads,
const unsigned int  subdomain_id = numbers::invalid_unsigned_int,
const unsigned int  material_id = numbers::invalid_unsigned_int 
) [inline, static]

Equivalent to the set of functions above, except that this one takes a quadrature collection for hp finite element dof handlers.

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::estimate_some ( const hp::MappingCollection< dim, spacedim > &  mapping,
const DH &  dof_handler,
const hp::QCollection< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const std::vector< const InputVector * > &  solutions,
const std::pair< const std::vector< bool > *, const Function< dim > * > &  component_mask_and_coefficients,
const std::pair< unsigned int, unsigned int this_thread,
typename FaceIntegrals< DH >::type &  face_integrals,
PerThreadData per_thread_data 
) [inline, static, private]

Computates the error on all cells of the domain with the number n, satisfying n=this_thread (mod n_threads) This enumeration is chosen to generate a random distribution of all cells.

This function is only needed in two or three dimensions. The error estimator in one dimension is implemented seperatly.

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::integrate_over_regular_face ( const DH &  dof_handler,
const hp::QCollection< dim-1 > &  quadrature,
const typename FunctionMap< dim >::type &  neumann_bc,
const std::vector< const InputVector * > &  solutions,
const std::vector< bool > &  component_mask,
const Function< dim > *  coefficients,
typename FaceIntegrals< DH >::type &  face_integrals,
PerThreadData per_thread_data,
const typename DH::active_cell_iterator &  cell,
const unsigned int  face_no,
hp::FEFaceValues< dim > &  fe_face_values_cell,
hp::FEFaceValues< dim > &  fe_face_values_neighbor 
) [inline, static, private]

Actually do the computation on a face which has no hanging nodes (it is regular), i.e. either on the other side there is nirvana (face is at boundary), or the other side's refinement level is the same as that of this side, then handle the integration of these both cases together.

The meaning of the parameters becomes clear when looking at the source code. This function is only externalized from estimate_error to avoid ending up with a function of 500 lines of code.

template<int dim, int spacedim = dim>
template<typename InputVector , class DH >
static void KellyErrorEstimator< dim, spacedim >::integrate_over_irregular_face ( const DH &  dof_handler,
const hp::QCollection< dim-1 > &  quadrature,
const std::vector< const InputVector * > &  solutions,
const std::vector< bool > &  component_mask,
const Function< dim > *  coefficients,
typename FaceIntegrals< DH >::type &  face_integrals,
PerThreadData per_thread_data,
const typename DH::active_cell_iterator &  cell,
const unsigned int  face_no,
hp::FEFaceValues< dim > &  fe_face_values,
hp::FESubfaceValues< dim > &  fe_subface_values 
) [inline, static, private]

The same applies as for the function above, except that integration is over face face_no of cell, where the respective neighbor is refined, so that the integration is a bit more complex.


Friends And Related Function Documentation

template<int dim, int spacedim = dim>
friend class PerThreadData [friend]

By the resolution of Defect Report 45 to the ISO C++ 1998 standard, nested classes automatically have access to members of the enclosing class. Nevertheless, some compilers don't implement this resolution yet, so we have to make them friend, which doesn't hurt on the other compilers as well.


The documentation for this class was generated from the following file:

deal.II documentation generated on Mon Nov 23 22:57:51 2009 by doxygen 1.6.1