1\def\jump#1{\left[\negthinspace\left[{#1}\right]\negthinspace\right]} 2\section{Implementation of error estimators}% 3\label{S:estimator}% 4\idx{error estimators|(} 5 6\subsection{Error estimator for elliptic problems}% 7\label{S:ellipt_est} 8 9\ALBERTA provides a residual type error estimator for non--linear 10elliptic problems of the type 11\begin{align*} 12 -\nabla \cdot A \nabla u(x) + f\Big(x,u(x),\nabla u(x)\Big) &= 0 13&&x \in \Omega,\\ 14 u(x) &= g_d &&x \in \Gamma_D,\\ 15 \nu\cdot A\nabla u(x) &= g_n &&x \in \Gamma_N, 16\end{align*} 17where $A \in \R^{n\times n}$ is a positive definite matrix and 18$\partial\Omega=\Gamma_D\cup\Gamma_N$. 19%% 20\ALBERTA implements for this kind of equations the $L^2$ and $H^1$ 21per-element estimators $\eta_{S,0}$ and $\eta_{S,1}$ ($S\in\tria$) 22\[ 23\begin{split} 24 \eta_{S,0}^2&:= 25 C_0^2\, h_S^4\, \|-\nabla \cdot A \nabla \uh + f(.,\uh,\nabla\uh)\|_{L^2(S)}^2\\ 26 &\qquad+ C_1^2\, \sum_{\Gamma\subset\partial S \cap\Omega} 27 h_S^3\, \|\jump{A\nabla \uh}\|_{L^2(\Gamma)}^2 28 + C_1^2\, \sum_{\Gamma\subset\partial S \cap\Gamma_N} 29 h_S^3\, \|\nu\cdot A\nabla \uh-g_n\|_{L^2(\Gamma)}^2, 30 \\ 31 \eta_{S,1}^2&:= 32 C_0^2\,h_S^2\, \|-\nabla \cdot A \nabla \uh + f(.,\uh,\nabla \uh)\|_{L^2(S)}^2 33 \\ 34 &\qquad + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Omega} 35 h_S\, \|\jump{A\nabla \uh}\|_{L^2(\Gamma)}^2 36 + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Gamma_N} 37 h_S\, \|\nu\cdot A\nabla \uh-g_n\|_{L^2(\Gamma)}^2, 38 \\ 39\end{split} 40\] 41where $\jump{.}$ denotes the jump of the normal component across an 42interior co-dimension $1$ sub-simplex (vertex/edge/face) $\Gamma 43\subset \partial S$. 44 45Verf\"urth proved for $g_d\equiv 0$ and $g_n\equiv 0$ in 46\cite{Verfuerth:94b} -- under suitable assumptions on $f$, 47$u$ and $\uh$ in the non-linear case -- the estimate 48\[ 49\|u - \uh\|_{H^1(\Omega)}^2\leq\sum_{S\in\tria} \eta_{S,1}^2, 50\] 51and B\"ansch and Siebert \cite{BaenschSiebert:95} proved a similar the 52$L^2$-estimate for the semi--linear case $f = f(x,u)$ and $g_d\equiv 530$ and $\Gamma_N=\emptyset$: 54\[ 55\|u - \uh\|_{L^2(\Omega)}^2\leq \sum_{S\in\tria}\eta_{S,0}^2. 56\] 57%% 58The following functions implement above estimators for scalar and 59vector-valued functions; the implementation works also for meshes with 60non-zero co-dimension as well as for periodic meshes. 61%% 62\fdx{ellipt_est()@{\code{ellipt\_est()}}}% 63\idx{error estimators!ellipt_est()@{\code{ellipt\_est()}}}% 64\fdx{ellipt_est_dow()@{\code{ellipt\_est\_dow()}}}% 65\idx{error estimators!ellipt_est_dow()@{\code{ellipt\_est\_dow()}}}% 66\fdx{ellipt_est_d()@{\code{ellipt\_est\_d()}}}% 67\idx{error estimators!ellipt_est_d()@{\code{ellipt\_est\_d()}}}% 68\bv\begin{verbatim} 69REAL ellipt_est(const DOF_REAL_VEC *uh, ADAPT_STAT *adapt, 70 REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), 71 int quad_deg, 72 NORM norm, REAL C[3], const REAL_DD A, 73 const BNDRY_FLAGS dirichlet_bndry, 74 REAL (*f)(const EL_INFO *el_info, 75 const QUAD *quad, int qp, 76 REAL uh_qp, const REAL_D grd_uh_gp), 77 FLAGS f_flags, 78 REAL (*gn)(const EL_INFO *el_info, 79 const QUAD *quad, int qp, 80 REAL uh_qp, const REAL_D normal), 81 FLAGS gn_flags); 82 83REAL ellipt_est_dow(const DOF_REAL_VEC_D *uh, ADAPT_STAT *adapt, 84 REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), 85 int quad_deg, 86 NORM norm, REAL C[3], 87 const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype, 88 bool sym_grad, 89 const BNDRY_FLAGS dirichlet_bndry, 90 const REAL *(*f)(REAL_D result, 91 const EL_INFO *el_info, 92 const QUAD *quad, int qp, 93 const REAL_D uh_qp, 94 const REAL_DD grd_uh_gp), 95 FLAGS f_flags, 96 const REAL *(*gn)(REAL_D result, 97 const EL_INFO *el_info, 98 const QUAD *quad, int qp, 99 const REAL_D uh_qp, 100 const REAL_D normal), 101 FLAGS gn_flags); 102 103REAL ellipt_est_d(const DOF_REAL_D_VEC *uh, ADAPT_STAT *adapt, 104 REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), 105 int quad_deg, 106 NORM norm, REAL C[3], 107 const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype, 108 bool sym_grad, 109 const BNDRY_FLAGS dirichlet_bndry, 110 const REAL *(*f)(REAL_D result, 111 const EL_INFO *el_info, 112 const QUAD *quad, int qp, 113 const REAL_D uh_qp, 114 const REAL_DD grd_uh_gp), 115 FLAGS f_flags, 116 const REAL *(*gn)(REAL_D result, 117 const EL_INFO *el_info, 118 const QUAD *quad, int qp, 119 const REAL_D uh_qp, 120 const REAL_D normal), 121 FLAGS gn_flags); 122\end{verbatim}\ev 123%% 124Description: 125\begin{descr} 126\kitem{ellipt\_est(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,} 127\kitem{~~~~~~~~~~~A, dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}~\hfill 128 129 computes an error estimate of the above type for the 130 $H^1$ or $L^2$ norm; the return value is an 131 approximation of the estimate $\|u - \uh\|$ by quadrature. 132 133 \begin{descr} 134 135 \kitem{uh} is a vector storing the coefficients of the 136 discrete solution; if \code{uh} is a \nil pointer, nothing is 137 done, the return value is \code{.0}. 138 139 \kitem{adapt} is a pointer to an \code{ADAPT\_STAT} structure; 140 if not \nil, the entries \code{adapt->p=2}, \code{err\_sum}, 141 and \code{err\_max} of \code{adapt} are set by 142 \code{ellipt\_est()} (compare 143 \secref{S:adapt_stat_in_ALBERTA}). 144 145 \kitem{rw\_el\_est} is a function for writing the local error 146 indicator for a single element (usually to some location 147 inside \code{leaf\_data}, compare \secref{S:leaf_data_info}); 148 if this function is \nil, only the global estimate is 149 computed, no local indicators are stored. 150 \code{rw\_el\_est(el)} returns for each leaf element 151 \code{el} a pointer to a \code{REAL} for storing the square 152 of the element indicator, which can directly be used in the 153 adaptive method, compare the \code{get\_el\_est()} function 154 pointer in the \code{ADAPT\_STAT} structure (compare 155 \secref{S:adapt_stat_in_ALBERTA}). 156 157 \kitem{rw\_el\_estc} is a function for writing the local 158 coarsening error indicator for a single element (usually to 159 some location inside \code{leaf\_data}, compare 160 \secref{S:leaf_data_info}); if this function is \nil, no 161 coarsening error indicators are computed and stored; 162 \code{rw\_el\_estc(el)} returns for each leaf element 163 \code{el} a pointer to a \code{REAL} for storing the square 164 of the element coarsening error indicator. 165 166 \kitem{quad\_deg} is the degree of the quadrature that should 167 be used for the approximation of the norms on the elements 168 and edges/faces; if \code{degree} is less than zero a 169 quadrature which is exact of degree 170 \code{2*uh->fe\_space->bas\_fcts->degree} is used. 171 172 \kitem{norm} can be either 173 \code{H1\_NORM}\cdx{H1_NORM@{\code{H1\_NORM}}} or 174 \code{L2\_NORM}\cdx{L2_NORM@{\code{L2\_NORM}}} (which are 175 defined as symbolic constants in \code{alberta.h}) to 176 indicate that the $H^1$ or $L^2$ error estimate has to be 177 calculated. 178 179 \kitem{C[0], C[1], C[2]} are the constants in 180 front of the element residual, wall residual, and coarsening 181 term respectively. If \code{C} is \nil, then all constants 182 are set to $1.0$. 183 184 \kitem{A} is the constant matrix of the second order term. 185 186 \kitem{dirichlet\_bndry} A bit-mask marking those parts of the 187 boundary which are subject to Dirichlet boundary conditions, see 188 \secref{S:boundary}. 189 190 \kitem{f} is a pointer to a function for the evaluation of the 191 lower order terms at all quadrature nodes, i.e. 192 $f(x(\lambda), u(\lambda), \nabla u(\lambda))$ ; if \code{f} 193 is a \nil pointer, $f\equiv0$ is assumed; 194 195 \code{f(el\_info, quad, qp, uh\_qp, grd\_uh\_qp)} returns the 196 value of the lower oder terms on element \code{el\_info->el} 197 at the quadrature node \code{quad->lambda[qp]}, where 198 \code{uh\_qp} is the value and \code{grd\_uh\_qp} the 199 gradient (with respect to the Cartesian coordinates) of the 200 discrete solution at that quadrature point. See also 201 \code{f\_flag} below: 202 203 \kitem{f\_flag} specifies whether the function \code{f()} 204 actually needs values of \code{uh\_qp} or \code{grd\_uh\_qp}, 205 \code{f\_flag} may be $0$ or 206 \code{INIT\_UH}\cdx{INIT_UH@{\code{INIT\_UH}}} or 207 \code{INIT\_GRD\_UH}\cdx{INIT_GRD_UH@{\code{INIT\_GRD\_UH}}} 208 or their bitwise composition (\code{|}). The arguments 209 \code{uh\_qp} and \code{grd\_uh\_qp} of \code{f()} only hold 210 valid information if the flags \code{INIT\_UH} respectively 211 \code{INIT\_GRD\_UH} are set. 212 213 \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to a 214 function for the evaluation of non-homogeneous Neumann boundary 215 data. \code{gn} may be \nil, in which case zero Neumann boundary 216 conditions are assumed. The argument \code{normal} always 217 contains the normal of the Neumann boundary facet. In the case of 218 non-vanishing co-dimension \code{normal} lies in the 219 lower-dimensional space which is spanned by the mesh simplex 220 defined by \code{el\_info}. \code{gn()} is evaluated on those 221 parts of the boundary which are \emph{not} flagged as 222 Dirichlet-boundaries by the argument \code{dirichlet\_bndry}. 223 224 \kitem{gn\_flag} controls whether the argument \code{uh\_qp} 225 of the function \code{gn()} actually contains the value of 226 \code{uh} at the quadrature point \code{qp}. Note that the 227 argument \code{normal} always contains valid data. 228 \end{descr} 229 230 The estimate is computed by traversing all leaf elements of 231 \code{uh->fe\_space->mesh}, using the quadrature for the 232 approximation of the residuals and storing the square of the 233 element indicators on the elements (if \code{rw\_el\_est} and 234 \code{rw\_el\_estc }are not \nil). 235 236\kitem{ellipt\_est\_d(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,} 237\kitem{~~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} 238\kitem{~~~~~~~~~~~~~dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}\hfill 239\kitem{ellipt\_est\_dow(uh, adapt, rw\_est, rw\_estc, quad\_deg, norm, C,} 240\kitem{~~~~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} 241\kitem{~~~~~~~~~~~~~~~dirichlet\_bndry, f, f\_flags, gn, gn\_flags)}\hfill 242 243Similar function for a (coupled) vector valued elliptic problem. We 244document only the arguments which are different from the arguments of 245\code{ellipt\_est()}: 246 247\begin{descr} 248 \kitem{A} 249 now represents a tensor 250 $(A_{ij}^{\mu\nu}\in\R^{n\times n,n\times n}$, 251 $i,j,\mu,\nu=0,\dots,n-1$. The indexing is 252 \begin{equation*} 253 \code{A[i][j][mu][nu]} = A^{\code{mu},\code{nu}}_{\code{ij}}, 254 \end{equation*} 255 with \code{i,j,mu,nu==0,\dots,\code{DIM\_OF\_WORLD-1}}, see 256 \secref{book:S:DisCoupled}. \code{A} describes the coefficients of the 257 principal part of a coupled system of elliptical equations: 258 \[ 259 -\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). 260 \] 261 The \code{quasi-stokes.c} demo-program contains an example. 262 263\kitem{A\_blocktype} must be one of \code{MATENT\_REAL}, 264 \code{MATENT\_REAL\_D} or \code{MATENT\_REAL\_DD}. It specifies the 265 symmetry type for coupling of the PDE system. Note that the storage 266 layout of \code{A} is determined by the argument 267 \code{A\_blocktype}: 268 \begin{descr} 269 \kitem{MATENT\_REAL:~~~} \code{REAL A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];} 270 \kitem{MATENT\_REAL\_D:~} \code{REAL\_D A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];} 271 \kitem{MATENT\_REAL\_DD:} \code{REAL\_DD A[DIM\_OF\_WORLD][DIM\_OF\_WORLD];} 272 \end{descr} 273 \code{A\_blocktype == MATENT\_REAL} or 274 \code{A\_blocktype == MATENT\_REAL\_D} means that the system is 275 actually decoupled. 276 277\kitem{A\_type} must be one of \code{MATENT\_REAL}, 278 \code{MATENT\_REAL\_D} or \code{MATENT\_REAL\_DD}. It specifies the 279 symmetry type of \code{A} with respect to the first two indices. For 280 a Laplacian, for example, one would use \code{DOWBM\_SCAL}. Note 281 that the value of \code{A\_type} does \emph{not} change the storage 282 layout of the array \code{A}. 283 284 \kitem{sym\_grad} If set to \code{true} then it is assumed that the 285 symmetric gradient has to be used for the computation of the jump- 286 and Neumann-residuals. The demo-program \code{quasi-stokes.c} uses 287 this feature to implement an error estimator for the Stokes equation 288 with stress boundary conditions. 289 290 \kitem{f} If the first argument of the function pointer 291 \code{f(result,\dots)} is not \nil then the result \emph{must} be 292 stored in the argument \code{result} and \code{f()} must return the 293 base address of the array \code{result}. If \code{result} is \nil, 294 then \code{f()} must store the result in a non-volatile storage area 295 and return the address of that area. 296 297 \kitem{dirichlet\_bndry} A bit-mask marking those parts of the 298 boundary which are subject to Dirichlet boundary conditions, see 299 \secref{S:boundary}. 300 301 \kitem{f} is a pointer to a function for the evaluation of the 302 lower order terms at all quadrature nodes, i.e. 303 $f(x(\lambda), u(\lambda), \nabla u(\lambda))$ ; if \code{f} 304 is a \nil pointer, $f\equiv0$ is assumed; 305 306 \code{f(el\_info, quad, qp, uh\_qp, grd\_uh\_qp)} returns the 307 value of the lower oder terms on element \code{el\_info->el} 308 at the quadrature node \code{quad->lambda[qp]}, where 309 \code{uh\_qp} is the value and \code{grd\_uh\_qp} the 310 gradient (with respect to the Cartesian coordinates) of the 311 discrete solution at that quadrature point. See also 312 \code{f\_flag} below: 313 314 \kitem{f\_flag} specifies whether the function \code{f()} 315 actually needs values of \code{uh\_qp} or \code{grd\_uh\_qp}, 316 \code{f\_flag} may be $0$ or 317 \code{INIT\_UH}\cdx{INIT_UH@{\code{INIT\_UH}}} or 318 \code{INIT\_GRD\_UH}\cdx{INIT_GRD_UH@{\code{INIT\_GRD\_UH}}} 319 or their bitwise composition (\code{|}). The arguments 320 \code{uh\_qp} and \code{grd\_uh\_qp} of \code{f()} only hold 321 valid information if the flags \code{INIT\_UH} respectively 322 \code{INIT\_GRD\_UH} are set. 323 324 \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to a 325 function for the evaluation of non-homogeneous Neumann boundary 326 data. \code{gn} may be \nil, in which case zero Neumann boundary 327 conditions are assumed. The argument \code{normal} always 328 contains the normal of the Neumann boundary facet. In the case of 329 non-vanishing co-dimension \code{normal} lies in the 330 lower-dimensional space which is spanned by the mesh simplex 331 defined by \code{el\_info}. \code{gn()} is evaluated on those 332 parts of the boundary which are \emph{not} flagged as 333 Dirichlet-boundaries by the argument \code{dirichlet\_bndry}. 334 335 \kitem{gn\_flag} controls whether the argument \code{uh\_qp} 336 of the function \code{gn()} actually contains the value of 337 \code{uh} at the quadrature point \code{qp}. Note that the 338 argument \code{normal} always contains valid data. 339\end{descr} 340\end{descr} 341 342\begin{example}[Linear problem]\label{E:est-impl} 343Consider the scalar linear model problem \mathref{book:E:strong} with constant 344coefficients $A$, $b$, and $c$: 345\begin{alignat*}{2} 346-\nabla \cdot A \nabla u + b \cdot \nabla u + c\, u &= r \qquad & 347&\mbox{in } \Omega,\\ 348u &= 0 & &\mbox{on } \partial\Omega. 349\end{alignat*} 350Let \code{A} be a \code{REAL\_DD} matrix storing $A$, which is then the 351eighth argument of \code{ellipt\_est()}. Assume that 352\code{const REAL *b(const REAL\_D)} is a function returning a pointer 353to a vector storing $b$, \code{REAL c(REAL\_D)} returns the 354value of $c$ and \code{REAL r(const REAL\_D)} returns the value of the 355right hand side $r$ of \mathref{book:E:strong} at some point in world 356coordinates. The implementation of the function \code{f} is: 357\bv\begin{verbatim} 358static REAL f(const EL_INFO *el_info, const QUAD *quad, int iq, REAL uh_iq, 359 const REAL_D grd_uh_iq) 360{ 361 FUNCNAME("f"); 362 const REAL *bx, *x; 363 extern const REAL b(const REAL_D); 364 extern REAL c(const REAL_D), r(const REAL_D); 365 366 x = coord_to_world(el_info, quad->lambda[iq], nil); 367 bx = b(x); 368 369 return(SCP_DOW(bx, grd_uh_iq) + c(x)*uh_iq - r(x)); 370} 371\end{verbatim}\ev 372As both \code{uh\_iq} and \code{grd\_uh\_iq} are used, the estimator parameter 373\code{f\_flag} must be given as \code{INIT\_UH|INIT\_GRD\_UH}. 374\end{example} 375 376%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 377 378\subsection{Error estimator for parabolic problems}\label{S:para_est} 379 380Similar to the stationary case, the \ALBERTA library provides 381an error estimator for the non--linear parabolic problem 382\begin{align*} 383 \partial_t u -\nabla \cdot A \nabla u(x) 384 + f\Big(x, t, u(x),\nabla u(x)\Big) &= 0 &&x \in \Omega, t>0,\\ 385 u(x,t) &= g_d && x \in \Gamma_D, t>0,\\ 386 \nu\cdot A\nabla u(x,t) &= g_n && x \in \Gamma_N, t>0,\\ 387 u(x,0) &= u_0 && x \in \Omega, 388\end{align*} 389where $A \in \R^{d\times d}$ is a positive definite matrix and 390$\partial\Omega=\Gamma_D\cup\Gamma_N$. The estimator is split in 391several parts, where the initial error 392\[ 393\eta_0 = \|u_0 - U_0\|_{L^2(\Omega)} 394\] 395can be approximated by the function \code{L2\_err()}, e.g. (compare 396\secref{S:error_cal}). 397 398For the estimation of the spatial discretization error, the coarsening 399error, and the time discretization error, the \ALBERTA estimator 400implements the following (local) indicators 401\begin{align*} 402 \eta_S^2 &= 403 C_0^2\, h_S^4\, \left\| \frac{U_{n+1} - I_{n+1}U_{n}}{\tau_{n+1}} 404 -\nabla \cdot A \nabla U_{n+1} + f(.,t_{n+1},U_{n+1},\nabla U_{n+1}) 405 \right\|_{L^2(S)}^2\\ 406 &\qquad + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Omega} 407 h_S^3\, \|\jump{A\nabla U_{n+1}}\|_{L^2(\Gamma)}^2 408 + C_1^2\, \sum_{\Gamma\subset\partial S \cap \Gamma_N} 409 h_S^3\, \|\nu\cdot A\nabla U_{n+1}-g_n\|_{L^2(\Gamma)}^2, 410 \\[2mm] 411 \eta_{S,c}^2 &= C_2^2\, 412 h_S^3\, \|\jump{\nabla U_{n}}\|_{L^2(\Gamma_c)}^2 413% \|U_{n} - I_{n+1}U_{n}\|_{L^2(S)}^2 414 \\[2mm] 415 \eta_\tau &= C_3 \|U_{n+1} - I_{n+1}U_{n}\|_{L^2(\Omega)}. 416\end{align*} 417The coarsening indicator is motivated by the fact that for piecewise 418linear Lagrange finite element functions it holds $\|U_{n} - 419I_{n+1}U_{n}\|_{L^2(S)}^2 = \eta_{S,c}^2$ with $C_2=C_2(d)$ and 420$\Gamma_c$ the face that would be removed during a coarsening operation. 421%% 422The implementation is done by the functions 423\fdx{heat_est()@{\code{heat\_est()}}}% 424\fdx{heat_est_dow()@{\code{heat\_est\_dow()}}}% 425\fdx{heat_est_d()@{\code{heat\_est\_d()}}}% 426\idx{error estimators!heat_est()@{\code{heat\_est()}}}% 427\idx{error estimators!heat_est_dow()@{\code{heat\_est\_dow()}}}% 428\idx{error estimators!heat_est_d()@{\code{heat\_est\_d()}}}% 429\bv\begin{verbatim} 430REAL heat_est(const DOF_REAL_VEC *uh, ADAPT_INSTAT *adapt, 431 REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), 432 int quad_degree, REAL C[4], const DOF_REAL_VEC *uh_old, 433 const REAL_DD A, const BNDRY_FLAGS dirichlet_bndry, 434 REAL (*f)(const EL_INFO *el_info, const QUAD *quad, int qp, 435 REAL uh_qp, const REAL_D grd_uh_gp, REAL time), 436 FLAGS f_flags, 437 REAL (*gn)(const EL_INFO *el_info, const QUAD *quad, int qp, 438 REAL uh_qp, const REAL_D normal, REAL time), 439 FLAGS gn_flags); 440REAL heat_est_dow(const DOF_REAL_D_VEC *uh, ADAPT_INSTAT *adapt, 441 REAL *(*rw_est)(EL *), REAL *(*rw_estc)(EL *), 442 int quad_degree, REAL C[4], const DOF_REAL_D_VEC *uh_old, 443 const void *A, MATENT_TYPE A_type, MATENT_TYPE A_blocktype, 444 bool sym_grad, 445 BNDRY_FLAGS dirichlet_bndry, 446 const REAL *(*f)(REAL_D result, 447 const EL_INFO *el_info, 448 const QUAD *quad, int qp, 449 const REAL_D uh_qp, 450 const REAL_DD grd_uh_gp, 451 REAL time), 452 FLAGS f_flags, 453 const REAL *(*gn)(REAL_D result, 454 const EL_INFO *el_info, 455 const QUAD *quad, int qp, 456 const REAL_D uh_qp, 457 const REAL_D normal, 458 REAL time), 459 FLAGS gn_flags); 460REAL heat_est_d(const DOF_REAL_D_VEC *uh, 461 const DOF_REAL_D_VEC *uh_old, 462 ADAPT_INSTAT *adapt, 463 REAL *(*rw_est)(EL *), 464 REAL *(*rw_estc)(EL *), 465 int quad_degree, 466 REAL C[4], 467 const void *A, 468 MATENT_TYPE A_type, 469 MATENT_TYPE A_blocktype, 470 bool sym_grad, 471 const BNDRY_FLAGS dirichlet_bndry, 472 const REAL *(*f)(REAL_D result, 473 const EL_INFO *el_info, 474 const QUAD *quad, 475 int qp, 476 const REAL_D uh_qp, 477 const REAL_DD grd_uh_gp, 478 REAL time), 479 FLAGS f_flags, 480 const REAL *(*gn)(REAL_D result, 481 const EL_INFO *el_info, 482 const QUAD *quad, 483 int qp, 484 const REAL_D uh_qp, 485 const REAL_D normal, 486 REAL time), 487 FLAGS gn_flags); 488\end{verbatim}\ev 489Description: 490\begin{descr} 491\kitem{heat\_est(uh, adapt, rw\_el\_est, rw\_el\_estc, degree, C, uh\_old,} 492\kitem{~~~~~~~~~A, dirichlet\_bndry, f, f\_flag, gn, gn\_flag)}\hfill 493 494 computes an error estimate of the above type, the local and 495 global space discretization estimators are stored in 496 \code{adapt->adapt\_space} and via the \code{rw\_...} pointers; 497 the return value is the time discretization indicator $\eta_\tau$. 498 499 \begin{descr} 500 \kitem{uh} is a vector storing the coefficients of the 501 discrete solution $U_{n+1}$; if \code{uh} is a \nil pointer, 502 nothing is done, the return value is \code{0.0}. 503 504 \kitem{adapt} is a pointer to an \code{ADAPT\_INSTAT} 505 structure; if it is not \nil, then the entries 506 \code{adapt\_space->p=2}, \code{adapt\_space->err\_sum} and 507 \code{adapt\_space->err\_max} of \code{adapt} are set by 508 \code{heat\_est()} (compare 509 \secref{S:adapt_stat_in_ALBERTA}). 510 511 \kitem{rw\_el\_est} is a function for writing the local error 512 indicator $\eta_S^2$ for a single element (usually to some 513 location inside \code{leaf\_data}, compare 514 \secref{S:leaf_data_info}); if this function is \nil, only 515 the global estimate is computed, no local indicators are 516 stored. \code{rw\_el\_est(el)} returns for each leaf element 517 \code{el} a pointer to a \code{REAL} for storing the square 518 of the element indicator, which can directly be used in the 519 adaptive method, compare the \code{get\_el\_est()} function 520 pointer in the \code{ADAPT\_STAT} structure (compare 521 \secref{S:adapt_stat_in_ALBERTA}). 522 523 \kitem{rw\_el\_estc} is a function for writing the local 524 coarsening error indicator $\eta_{S,c}^2$ for a single 525 element (usually to some location inside \code{leaf\_data}, 526 compare \secref{S:leaf_data_info}); if this function is \nil, 527 no coarsening error indicators are computed and stored; 528 \code{rw\_el\_estc(el)} returns for each leaf element 529 \code{el} a pointer to a \code{REAL} for storing the square 530 of the element coarsening error indicator. The coarsening 531 indicator is not used at the moment. 532 533 \kitem{degree} is the degree of the quadrature that should be 534 used for the approximation of the norms on the elements and 535 edges/faces; if \code{degree} is less than zero a quadrature 536 which is exact of degree 537 \code{2*uh->fe\_space->bas\_fcts->degree} is used. 538 539 \kitem{C[0]}, \code{C[1]}, \code{C[2]}, \code{C[3]} are the 540 constants in front of the element residual, wall residual, 541 coarsening term, and time residual, respectively. If \code{C} 542 is \nil, then all constants are set to $1.0$. 543 544 \kitem{uh\_old} is a vector storing the coefficients of the 545 discrete solution $U_{n}$ from previous time step; if 546 \code{uh\_old} is a \nil pointer, nothing is done, the return 547 value is \code{0.0}. 548 549 \kitem{A} is the constant matrix of the second order term. 550 551 \kitem{dirichlet\_bndry} A bit mask marking those parts of the 552 boundary which are subject to Dirichlet boundary conditions. See 553 \secref{S:boundary}. 554 555 \kitem{f} is a pointer to a function for the evaluation of the 556 lower order terms at all quadrature nodes, i.e. 557 $f(x(\lambda), t, u(\lambda), \nabla u(\lambda))$ ; if 558 \code{f} is a \nil pointer, $f\equiv0$ is assumed; 559 560 \code{f(el\_info, quad, iq, t, uh\_iq, grd\_uh\_iq)} returns 561 the value of the lower oder terms on element 562 \code{el\_info->el} at the quadrature node 563 \code{quad->lambda[iq]}, where \code{uh\_iq} is the value and 564 \code{grd\_uh\_iq} the gradient (with respect to the world 565 coordinates) of the discrete solution at that quadrature 566 node. 567 568 \kitem{f\_flag} specifies whether the function \code{f()} 569 actually needs values of \code{uh\_iq} or \code{grd\_uh\_iq}. 570 This flag may hold zero, the predefined values 571 \code{INIT\_UH} or \code{INIT\_GRD\_UH}, or their composition 572 \code{INIT\_UH|INIT\_GRD\_UH}; the arguments \code{uh\_iq} 573 and \code{grd\_uh\_iq} of \code{f()} only hold valid 574 information, if the flags \code{INIT\_UH} respectively 575 \code{INIT\_GRD\_UH} are set. 576 577 \kitem{gn(el\_info, quad, qp, uh\_qp, normal)} is a pointer to 578 a function for the evaluation of non-homogeneous Neumann 579 boundary data. \code{gn} may be \nil, in which case zero 580 Neumann boundary conditions are assumed. The argument 581 \code{normal} always contains the normal of the Neumann 582 boundary facet. In the case of non-vanishing co-dimension 583 \code{normal} lies in the lower-dimensional space which is 584 spanned by the mesh simplex defined by \code{el\_info}. 585 586 \kitem{gn\_flag} controls whether the argument \code{uh\_qp} 587 of the function \code{gn()} actually contains the value of 588 \code{uh} at the quadrature point \code{qp}. Note that the 589 argument \code{normal} always contains valid data. 590 \end{descr} 591 \smallskip 592 593 The estimate is computed by traversing all leaf elements of 594 \code{uh->fe\_space->mesh}, using the quadrature for the 595 approximation of the residuals and storing the square of the 596 element indicators on the elements (if \code{rw\_el\_est} and 597 \code{rw\_el\_estc }are not \nil). 598 599\kitem{heat\_est\_d(uh, adapt, rw, rwc, deg, C, uh\_old,} 600\kitem{~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} 601\kitem{~~~~~~~~~~dirichlet\_bndry, f, f\_flag)}\hfill 602\kitem{heat\_est\_dow(uh, adapt, rw, rwc, deg, C, uh\_old,} 603\kitem{~~~~~~~~~~~~A, A\_type, A\_blocktype, sym\_grad,} 604\kitem{~~~~~~~~~~~~dirichlet\_bndry, f, f\_flag)}\hfill 605 606 Coupled vector valued version. See \code{ellipt\_est\_dow()} above. 607\end{descr} 608 609There are also some less high-level support functions which allow for 610custom contributions to the per-element error estimates. We will not 611document this in detail, but rather refer the reader to the 612\code{stokes.c} and \code{quasi-stokes.c} demo-programs. 613%% 614\fdx{ellipt_est_init()@{\code{ellipt\_est\_init()}}} 615\fdx{heat_est_init()@{\code{heat\_est\_init()}}} 616\fdx{element_est()@{\code{element\_est()}}} 617\fdx{element_est_finish()@{\code{element\_est\_finish()}}} 618\fdx{element_est_uh()@{\code{element\_est\_uh()}}} 619\fdx{element_est_grd_uh()@{\code{element\_est\_grd\_uh()}}} 620\fdx{ellipt_est_finish()@{\code{ellipt\_est\_finish()}}} 621\fdx{heat_est_finish()@{\code{heat\_est\_finish()}}} 622%% 623\fdx{ellipt_est_dow_init()@{\code{ellipt\_est\_dow\_init()}}} 624\fdx{heat_est_dow_init()@{\code{heat\_est\_dow\_init()}}} 625\fdx{element_est_dow()@{\code{element\_est\_dow()}}} 626\fdx{element_est_dow_finish()@{\code{element\_est\_dow\_finish()}}} 627\fdx{element_est_uh_dow()@{\code{element\_est\_uh\_dow()}}} 628\fdx{element_est_grd_uh_dow()@{\code{element\_est\_grd\_uh\_dow()}}} 629\fdx{ellipt_est_dow_finish()@{\code{ellipt\_est\_dow\_finish()}}} 630\fdx{heat_est_dow_finish()@{\code{heat\_est\_dow\_finish()}}} 631%% 632\bv\begin{lstlisting} 633const void *ellipt_est_init(const DOF_REAL_VEC *uh, 634 ADAPT_STAT *adapt, 635 REAL *(*rw_est)(EL *), 636 REAL *(*rw_estc)(EL *), 637 const QUAD *quad, 638 const WALL_QUAD *wall_quad, 639 NORM norm, 640 REAL C[3], 641 const REAL_DD A, 642 const BNDRY_FLAGS dirichlet_bndry, 643 REAL (*f)(const EL_INFO *el_info, 644 const QUAD *quad, 645 int qp, 646 REAL uh_qp, 647 const REAL_D grd_uh_gp), 648 FLAGS f_flags, 649 REAL (*gn)(const EL_INFO *el_info, 650 const QUAD *quad, 651 int qp, 652 REAL uh_qp, 653 const REAL_D normal), 654 FLAGS gn_flags); 655const void *heat_est_init(const DOF_REAL_VEC *uh, 656 const DOF_REAL_VEC *uh_old, 657 ADAPT_INSTAT *adapt, 658 REAL *(*rw_est)(EL *), 659 REAL *(*rw_estc)(EL *), 660 const QUAD *quad, 661 const WALL_QUAD *wall_quad, 662 REAL C[4], 663 const REAL_DD A, 664 const BNDRY_FLAGS dirichlet_bndry, 665 REAL (*f)(const EL_INFO *el_info, 666 const QUAD *quad, 667 int qp, 668 REAL uh_qp, 669 const REAL_D grd_uh_gp, 670 REAL time), 671 FLAGS f_flags, 672 REAL (*gn)(const EL_INFO *el_info, 673 const QUAD *quad, 674 int qp, 675 REAL uh_qp, 676 const REAL_D normal, 677 REAL time), 678 FLAGS gn_flags); 679 680REAL element_est(const EL_INFO *el_info, const void *est_handle); 681void element_est_finish(const EL_INFO *el_info, 682 REAL est_el, const void *est_handle); 683const REAL *element_est_uh(const void *est_handle); 684const REAL_D *element_est_grd_uh(const void *est_handle); 685REAL ellipt_est_finish(ADAPT_STAT *adapt, const void *est_handle); 686REAL heat_est_finish(ADAPT_INSTAT *adapt, const void *est_handle); 687\end{lstlisting}\ev 688%% 689There are similar proto-types for the vector-valued case. Now, what 690are these functions good for? The \code{stokes.c} program makes use of 691this framework to add a contribution concerning the divergence 692constraint. Of course, this is an ad-hoc error indicator, and only 693meant to demonstrate the programming frame-work. The functions 694\code{element\_est\_uh[\_dow]()} and 695\code{element\_est\_grd\_uh[\_dow]()} give the application access to 696the values of the discrete solution at the quadrature points 697(respectively to its Jaocbians). Otherwise, the general layout is like 698follows: 699%% 700\bv\begin{lstlisting} 701void *est_handle = ellipt_est_init(...); 702TRAVERSE_FIRST(mesh, -1, <suitable fill-flags>) { 703 REAL est_el = element_est(el_info, est_handle); 704 705 ... /* add whatever you like to est_el */ 706 707 element_est_finish(el_info, est_el, est_handle); 708} TRAVERSE_NEXT(); 709REAL est = ellipt_est_finish(adapt, est_handle); 710\end{lstlisting}\ev 711%% 712The relevant excerpt from \code{stokes.c} reads as follows: 713%% 714\bv\begin{lstlisting} 715 est_handle = ellipt_est_dow_init(u_h, adapt, rw_el_est, NULL /* rw_estc */, 716 quad, NULL /* wall_quad */, 717 H1_NORM, C, 718 A, MATENT_REAL, MATENT_REAL, 719 false /* !sym_grad */, 720 dirichlet_mask, 721 r, INIT_GRD_UH, 722 NULL /* inhomog. Neumann res. */, 0); 723 724 fill_flags = FILL_NEIGH|FILL_COORDS|FILL_OPP_COORDS|FILL_BOUND|CALL_LEAF_EL; 725 fill_flags |= u_fe_space->bas_fcts->fill_flags; 726 fill_flags |= p_fe_space->bas_fcts->fill_flags; 727 TRAVERSE_FIRST(mesh, -1, fill_flags) { 728 const EL_GEOM_CACHE *elgc; 729 const QUAD_EL_CACHE *qelc; 730 REAL est_el; 731 732 est_el = element_est_dow(el_info, est_handle); 733 734 if (C[3]) { 735 REAL div_uh_el, div_uh_qp; 736 const REAL_DD *grd_uh_qp; 737 int qp, i; 738 739 grd_uh_qp = element_est_grd_uh_d(est_handle); 740 div_uh_el = 0.0; 741 if (!(el_info->fill_flag & FILL_COORDS)) { 742 qelc = fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_DET); 743 744 for (qp = 0; qp < quad->n_points; qp++) { 745 div_uh_qp = 0; 746 for (i = 0; i < DIM_OF_WORLD; i++) { 747 div_uh_qp += grd_uh_qp[qp][i][i]; 748 } 749 div_uh_el += qelc->param.det[qp]*quad->w[qp]*SQR(div_uh_qp); 750 } 751 } else { 752 elgc = fill_el_geom_cache(el_info, FILL_EL_DET); 753 754 for (qp = 0; qp < quad->n_points; qp++) { 755 div_uh_qp = 0; 756 for (i = 0; i < DIM_OF_WORLD; i++) { 757 div_uh_qp += grd_uh_qp[qp][i][i]; 758 } 759 div_uh_el += quad->w[qp]*SQR(div_uh_qp); 760 } 761 div_uh_el *= elgc->det; 762 } 763 764 est_el += C[3] * div_uh_el; 765 } 766 767 element_est_dow_finish(el_info, est_el, est_handle); 768 } TRAVERSE_NEXT(); 769 est = ellipt_est_dow_finish(adapt, est_handle); 770\end{lstlisting}\ev 771 772\idx{error estimators|)} 773 774 775%%% Local Variables: 776%%% mode: latex 777%%% TeX-master: "alberta-man" 778%%% End: 779