\def\jump#1{\left[\negthinspace\left[{#1}\right]\negthinspace\right]} \section{Implementation of error estimators}% \label{S:estimator}% \idx{error estimators|(} \subsection{Error estimator for elliptic problems}% \label{S:ellipt_est} \ALBERTA provides a residual type error estimator for non--linear elliptic problems of the type \begin{align*} -\nabla \cdot A \nabla u(x) + f\Big(x,u(x),\nabla u(x)\Big) &= 0 &&x \in \Omega,\\ u(x) &= g_d &&x \in \Gamma_D,\\ \nu\cdot A\nabla u(x) &= g_n &&x \in \Gamma_N, \end{align*} where $A \in \R^{n\times n}$ is a positive definite matrix and $\partial\Omega=\Gamma_D\cup\Gamma_N$. %% \ALBERTA implements for this kind of equations the $L^2$ and $H^1$ per-element estimators $\eta_{S,0}$ and $\eta_{S,1}$ ($S\in\tria$) \[ \begin{split} \eta_{S,0}^2&:= C_0^2\, h_S^4\, \|-\nabla \cdot A \nabla \uh + f(.,\uh,\nabla\uh)\|_{L^2(S)}^2\\ &\qquad+ C_1^2\, \sum_{\Gamma\subset\partial S \cap\Omega} h_S^3\, \|\jump{A\nabla \uh}\|_{L^2(\Gamma)}^2 + C_1^2\, \sum_{\Gamma\subset\partial S \cap\Gamma_N} h_S^3\, \|\nu\cdot A\nabla \uh-g_n\|_{L^2(\Gamma)}^2, \\ \eta_{S,1}^2&:= C_0^2\,h_S^2\, \|-\nabla \cdot A \nabla \uh + f(.,\uh,\nabla \uh)\|_{L^2(S)}^2 \\ &\qquad + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Omega} h_S\, \|\jump{A\nabla \uh}\|_{L^2(\Gamma)}^2 + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Gamma_N} h_S\, \|\nu\cdot A\nabla \uh-g_n\|_{L^2(\Gamma)}^2, \\ \end{split} \] where $\jump{.}$ denotes the jump of the normal component across an interior co-dimension $1$ sub-simplex (vertex/edge/face) $\Gamma \subset \partial S$. Verf\"urth proved for $g_d\equiv 0$ and $g_n\equiv 0$ in \cite{Verfuerth:94b} -- under suitable assumptions on $f$, $u$ and $\uh$ in the non-linear case -- the estimate \[ \|u - \uh\|_{H^1(\Omega)}^2\leq\sum_{S\in\tria} \eta_{S,1}^2, \] and B\"ansch and Siebert \cite{BaenschSiebert:95} proved a similar the $L^2$-estimate for the semi--linear case $f = f(x,u)$ and $g_d\equiv 0$ and $\Gamma_N=\emptyset$: \[ \|u - \uh\|_{L^2(\Omega)}^2\leq \sum_{S\in\tria}\eta_{S,0}^2. \] %% The following functions implement above estimators for scalar and vector-valued functions; the implementation works also for meshes with non-zero co-dimension as well as for periodic meshes. %% \fdx{ellipt_est()@{\code{ellipt\_est()}}}% \idx{error estimators!ellipt_est()@{\code{ellipt\_est()}}}% \fdx{ellipt_est_dow()@{\code{ellipt\_est\_dow()}}}% \idx{error estimators!ellipt_est_dow()@{\code{ellipt\_est\_dow()}}}% \fdx{ellipt_est_d()@{\code{ellipt\_est\_d()}}}% \idx{error estimators!ellipt_est_d()@{\code{ellipt\_est\_d()}}}% \bv\begin{verbatim} REAL ellipt_est(const DOF_REAL_VEC *uh, ADAPT_STAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), int quad_deg, NORM norm, REAL C[3], const REAL_DD A, const BNDRY_FLAGS dirichlet_bndry, REAL (*f)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D grd_uh_gp), FLAGS f_flags, REAL (*gn)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D normal), FLAGS gn_flags); REAL ellipt_est_dow(const DOF_REAL_VEC_D *uh, ADAPT_STAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), int quad_deg, NORM norm, REAL C[3], const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype, bool sym_grad, const BNDRY_FLAGS dirichlet_bndry, const REAL *(*f)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_DD grd_uh_gp), FLAGS f_flags, const REAL *(*gn)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_D normal), FLAGS gn_flags); REAL ellipt_est_d(const DOF_REAL_D_VEC *uh, ADAPT_STAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), int quad_deg, NORM norm, REAL C[3], const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype, bool sym_grad, const BNDRY_FLAGS dirichlet_bndry, const REAL *(*f)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_DD grd_uh_gp), FLAGS f_flags, const REAL *(*gn)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_D normal), FLAGS gn_flags); \end{verbatim}\ev %% Description: \begin{descr} \kitem{ellipt\_est(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,} \kitem{~~~~~~~~~~~A, dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}~\hfill computes an error estimate of the above type for the $H^1$ or $L^2$ norm; the return value is an approximation of the estimate $\|u - \uh\|$ by quadrature. \begin{descr} \kitem{uh} is a vector storing the coefficients of the discrete solution; if \code{uh} is a \nil pointer, nothing is done, the return value is \code{.0}. \kitem{adapt} is a pointer to an \code{ADAPT\_STAT} structure; if not \nil, the entries \code{adapt->p=2}, \code{err\_sum}, and \code{err\_max} of \code{adapt} are set by \code{ellipt\_est()} (compare \secref{S:adapt_stat_in_ALBERTA}). \kitem{rw\_el\_est} is a function for writing the local error indicator for a single element (usually to some location inside \code{leaf\_data}, compare \secref{S:leaf_data_info}); if this function is \nil, only the global estimate is computed, no local indicators are stored. \code{rw\_el\_est(el)} returns for each leaf element \code{el} a pointer to a \code{REAL} for storing the square of the element indicator, which can directly be used in the adaptive method, compare the \code{get\_el\_est()} function pointer in the \code{ADAPT\_STAT} structure (compare \secref{S:adapt_stat_in_ALBERTA}). \kitem{rw\_el\_estc} is a function for writing the local coarsening error indicator for a single element (usually to some location inside \code{leaf\_data}, compare \secref{S:leaf_data_info}); if this function is \nil, no coarsening error indicators are computed and stored; \code{rw\_el\_estc(el)} returns for each leaf element \code{el} a pointer to a \code{REAL} for storing the square of the element coarsening error indicator. \kitem{quad\_deg} is the degree of the quadrature that should be used for the approximation of the norms on the elements and edges/faces; if \code{degree} is less than zero a quadrature which is exact of degree \code{2*uh->fe\_space->bas\_fcts->degree} is used. \kitem{norm} can be either \code{H1\_NORM}\cdx{H1_NORM@{\code{H1\_NORM}}} or \code{L2\_NORM}\cdx{L2_NORM@{\code{L2\_NORM}}} (which are defined as symbolic constants in \code{alberta.h}) to indicate that the $H^1$ or $L^2$ error estimate has to be calculated. \kitem{C[0], C[1], C[2]} are the constants in front of the element residual, wall residual, and coarsening term respectively. If \code{C} is \nil, then all constants are set to $1.0$. \kitem{A} is the constant matrix of the second order term. \kitem{dirichlet\_bndry} A bit-mask marking those parts of the boundary which are subject to Dirichlet boundary conditions, see \secref{S:boundary}. \kitem{f} is a pointer to a function for the evaluation of the lower order terms at all quadrature nodes, i.e. $f(x(\lambda), u(\lambda), \nabla u(\lambda))$ ; if \code{f} is a \nil pointer, $f\equiv0$ is assumed; \code{f(el\_info, quad, qp, uh\_qp, grd\_uh\_qp)} returns the value of the lower oder terms on element \code{el\_info->el} at the quadrature node \code{quad->lambda[qp]}, where \code{uh\_qp} is the value and \code{grd\_uh\_qp} the gradient (with respect to the Cartesian coordinates) of the discrete solution at that quadrature point. See also \code{f\_flag} below: \kitem{f\_flag} specifies whether the function \code{f()} actually needs values of \code{uh\_qp} or \code{grd\_uh\_qp}, \code{f\_flag} may be $0$ or \code{INIT\_UH}\cdx{INIT_UH@{\code{INIT\_UH}}} or \code{INIT\_GRD\_UH}\cdx{INIT_GRD_UH@{\code{INIT\_GRD\_UH}}} or their bitwise composition (\code{|}). The arguments \code{uh\_qp} and \code{grd\_uh\_qp} of \code{f()} only hold valid information if the flags \code{INIT\_UH} respectively \code{INIT\_GRD\_UH} are set. \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to a function for the evaluation of non-homogeneous Neumann boundary data. \code{gn} may be \nil, in which case zero Neumann boundary conditions are assumed. The argument \code{normal} always contains the normal of the Neumann boundary facet. In the case of non-vanishing co-dimension \code{normal} lies in the lower-dimensional space which is spanned by the mesh simplex defined by \code{el\_info}. \code{gn()} is evaluated on those parts of the boundary which are \emph{not} flagged as Dirichlet-boundaries by the argument \code{dirichlet\_bndry}. \kitem{gn\_flag} controls whether the argument \code{uh\_qp} of the function \code{gn()} actually contains the value of \code{uh} at the quadrature point \code{qp}. Note that the argument \code{normal} always contains valid data. \end{descr} The estimate is computed by traversing all leaf elements of \code{uh->fe\_space->mesh}, using the quadrature for the approximation of the residuals and storing the square of the element indicators on the elements (if \code{rw\_el\_est} and \code{rw\_el\_estc }are not \nil). \kitem{ellipt\_est\_d(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,} \kitem{~~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} \kitem{~~~~~~~~~~~~~dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}\hfill \kitem{ellipt\_est\_dow(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,} \kitem{~~~~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} \kitem{~~~~~~~~~~~~~~~dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}\hfill Similar function for a (coupled) vector valued elliptic problem. We document only the arguments which are different from the arguments of \code{ellipt\_est()}: \begin{descr} \kitem{A} now represents a tensor $(A_{ij}^{\mu\nu}\in\R^{n\times n,n\times n}$, $i,j,\mu,\nu=0,\dots,n-1$. The indexing is \begin{equation*} \code{A[i][j][mu][nu]} = A^{\code{mu},\code{nu}}_{\code{ij}}, \end{equation*} with \code{i,j,mu,nu==0,\dots,\code{DIM\_OF\_WORLD-1}}, see \secref{book:S:DisCoupled}. \code{A} describes the coefficients of the principal part of a coupled system of elliptical equations: \[ -\sum_{\nu,i,j=0}^{n-1}\partial_i A_{ij}^{\mu\nu}\partial_j u^\nu +\text{lower order terms} = f^\mu\quad(\mu = 0,\,\dots,\,n-1). \] The \code{quasi-stokes.c} demo-program contains an example. \kitem{A\_blocktype} must be one of \code{MATENT\_REAL}, \code{MATENT\_REAL\_D} or \code{MATENT\_REAL\_DD}. It specifies the symmetry type for coupling of the PDE system. Note that the storage layout of \code{A} is determined by the argument \code{A\_blocktype}: \begin{descr} \kitem{MATENT\_REAL:~~~} \code{REAL A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];} \kitem{MATENT\_REAL\_D:~} \code{REAL\_D A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];} \kitem{MATENT\_REAL\_DD:} \code{REAL\_DD A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];} \end{descr} \code{A\_blocktype == MATENT\_REAL} or \code{A\_blocktype == MATENT\_REAL\_D} means that the system is actually decoupled. \kitem{A\_type} must be one of \code{MATENT\_REAL}, \code{MATENT\_REAL\_D} or \code{MATENT\_REAL\_DD}. It specifies the symmetry type of \code{A} with respect to the first two indices. For a Laplacian, for example, one would use \code{DOWBM\_SCAL}. Note that the value of \code{A\_type} does \emph{not} change the storage layout of the array \code{A}. \kitem{sym\_grad} If set to \code{true} then it is assumed that the symmetric gradient has to be used for the computation of the jump- and Neumann-residuals. The demo-program \code{quasi-stokes.c} uses this feature to implement an error estimator for the Stokes equation with stress boundary conditions. \kitem{f} If the first argument of the function pointer \code{f(result,\dots)} is not \nil then the result \emph{must} be stored in the argument \code{result} and \code{f()} must return the base address of the array \code{result}. If \code{result} is \nil, then \code{f()} must store the result in a non-volatile storage area and return the address of that area. \kitem{dirichlet\_bndry} A bit-mask marking those parts of the boundary which are subject to Dirichlet boundary conditions, see \secref{S:boundary}. \kitem{f} is a pointer to a function for the evaluation of the lower order terms at all quadrature nodes, i.e. $f(x(\lambda), u(\lambda), \nabla u(\lambda))$ ; if \code{f} is a \nil pointer, $f\equiv0$ is assumed; \code{f(el\_info, quad, qp, uh\_qp, grd\_uh\_qp)} returns the value of the lower oder terms on element \code{el\_info->el} at the quadrature node \code{quad->lambda[qp]}, where \code{uh\_qp} is the value and \code{grd\_uh\_qp} the gradient (with respect to the Cartesian coordinates) of the discrete solution at that quadrature point. See also \code{f\_flag} below: \kitem{f\_flag} specifies whether the function \code{f()} actually needs values of \code{uh\_qp} or \code{grd\_uh\_qp}, \code{f\_flag} may be $0$ or \code{INIT\_UH}\cdx{INIT_UH@{\code{INIT\_UH}}} or \code{INIT\_GRD\_UH}\cdx{INIT_GRD_UH@{\code{INIT\_GRD\_UH}}} or their bitwise composition (\code{|}). The arguments \code{uh\_qp} and \code{grd\_uh\_qp} of \code{f()} only hold valid information if the flags \code{INIT\_UH} respectively \code{INIT\_GRD\_UH} are set. \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to a function for the evaluation of non-homogeneous Neumann boundary data. \code{gn} may be \nil, in which case zero Neumann boundary conditions are assumed. The argument \code{normal} always contains the normal of the Neumann boundary facet. In the case of non-vanishing co-dimension \code{normal} lies in the lower-dimensional space which is spanned by the mesh simplex defined by \code{el\_info}. \code{gn()} is evaluated on those parts of the boundary which are \emph{not} flagged as Dirichlet-boundaries by the argument \code{dirichlet\_bndry}. \kitem{gn\_flag} controls whether the argument \code{uh\_qp} of the function \code{gn()} actually contains the value of \code{uh} at the quadrature point \code{qp}. Note that the argument \code{normal} always contains valid data. \end{descr} \end{descr} \begin{example}[Linear problem]\label{E:est-impl} Consider the scalar linear model problem \mathref{book:E:strong} with constant coefficients $A$, $b$, and $c$: \begin{alignat*}{2} -\nabla \cdot A \nabla u + b \cdot \nabla u + c\, u &= r \qquad & &\mbox{in } \Omega,\\ u &= 0 & &\mbox{on } \partial\Omega. \end{alignat*} Let \code{A} be a \code{REAL\_DD} matrix storing $A$, which is then the eighth argument of \code{ellipt\_est()}. Assume that \code{const REAL *b(const REAL\_D)} is a function returning a pointer to a vector storing $b$, \code{REAL c(REAL\_D)} returns the value of $c$ and \code{REAL r(const REAL\_D)} returns the value of the right hand side $r$ of \mathref{book:E:strong} at some point in world coordinates. The implementation of the function \code{f} is: \bv\begin{verbatim} static REAL f(const EL_INFO *el_info, const QUAD *quad, int iq, REAL uh_iq, const REAL_D grd_uh_iq) { FUNCNAME("f"); const REAL *bx, *x; extern const REAL b(const REAL_D); extern REAL c(const REAL_D), r(const REAL_D); x = coord_to_world(el_info, quad->lambda[iq], nil); bx = b(x); return(SCP_DOW(bx, grd_uh_iq) + c(x)*uh_iq - r(x)); } \end{verbatim}\ev As both \code{uh\_iq} and \code{grd\_uh\_iq} are used, the estimator parameter \code{f\_flag} must be given as \code{INIT\_UH|INIT\_GRD\_UH}. \end{example} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Error estimator for parabolic problems}\label{S:para_est} Similar to the stationary case, the \ALBERTA library provides an error estimator for the non--linear parabolic problem \begin{align*} \partial_t u -\nabla \cdot A \nabla u(x) + f\Big(x, t, u(x),\nabla u(x)\Big) &= 0 &&x \in \Omega, t>0,\\ u(x,t) &= g_d && x \in \Gamma_D, t>0,\\ \nu\cdot A\nabla u(x,t) &= g_n && x \in \Gamma_N, t>0,\\ u(x,0) &= u_0 && x \in \Omega, \end{align*} where $A \in \R^{d\times d}$ is a positive definite matrix and $\partial\Omega=\Gamma_D\cup\Gamma_N$. The estimator is split in several parts, where the initial error \[ \eta_0 = \|u_0 - U_0\|_{L^2(\Omega)} \] can be approximated by the function \code{L2\_err()}, e.g. (compare \secref{S:error_cal}). For the estimation of the spatial discretization error, the coarsening error, and the time discretization error, the \ALBERTA estimator implements the following (local) indicators \begin{align*} \eta_S^2 &= C_0^2\, h_S^4\, \left\| \frac{U_{n+1} - I_{n+1}U_{n}}{\tau_{n+1}} -\nabla \cdot A \nabla U_{n+1} + f(.,t_{n+1},U_{n+1},\nabla U_{n+1}) \right\|_{L^2(S)}^2\\ &\qquad + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Omega} h_S^3\, \|\jump{A\nabla U_{n+1}}\|_{L^2(\Gamma)}^2 + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Gamma_N} h_S^3\, \|\nu\cdot A\nabla U_{n+1}-g_n\|_{L^2(\Gamma)}^2, \\[2mm] \eta_{S,c}^2 &= C_2^2\, h_S^3\, \|\jump{\nabla U_{n}}\|_{L^2(\Gamma_c)}^2 % \|U_{n} - I_{n+1}U_{n}\|_{L^2(S)}^2 \\[2mm] \eta_\tau &= C_3 \|U_{n+1} - I_{n+1}U_{n}\|_{L^2(\Omega)}. \end{align*} The coarsening indicator is motivated by the fact that for piecewise linear Lagrange finite element functions it holds $\|U_{n} - I_{n+1}U_{n}\|_{L^2(S)}^2 = \eta_{S,c}^2$ with $C_2=C_2(d)$ and $\Gamma_c$ the face that would be removed during a coarsening operation. %% The implementation is done by the functions \fdx{heat_est()@{\code{heat\_est()}}}% \fdx{heat_est_dow()@{\code{heat\_est\_dow()}}}% \fdx{heat_est_d()@{\code{heat\_est\_d()}}}% \idx{error estimators!heat_est()@{\code{heat\_est()}}}% \idx{error estimators!heat_est_dow()@{\code{heat\_est\_dow()}}}% \idx{error estimators!heat_est_d()@{\code{heat\_est\_d()}}}% \bv\begin{verbatim} REAL heat_est(const DOF_REAL_VEC *uh, ADAPT_INSTAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), int quad_degree, REAL C[4], const DOF_REAL_VEC *uh_old, const REAL_DD A, const BNDRY_FLAGS dirichlet_bndry, REAL (*f)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D grd_uh_gp, REAL time), FLAGS f_flags, REAL (*gn)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D normal, REAL time), FLAGS gn_flags); REAL heat_est_dow(const DOF_REAL_D_VEC *uh, ADAPT_INSTAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), int quad_degree, REAL C[4], const DOF_REAL_D_VEC *uh_old, const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype, bool sym_grad, BNDRY_FLAGS dirichlet_bndry, const REAL *(*f)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_DD grd_uh_gp, REAL time), FLAGS f_flags, const REAL *(*gn)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_D normal, REAL time), FLAGS gn_flags); REAL heat_est_d(const DOF_REAL_D_VEC *uh, const DOF_REAL_D_VEC *uh_old, ADAPT_INSTAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), int quad_degree, REAL C[4], const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype, bool sym_grad, const BNDRY_FLAGS dirichlet_bndry, const REAL *(*f)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_DD grd_uh_gp, REAL time), FLAGS f_flags, const REAL *(*gn)(REAL_D result, const EL_INFO *el_info, const QUAD *quad, int qp, const REAL_D uh_qp, const REAL_D normal, REAL time), FLAGS gn_flags); \end{verbatim}\ev Description: \begin{descr} \kitem{heat\_est(uh, adapt, rw\_el\_est, rw\_el\_estc, degree, C, uh\_old,} \kitem{~~~~~~~~~A, dirichlet\_bndry, f, f\_flag, gn, gn\_flag)}\hfill computes an error estimate of the above type, the local and global space discretization estimators are stored in \code{adapt->adapt\_space} and via the \code{rw\_...} pointers; the return value is the time discretization indicator $\eta_\tau$. \begin{descr} \kitem{uh} is a vector storing the coefficients of the discrete solution $U_{n+1}$; if \code{uh} is a \nil pointer, nothing is done, the return value is \code{0.0}. \kitem{adapt} is a pointer to an \code{ADAPT\_INSTAT} structure; if it is not \nil, then the entries \code{adapt\_space->p=2}, \code{adapt\_space->err\_sum} and \code{adapt\_space->err\_max} of \code{adapt} are set by \code{heat\_est()} (compare \secref{S:adapt_stat_in_ALBERTA}). \kitem{rw\_el\_est} is a function for writing the local error indicator $\eta_S^2$ for a single element (usually to some location inside \code{leaf\_data}, compare \secref{S:leaf_data_info}); if this function is \nil, only the global estimate is computed, no local indicators are stored. \code{rw\_el\_est(el)} returns for each leaf element \code{el} a pointer to a \code{REAL} for storing the square of the element indicator, which can directly be used in the adaptive method, compare the \code{get\_el\_est()} function pointer in the \code{ADAPT\_STAT} structure (compare \secref{S:adapt_stat_in_ALBERTA}). \kitem{rw\_el\_estc} is a function for writing the local coarsening error indicator $\eta_{S,c}^2$ for a single element (usually to some location inside \code{leaf\_data}, compare \secref{S:leaf_data_info}); if this function is \nil, no coarsening error indicators are computed and stored; \code{rw\_el\_estc(el)} returns for each leaf element \code{el} a pointer to a \code{REAL} for storing the square of the element coarsening error indicator. The coarsening indicator is not used at the moment. \kitem{degree} is the degree of the quadrature that should be used for the approximation of the norms on the elements and edges/faces; if \code{degree} is less than zero a quadrature which is exact of degree \code{2*uh->fe\_space->bas\_fcts->degree} is used. \kitem{C[0]}, \code{C[1]}, \code{C[2]}, \code{C[3]} are the constants in front of the element residual, wall residual, coarsening term, and time residual, respectively. If \code{C} is \nil, then all constants are set to $1.0$. \kitem{uh\_old} is a vector storing the coefficients of the discrete solution $U_{n}$ from previous time step; if \code{uh\_old} is a \nil pointer, nothing is done, the return value is \code{0.0}. \kitem{A} is the constant matrix of the second order term. \kitem{dirichlet\_bndry} A bit mask marking those parts of the boundary which are subject to Dirichlet boundary conditions. See \secref{S:boundary}. \kitem{f} is a pointer to a function for the evaluation of the lower order terms at all quadrature nodes, i.e. $f(x(\lambda), t, u(\lambda), \nabla u(\lambda))$ ; if \code{f} is a \nil pointer, $f\equiv0$ is assumed; \code{f(el\_info, quad, iq, t, uh\_iq, grd\_uh\_iq)} returns the value of the lower oder terms on element \code{el\_info->el} at the quadrature node \code{quad->lambda[iq]}, where \code{uh\_iq} is the value and \code{grd\_uh\_iq} the gradient (with respect to the world coordinates) of the discrete solution at that quadrature node. \kitem{f\_flag} specifies whether the function \code{f()} actually needs values of \code{uh\_iq} or \code{grd\_uh\_iq}. This flag may hold zero, the predefined values \code{INIT\_UH} or \code{INIT\_GRD\_UH}, or their composition \code{INIT\_UH|INIT\_GRD\_UH}; the arguments \code{uh\_iq} and \code{grd\_uh\_iq} of \code{f()} only hold valid information, if the flags \code{INIT\_UH} respectively \code{INIT\_GRD\_UH} are set. \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to a function for the evaluation of non-homogeneous Neumann boundary data. \code{gn} may be \nil, in which case zero Neumann boundary conditions are assumed. The argument \code{normal} always contains the normal of the Neumann boundary facet. In the case of non-vanishing co-dimension \code{normal} lies in the lower-dimensional space which is spanned by the mesh simplex defined by \code{el\_info}. \kitem{gn\_flag} controls whether the argument \code{uh\_qp} of the function \code{gn()} actually contains the value of \code{uh} at the quadrature point \code{qp}. Note that the argument \code{normal} always contains valid data. \end{descr} \smallskip The estimate is computed by traversing all leaf elements of \code{uh->fe\_space->mesh}, using the quadrature for the approximation of the residuals and storing the square of the element indicators on the elements (if \code{rw\_el\_est} and \code{rw\_el\_estc }are not \nil). \kitem{heat\_est\_d(uh, adapt, rw, rwc, deg, C, uh\_old,} \kitem{~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} \kitem{~~~~~~~~~~dirichlet\_bndry, f, f\_flag)}\hfill \kitem{heat\_est\_dow(uh, adapt, rw, rwc, deg, C, uh\_old,} \kitem{~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} \kitem{~~~~~~~~~~~~dirichlet\_bndry, f, f\_flag)}\hfill Coupled vector valued version. See \code{ellipt\_est\_dow()} above. \end{descr} There are also some less high-level support functions which allow for custom contributions to the per-element error estimates. We will not document this in detail, but rather refer the reader to the \code{stokes.c} and \code{quasi-stokes.c} demo-programs. %% \fdx{ellipt_est_init()@{\code{ellipt\_est\_init()}}} \fdx{heat_est_init()@{\code{heat\_est\_init()}}} \fdx{element_est()@{\code{element\_est()}}} \fdx{element_est_finish()@{\code{element\_est\_finish()}}} \fdx{element_est_uh()@{\code{element\_est\_uh()}}} \fdx{element_est_grd_uh()@{\code{element\_est\_grd\_uh()}}} \fdx{ellipt_est_finish()@{\code{ellipt\_est\_finish()}}} \fdx{heat_est_finish()@{\code{heat\_est\_finish()}}} %% \fdx{ellipt_est_dow_init()@{\code{ellipt\_est\_dow\_init()}}} \fdx{heat_est_dow_init()@{\code{heat\_est\_dow\_init()}}} \fdx{element_est_dow()@{\code{element\_est\_dow()}}} \fdx{element_est_dow_finish()@{\code{element\_est\_dow\_finish()}}} \fdx{element_est_uh_dow()@{\code{element\_est\_uh\_dow()}}} \fdx{element_est_grd_uh_dow()@{\code{element\_est\_grd\_uh\_dow()}}} \fdx{ellipt_est_dow_finish()@{\code{ellipt\_est\_dow\_finish()}}} \fdx{heat_est_dow_finish()@{\code{heat\_est\_dow\_finish()}}} %% \bv\begin{lstlisting} const void *ellipt_est_init(const DOF_REAL_VEC *uh, ADAPT_STAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), const QUAD *quad, const WALL_QUAD *wall_quad, NORM norm, REAL C[3], const REAL_DD A, const BNDRY_FLAGS dirichlet_bndry, REAL (*f)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D grd_uh_gp), FLAGS f_flags, REAL (*gn)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D normal), FLAGS gn_flags); const void *heat_est_init(const DOF_REAL_VEC *uh, const DOF_REAL_VEC *uh_old, ADAPT_INSTAT *adapt, REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), const QUAD *quad, const WALL_QUAD *wall_quad, REAL C[4], const REAL_DD A, const BNDRY_FLAGS dirichlet_bndry, REAL (*f)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D grd_uh_gp, REAL time), FLAGS f_flags, REAL (*gn)(const EL_INFO *el_info, const QUAD *quad, int qp, REAL uh_qp, const REAL_D normal, REAL time), FLAGS gn_flags); REAL element_est(const EL_INFO *el_info, const void *est_handle); void element_est_finish(const EL_INFO *el_info, REAL est_el, const void *est_handle); const REAL *element_est_uh(const void *est_handle); const REAL_D *element_est_grd_uh(const void *est_handle); REAL ellipt_est_finish(ADAPT_STAT *adapt, const void *est_handle); REAL heat_est_finish(ADAPT_INSTAT *adapt, const void *est_handle); \end{lstlisting}\ev %% There are similar proto-types for the vector-valued case. Now, what are these functions good for? The \code{stokes.c} program makes use of this framework to add a contribution concerning the divergence constraint. Of course, this is an ad-hoc error indicator, and only meant to demonstrate the programming frame-work. The functions \code{element\_est\_uh[\_dow]()} and \code{element\_est\_grd\_uh[\_dow]()} give the application access to the values of the discrete solution at the quadrature points (respectively to its Jaocbians). Otherwise, the general layout is like follows: %% \bv\begin{lstlisting} void *est_handle = ellipt_est_init(...); TRAVERSE_FIRST(mesh, -1, ) { REAL est_el = element_est(el_info, est_handle); ... /* add whatever you like to est_el */ element_est_finish(el_info, est_el, est_handle); } TRAVERSE_NEXT(); REAL est = ellipt_est_finish(adapt, est_handle); \end{lstlisting}\ev %% The relevant excerpt from \code{stokes.c} reads as follows: %% \bv\begin{lstlisting} est_handle = ellipt_est_dow_init(u_h, adapt, rw_el_est, NULL /* rw_estc */, quad, NULL /* wall_quad */, H1_NORM, C, A, MATENT_REAL, MATENT_REAL, false /* !sym_grad */, dirichlet_mask, r, INIT_GRD_UH, NULL /* inhomog. Neumann res. */, 0); fill_flags = FILL_NEIGH|FILL_COORDS|FILL_OPP_COORDS|FILL_BOUND|CALL_LEAF_EL; fill_flags |= u_fe_space->bas_fcts->fill_flags; fill_flags |= p_fe_space->bas_fcts->fill_flags; TRAVERSE_FIRST(mesh, -1, fill_flags) { const EL_GEOM_CACHE *elgc; const QUAD_EL_CACHE *qelc; REAL est_el; est_el = element_est_dow(el_info, est_handle); if (C[3]) { REAL div_uh_el, div_uh_qp; const REAL_DD *grd_uh_qp; int qp, i; grd_uh_qp = element_est_grd_uh_d(est_handle); div_uh_el = 0.0; if (!(el_info->fill_flag & FILL_COORDS)) { qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_DET); for (qp = 0; qp < quad->n_points; qp++) { div_uh_qp = 0; for (i = 0; i < DIM_OF_WORLD; i++) { div_uh_qp += grd_uh_qp[qp][i][i]; } div_uh_el += qelc->param.det[qp]*quad->w[qp]*SQR(div_uh_qp); } } else { elgc = fill_el_geom_cache(el_info, FILL_EL_DET); for (qp = 0; qp < quad->n_points; qp++) { div_uh_qp = 0; for (i = 0; i < DIM_OF_WORLD; i++) { div_uh_qp += grd_uh_qp[qp][i][i]; } div_uh_el += quad->w[qp]*SQR(div_uh_qp); } div_uh_el *= elgc->det; } est_el += C[3] * div_uh_el; } element_est_dow_finish(el_info, est_el, est_handle); } TRAVERSE_NEXT(); est = ellipt_est_dow_finish(adapt, est_handle); \end{lstlisting}\ev \idx{error estimators|)} %%% Local Variables: %%% mode: latex %%% TeX-master: "alberta-man" %%% End: