1\section{Finite element spaces and finite element discretization}% 2\label{S:FEM} 3 4\def\vI{\vec{I}} 5 6In the sequel we assume that $\Omega \subset \R^d$ is a bounded 7domain triangulated by $\tria$, i.e. 8\[ 9\bar \Omega \; = \; \bigcup_{S \in \tria} S. 10\] 11 12The following considerations are also valid for a triangulation of an 13immersed surface (with $n > d$). In this situation one has to exchange 14derivatives (those with respect to $x$) by \emph{tangential} 15derivatives (tangential to the actual element, derivatives are always 16taken element--wise) and the determinant of the parameterization's 17Jacobian has to be replaced by Gram's determinant of the 18parameterization. But for the sake of clearness and simplicity we 19restrict our considerations to the case $n = d$. 20 21The values of a finite element function or the values of its 22derivatives are uniquely defined by the values of its DOFs and the 23values of the basis functions or the derivatives of the basis 24functions connected with these DOFs. Usually, evaluation of those 25values is performed element--wise. On a single element the value of a 26finite element function at a point $x$ in this element is 27determined by the DOFs associated with this specific element and the 28values of the non vanishing basis functions at this point. 29 30We follow the concept of finite elements which are given on a 31single element $S$ in local coordinates. We distinguish two 32classes of finite elements: 33 34Finite element functions on an element $S$ defined by a finite 35dimensional function space $\Pbar$ on a \emph{reference element} 36$\Sbar$ and the (one to one) mapping $\lS$ from the reference element 37$\Sbar$ to $S$. For this class the dependency on the actual element 38$S$ is fully described by the mapping $\lS$. For example, all Lagrange 39finite elements belong to this class. 40 41Secondly, finite element functions depending on the actual element 42$S$. Hence, the basis functions are not fully described by $\Pbar$ and 43the one to one mapping $\lS$. But using an initialization of the 44actual element (which defines a finite dimensional function space 45$\Pbar$ with information about the actual element), we can implement 46this class in the same way as the first one. This class is needed for 47Hermite finite elements which are not affine equivalent to the 48reference element. Examples in 2d are the \emph{Hsieh--Clough--Tocher} 49or \emph{HCT element} or the \emph{Argyris element} where only 50the normal derivative at the midpoint of edges are used 51in the definition of finite element functions; both 52elements lead to functions which are globally $C^1(\bar\Omega)$. The 53concrete implementation for this class in \ALBERTA is future work. 54 55All in all, for a very general situation, we only need a vector of 56basis functions (and their derivatives) on $\Sbar$ and a function which 57connects each of these basis functions with its degree of freedom on 58the element. For the second class, we additionally need an 59initialization routine for the actual element. By such information, 60every finite element function is uniquely described on every element 61of the grid. 62 63\subsection{Barycentric coordinates}\label{S:bary_coord}% 64\idx{barycentric coordinates|(} 65 66For describing finite elements on simplicial grids, it is very 67convenient to use $d+1$ barycentric coordinates as a local coordinate 68system on an element of the triangulation. Using $d+1$ 69local coordinates, the \emph{reference simplex} $\Sbar$ is a subset of 70a hyper surface in $\R^{d+1}$: 71\begin{equation*}\label{S:Sbar} 72\Sbar := \big\{ (\lambda_0,\dots,\lambda_d) 73\in \R^{d+1}; \lambda_k \ge 0,\; \sum\limits_{k = 0}^d \lambda_k = 1 \big\} 74\end{equation*} 75On the other hand, for numerical integration on an 76element it is much more convenient to use the standard element 77$\Shat \in \R^d$ defined in \secref{S:ref.coarse} as 78\[ 79\Shat = \mbox{conv hull} 80\left\{\hat a_0 = 0, \hat a_1 = e_1, \dots, \hat a_d = e_d\right\} 81\] 82where $e_i$ are the unit vectors in $\R^d$; using $\Shat$ for 83the numerical integration, we only have to compute the determinant of 84the parameterization's Jacobian and not Gram's determinant. 85 86The relation between a given simplex $S$, the reference simplex $\Sbar$, 87and the standard simplex $\Shat$ is now described in detail. 88 89Let $S$ be an element of the triangulation with vertices 90$\{a_0,\dots,a_d\}$; let $F_S\colon \Shat \to S$ be the diffeomorphic 91parameterization of $S$ over $\Shat$ with regular Jacobian $DF_S$, such 92that 93\[ 94 F_S(\hat a_k) = a_k, \qquad k = 0,\dots,d 95\] 96holds. For a point $x \in S$ we set 97\begin{equation*}\label{E:xhat} 98\hat x = F_S^{-1}(x) \in \Shat. 99\end{equation*} 100For a simplex $S$ the easiest choice for $F_S$ is the unique affine 101mapping \mathref{E:FS} defined on page \pageref{E:FS}. For an affine mapping, 102$DF_S$ is constant. In the following, we assume that the 103parameterization $F_S$ of a simplex $S$ is affine. 104 105For a simplex $S$ the barycentric coordinates 106\[ 107\lS(x) \; = \; (\lS_0,\dots,\lS_d)(x) \in \R^{d+1} 108\] 109of some point $x \in \R^d$ are (uniquely) determined by the 110$(d+1)$ equations 111\[ 112\sum\limits_{k = 0}^d \lS_k(x)\; a_k = x 113\qquad\mbox{and}\qquad 114\sum\limits_{k = 0}^d \lS_k(x) = 1. 115\] 116The following relation holds: 117\[ 118x \in S \qquad \mbox{\iff} \qquad \lS_k(x) \in [0,1]\ \ \mbox{for all } 119k = 0,\dots,d \qquad \mbox{iff} \qquad \lS \in \Sbar. 120\] 121On the other hand, each $\lambda \in \Sbar$ defines a unique point 122$\xS \in S$ by 123\[ 124\xS(\lambda) = \sum\limits_{k = 0}^d \lambda_k\, a_k. 125\] 126Thus, $\xS\colon \Sbar \to S$ is an invertible mapping with inverse 127$\lS\colon S \to \Sbar$. The barycentric coordinates of $x$ on $S$ 128are the same as those of $\hat x$ on $\Shat$, i.e. $\lS(x) = 129\lShat(\hat x)$. 130 131In the general situation, when $F_S$ may not be affine, i.e.\ 132we have a parametric element, the barycentric coordinates $\lS$ are 133given by the inverse of the parameterization $F_S$ and 134the barycentric coordinates on $\Shat$: 135\begin{equation*}\label{E:lS} 136\lS(x) = \lShat(\hat x) = \lShat\left(F_S^{-1}(x)\right) 137\end{equation*} 138and the \emph{world coordinates} of a point $\xS \in S$ with barycentric 139coordinates $\lambda$ are given by 140\begin{equation*}\label{E:xS} 141\xS(\lambda) \; = \; F_S\left(\sum_{k = 0}^d \lambda_k \hat a_k\right) 142 \; = \; F_S\left(\xShat(\lambda)\right) 143\end{equation*} 144(see also Figure~\ref{F:S_Shat_Sb}). 145 146\begin{figure}[htbp] 147\centerline{\includegraphics[scale=1.0]{EPS/bary}} 148\caption{Definition of $\lS \colon S \to \Sbar$ via $F_S^{-1}$ and $\lShat$, 149and $\xS\colon \Sbar \to S$ via $\xShat$ and $F_S$}\label{F:S_Shat_Sb} 150\end{figure} 151 152Every function $f \colon S \to V$ defines (uniquely) two functions 153\[ 154\begin{array}{rcl} 155\bar f \colon\, \Sbar &\to & V\\ 156 \lambda & \mapsto & f(\xS(\lambda)) 157\end{array} 158\qquad\mbox{and}\qquad 159\begin{array}{rcl} 160\hat f \colon\, \Shat & \to &V\\ 161 \xhat & \mapsto & f(F_S(\xhat)). 162\end{array} 163\] 164Accordingly, $\hat f \colon\, \Shat \to V$ defines two functions 165$f \colon\, S \to V$ and $\bar f \colon\, \Sbar \to V$, and 166$\bar f \colon\, \Sbar \to V$ defines $f \colon\, S \to V$ and 167$\hat f \colon\, \Shat \to V$. 168 169Assuming that a function space $\Pbar \subset C^0(\Sbar)$ is given, 170it uniquely defines function spaces $\Phat$ and $\PS$ by 171\begin{equation}\label{Pbar_PS} 172\Phat = \left\{\phat \in C^0(\Shat);\, \pbar \in \Pbar\right\} 173\qquad \mbox{and} \qquad 174\PS = \left\{\vphi \in C^0(S);\, \pbar \in \Pbar\right\}. 175\end{equation} 176We can also assume that the function space $\Phat$ is given and this 177space uniquely defines $\Pbar$ and $\P_S$ in the same manner. In 178\ALBERTA, we use the function space $\Pbar$ on $\Sbar$; the 179implementation of a basis $\left\{\pbar^1,\dots,\pbar^m\right\}$ of 180$\Pbar$ is much simpler than the implementation of a basis 181$\left\{\phat^1,\dots,\phat^m\right\}$ of $\Phat$ as we are able to 182use symmetry properties of the barycentric coordinates $\lambda$. 183 184In the following we shall often drop the superscript $S$ of $\lS$ 185and $\xS$. The mappings $\lambda(x) = \lS(x)$ and $x(\lambda) = 186\xS(\lambda)$ are always defined with respect to the actual element 187$S \in \tria$. 188\idx{barycentric coordinates|)} 189 190\subsection{Finite element spaces}% 191\label{S:FES}% 192\idx{finite element spaces|(} 193 194\ALBERTA supports scalar and vector valued finite element functions. 195The basis functions are always real valued; thus, the coefficients 196of a finite element function in a representation by the basis functions 197are either scalars or vectors of length $n$. 198 199For a given function space $\Pbar$ and some given real valued function space 200$C$ on $\Omega$, a finite element space $X_h$ is defined by 201\begin{equation*}\label{E:Xh} 202X_h = X_h(\tria, \Pbar, C) = 203\left\{\vphi \in C; \; \vphi_{|S} \in \P_S \mbox{ for all } S \in \tria\right\} 204\end{equation*} 205for scalar finite elements. For vector valued finite elements, $X_h$ is 206given by 207\begin{align*}\label{E:Xhv} 208X_h &= X_h(\tria, \Pbar, C) 209%\notag\\ & 210= \left\{ 211\vphi = (\vphi_1,\dots,\vphi_n) \in C^n; 212\, \vphi_{i\,|S} \in \P_S \mbox{ for all } i=1,\dots,n,\, S \in \tria\right\}. 213\end{align*} 214The spaces $\P_S$ are defined by $\Pbar$ via \mathref{Pbar_PS}. 215 216For conforming finite element discretizations, $C$ is the continuous 217space $X$ (for 2nd order problems, $C = X = H^1(\Omega)$). For 218non--conforming finite element discretizations, $C$ may control the 219\emph{non conformity} of the ansatz space which has to be controlled in 220order to obtain optimal error estimates (for example, in the 221discretization of the incompressible Navier--Stokes equation by the 222non--conforming Crouzeix--Raviart element, the finite element 223functions are continuous only in the midpoints of the edges). 224 225The abstract definition of finite element spaces as a triple (mesh, 226local basis functions, and DOFs) now matches the mathematical 227definition as $X_h = X_h(\tria,\Pbar,C)$ in the following way: The 228mesh is the underlying triangulation $\tria$, the local basis 229functions are given by the function space $\Pbar$, and together with 230the relation between global and local degrees of freedom every finite 231element function satisfies the global continuity constraint 232$C$. This relation between global and local DOFs is now described in 233detail.\idx{finite element spaces|)} 234 235\subsection{Evaluation of finite element functions}% 236\label{S:eval_fe}% 237\idx{DOFs!relation global and local DOFs}% 238\idx{evaluation of finite element functions} 239 240Let $\left\{\pbar^1,\dots,\pbar^m\right\}$ be a basis of $\Pbar$ and 241let $\left\{\vphi_1,\dots,\vphi_N\right\}$ be a basis of $X_h$, $N = 242\dim X_h$, such that for every $S \in \tria$ and for all 243basis functions $\vphi_j$ which do not vanish on $S$ 244\[ 245\vphi_j|_S(x(\lambda)) = \pbar^i(\lambda) \qquad \mbox{for all } \lambda 246\in \Sbar 247\] 248holds with some $i \in \{1,\dots,m\}$ depending on $j$ and $S$. Thus, the 249numbering of the basis functions in $\Pbar$ and the mapping $\xS$ 250induces a local numbering of the non vanishing basis functions on 251$S$. Denoting by $J_S$ the index set of all non vanishing basis 252functions on $S$, the mapping $i_S: J_S \to \{1,\dots,m\}$ is one to one 253and uniquely connects the degrees of freedom of a finite element 254function on $S$ with the local numbering of basis functions. 255If $j_S: \{1,\dots,m\} \to J_S$ denotes the inverse mapping of $i_S$, 256the connection between global and local basis functions 257is uniquely determined on each element $S$ by 258\begin{subequations}\label{E:global-lokal} 259\begin{alignat}{2} 260\vphi_j(x(\lambda)) &= \pbar^{i_S(j)}(\lambda), 261&\qquad & \mbox{for all } \lambda \in \Sbar,\,j \in J_S, 262\label{E:global-lokal.a}\\ 263\vphi_{j_S(i)}(x(\lambda)) &= \pbar^{i} (\lambda), 264&\qquad & \mbox{for all } \lambda \in \Sbar,\,i \in \{1,\dots,m\}. 265\label{E:global-lokal.b} 266\end{alignat} 267\end{subequations} 268 269For $\uh \in X_h$ denote by $\left(u_1,\dots,u_N\right)$ the 270\emph{global coefficient vector} of the basis $\left\{\vphi_j\right\}$ with 271$u_j \in \R$ for scalar finite elements, $u_j \in \R^n$ for 272vector valued finite elements, i.e. 273\[ 274\uh(x) = \sum\limits_{j = 1}^N u_j \vphi_j(x) \qquad 275\mbox{for all } x \in \Omega, 276\] 277and the \emph{local coefficient vector} 278\begin{equation*}\label{E:local_coeff} 279\left(u_S^1,\dots,u_S^m\right) = 280\left(u_{j_S(1)},\dots,u_{j_S(m)}\right) 281\end{equation*} 282of $\uh$ on $S$ with respect to the local numbering of the non 283vanishing basis functions (local numbering is denoted by a superscript 284index and global numbering is denoted by a subscript index). Using the 285local coefficient vector and the local basis functions we obtain the 286\emph{local representation} of $u_h$ on $S$: 287\begin{equation*}\label{E:uh_on_Sx} 288\uh(x) = \sum\limits_{i = 1}^m u^i_S\, \pbar^i(\lambda(x)) 289\qquad \mbox{for all } x \in S. 290\end{equation*} 291In finite element codes, finite element functions are \emph{not} 292evaluated at \emph{world coordinates} $x$ as in \mathref{E:uh_on_Sx}, 293but they are evaluated on a single element $S$ at barycentric coordinates 294$\lambda$ on $S$ giving the value at the world coordinates $x(\lambda)$: 295\begin{equation*}\label{E:uh_on_S} 296\uh(x(\lambda)) = \sum\limits_{i = 1}^m u^i_S\, \pbar^i(\lambda) 297\qquad \mbox{for all } \lambda \in \Sbar. 298\end{equation*} 299The mapping $j_S$, which allows the access to the local 300coefficient vector from the global one, is the relation between 301local DOFs and global DOFs. Global DOFs for a finite element 302space are stored on the mesh elements at positions which are known 303to the DOF administration of the finite element space. Thus, 304the corresponding DOF administration will provide information 305for the implementation of the function $j_S$ and therefore 306information for reading/writing local coefficients from/to global 307coefficient vectors (compare Section~\ref{man:S:basfct_data}). 308 309\subsubsection{Evaluating derivatives of finite element functions}% 310\label{S:eval_Dfe}% 311\idx{evaluation of derivatives} 312 313Derivatives of finite element functions are also evaluated on single 314elements $S$ in barycentric coordinates. In the calculation of 315derivatives in terms of the basis functions $\pbar^i$, the Jacobian 316$\Lambda = \Lambda_S \in \R^{d\times\code{DIM\_OF\_WORLD}}$ 317of the barycentric coordinates on $S$ is involved (we consider here 318only the case $d=\code{DIM\_OF\_WORLD}=n$): 319\begin{equation*}\label{E:LambdaS} 320\Lambda(x) = 321\left( 322\begin{matrix} 323\lambda_{0,x_1}(x) &\lambda_{0,x_2}(x) & \cdots & \lambda_{0,x_n} (x)\\ 324\vdots & \vdots & & \vdots\\ 325\lambda_{d,x_1}(x) &\lambda_{d,x_2}(x) & \cdots & \lambda_{d,x_n} (x)\\ 326\end{matrix} 327\right) 328= 329\left( 330\begin{matrix} 331\mbox{--}&\nabla \lambda_0(x)^t&\mbox{--}\\ 332& \vdots &\\ 333\mbox{--}&\nabla \lambda_d(x)^t&\mbox{--} 334\end{matrix} 335\right). 336\end{equation*} 337Now, using the chain rule we obtain for every function $\vphi \in \PS$ 338\begin{equation*}\label{E:grad_phi} 339\nabla \vphi(x) = \nabla \left(\pbar(\lambda(x))\right) 340= \sum_{k = 0}^d \pbar_{,\lambda_k}(\lambda(x)) \nabla \lambda_k(x) 341= \Lambda^t(x) \nablal \pbar(\lambda(x)), \qquad x \in S 342\end{equation*} 343and 344\begin{equation*}\label{E:D2_phi} 345D^2 \vphi(x) = \Lambda^t(x) D^2_\lambda \pbar(\lambda(x)) \Lambda(x) 346 + \sum_{k = 0}^d D^2 \lambda_k(x) \, \pbar_{,\lambda_k}(\lambda(x)), 347\qquad x \in S. 348\end{equation*} 349For a simplex $S$ with an affine parameterization $F_S$, 350$\nabla \lambda_k$ is constant on $S$ and we get 351\[ 352\nabla \vphi(x) = \Lambda^t \nablal \pbar(\lambda(x)) 353\quad \mbox{and} \quad 354D^2 \vphi(x) = \Lambda^t D^2_\lambda \pbar(\lambda(x)) \Lambda, 355\qquad x \in S. 356\] 357Using these equations, we immediately obtain 358\begin{equation*}\label{E:grd_uh_on_Sx} 359\nabla u_h(x) = \Lambda^t(x) \sum_{i = 1}^m u_S^i \, \nablal 360\pbar^i(\lambda(x)), \qquad x \in S 361\end{equation*} 362and 363\begin{equation*}\label{E:D2_uh_on_Sx} 364D^2 u_h(x) = \Lambda^t(x) \sum_{i = 1}^m u_S^i \, 365D^2_\lambda \pbar^i(\lambda(x))\Lambda(x) 366 + 367\sum_{k = 0}^d D^2 \lambda_k(x) \sum_{i = 1}^m u_S^i \, 368 \pbar^i_{,\lambda_k} (\lambda(x)), \qquad x \in S. 369\end{equation*} 370Since the evaluation is actually done in barycentric coordinates, 371this turns on $S$ into 372\begin{equation*}\label{E:grd_uh_on_S} 373\nabla u_h(x(\lambda)) = \Lambda^t(x(\lambda)) \sum_{i = 1}^m u_S^i \, \nablal 374\pbar^i(\lambda), \qquad \lambda \in \Sbar 375\end{equation*} 376and 377\begin{equation*}\label{E:D2_uh_on_S} 378D^2 u_h(x(\lambda)) = \Lambda^t(x(\lambda)) \sum_{i = 1}^m u_S^i \, 379D^2_\lambda \pbar^i(\lambda)\Lambda(x(\lambda)) 380 + 381\sum_{k = 0}^d D^2 \lambda_k(x(\lambda)) \sum_{i = 1}^m u_S^i \, 382 \pbar^i_{,\lambda_k} (\lambda), \qquad \lambda \in \Sbar. 383\end{equation*} 384 385Once the values of the basis functions, its derivatives, and the local 386coefficient vector $(u_S^1, \dots, u_S^{m})$ are known, the 387evaluation of $u_h$ and its derivatives depends only on $\Lambda$ 388and can be performed by some general evaluation routine (compare 389Section~\ref{man:S:eval}). 390 391\subsection{Interpolation and restriction during refinement and coarsening}% 392\label{S:inter_restrict}% 393\idx{interpolation and restriction of DOF vectors|(} 394 395We assume the following situation: Let $S$ be a (non--parametric) 396simplex with children $S_0$ and $S_1$ generated by the bisection of 397$S$ (compare \algorithmref{A:recursive_refine}). Let $X_S$, 398$X_{S_0,S_1}$ be the finite element spaces restricted to $S$ and $S_0 399\cup S_1$ respectively. 400 401Throughout this section we denote by $\big\{\vphi^i\big\}_{i = 1, \dots, m}$ 402the basis of the coarse grid space $X_S$ and by 403$\big\{\psi^j\}_{j = 1, \dots, k}$ the basis functions of 404$X_{S_0 \cup S_1}$. For a function $\uh \in X_S$ we denote by 405$\vec{u}_\vphi = (u^1_\vphi, \dots, u^m_\vphi)^t$ 406the coefficient vector with respect to the basis $\big\{\vphi^i\big\}$ 407and for a function $\vh \in X_{S_0 \cup S_1}$ by 408$\vec{v}_\psi = (v^1_\psi, \dots, v^k_\psi)^t$ the 409coefficient vector with respect to $\big\{\psi^j\big\}$. 410 411We now derive equations for the transformation of local coefficient 412vectors for finite element functions that are interpolated during 413refinement and coarsening, and for vectors storing values of a 414linear functional applied to the basis functions on the fine grid 415which are restricted to the coarse functions during coarsening. 416 417Let 418\[ 419\I^S_{S_0,S_1} \colon\, X_S \to X_{S_0 \cup S_1} 420\] 421be an interpolation operator. For nested finite element spaces, i.e. 422$X_S \subset X_{S_0 \cup S_1}$, every coarse grid function $\uh \in 423X_S$ belongs also to $X_{S_0 \cup S_1}$, so the natural choice is 424$\I^S_{S_0,S_1} = id$ on $X_S$ (for example, Lagrange finite elements 425are nested). The interpolants $\I^S_{S_0,S_1} \vphi^i$ can be written 426in terms of the fine grid basis functions 427\[ 428\I^S_{S_0,S_1} \vphi^i = \sum_{j = 1}^k a_{ij} \psi^j 429\] 430defining the $(m \times k)$--matrix 431\begin{equation}\label{E:a_matrix} 432A = (a_{ij})_{\atop{i = 1,\dots,m}{j = 1,\dots,k}}. 433\end{equation} 434This matrix $A$ is involved in the interpolation during refinement and the 435transformation of a linear functional during coarsening. 436 437For the interpolation of functions during coarsening, we need an 438interpolation operator $\I^{S_0,S_1}_S \colon\, X_{S_0 \cup S_1} \to X_S$. 439The interpolants $\I^{S_0,S_1}_S \psi^j$ of the fine grid 440functions can now be represented by the coarse grid basis 441\[ 442 \I^{S_0,S_1}_S \psi^j = \sum_{i = 1}^m b_{ij}\, \vphi^i 443\] 444defining the $(m \times k)$--matrix 445\begin{equation}\label{E:b_matrix} 446B = (b_{ij})_{\atop{i = 1,\dots,m}{j = 1,\dots,k}}. 447\end{equation} 448This matrix $B$ is used for the interpolation during coarsening. 449 450Both matrices depend only on the set of local basis functions on 451parent and children. Thus, they depend on the reference element 452$\Sbar$ and one single bisection of the reference element into 453$\Sbar_0$, $\Sbar_1$. The matrices do depend on the local numbering 454of the basis functions on the children with respect to the 455parent. Thus, in 3d the matrices depend on the element type of $S$ 456also (for an element of type \code{0} the numbering of basis functions 457on $\Sbar_1$ differs from the numbering on $\Sbar_1$ for an element of 458type \code{1}, \code{2}). But all matrices can be calculated by the 459local set of basis functions on the reference element. 460 461DOFs can be shared by several elements, compare Section~\ref{S:DOFs1}. 462Every DOF is connected to a basis function which has a support on all 463elements sharing this DOF\@. Each DOF refers to one coefficient of a 464finite element function, and this coefficient has to be calculated 465only once during interpolation. During the restriction of a linear 466functional, contributions from several basis functions are added to 467the coefficient of another basis function. Here we have to control 468that for two DOFs, both shared by several elements, the contribution 469of the basis function at one DOF is only added once to the other DOF 470and vice versa. This can only be done by performing the interpolation 471and restriction on the whole refinement/coarsening patch at the same 472time. 473 474\subsubsection{Interpolation during refinement}% 475\idx{refinement!interpolation of DOF vectors}% 476\idx{interpolation of DOF vectors} 477 478Let $\vec{u}_\vphi = (u^1_\vphi, \dots, u^m_\vphi)^t$ be the coefficient 479vector of a finite element function $\uh \in X_S$ with respect to 480$\big\{\vphi^i\big\}$, 481and let $\vec{u}_\psi = (u^1_\psi, \dots, u^k_\psi)^t$ the 482coefficient vector of $\I^S_{S_0,S_1} \uh$ with respect to 483$\big\{\psi^j\big\}$. Using matrix $A$ defined in \mathref{E:a_matrix} 484we conclude 485\[ 486\I^S_{S_0,S_1} \uh = \sum_{i = 1}^m u^i_\vphi\,\I^S_{S_0,S_1} \vphi^i 487= \sum_{i = 1}^m u^i_\vphi \, \sum_{j = 1}^k a_{ij} \, \psi^j 488 = \sum_{j = 1}^k \left(A^t \vec{u}_\vphi \right)_j \psi^j, 489\] 490or equivalently 491\[ 492\vec{u}_\psi = A^t \vec{u}_\vphi. 493\] 494A subroutine which interpolates a finite element function during 495refinement is an efficient implementation of this matrix--vector 496multiplication. 497 498\subsubsection{Restriction during coarsening}% 499\idx{coarsening!restriction of DOF vectors}% 500\idx{restriction of DOF vectors} 501 502In an (Euler, e.g.) discretization of a time dependent problem, the 503term $(\uh^\mathrm{old}, \vphi_i)_{L^2(\Omega)}$ appears on the right 504hand side of the discrete system, where $\uh^\mathrm{old}$ is the 505solution from the last time step. Such an expression can be calculated 506exactly, if the grid does not change from one time step to another. 507Assuming that the finite element spaces are nested, it is also 508possible to calculate this expression exactly when the grid was 509refined, since $\uh^\mathrm{old}$ belongs to the fine grid finite 510element space also. Usually, during coarsening information is lost, 511since we can not represent $\uh^\mathrm{old}$ exactly on a coarser 512grid. But we can calculate $(\uh^\mathrm{old}, \psi_i)_{L^2(\Omega)}$ 513exactly on the fine grid; using the representation of the 514coarse grid basis functions $\vphi_i$ by the fine grid functions 515$\psi_i$, we can transform data during coarsening such that 516$(\uh^\mathrm{old}, \vphi_i)_{L^2(\Omega)}$ is calculated exactly 517for the coarse grid functions too. 518 519More general, assume that the finite element spaces are nested and 520that we can evaluate a linear functional $f$ exactly for all basis 521functions of the fine grid. Knowing the values 522$\vec{f}_\psi = (\langle f,\,\psi^1\rangle, \dots, 523\langle f,\,\psi^k\rangle)^t$ 524for the fine grid functions, we obtain with matrix $A$ 525from \mathref{E:a_matrix} for the values 526$\vec{f}_\vphi = (\langle f,\,\vphi^1\rangle, \dots, \langle 527f,\,\vphi^m\rangle)^t $ on the coarse grid 528\[ 529\vec{f}_\vphi = A \vec{f}_\psi 530\] 531since 532\[ 533 \langle f, \, \vphi^i \rangle 534 = 535 \langle f,\,\sum_{j = 1}^k a_{ij}\psi^j \rangle 536 = 537 \sum_{j = 1}^k a_{ij} \langle f,\,\psi^j \rangle 538\] 539holds (here we used the fact, that $\I^S_{S_0,S_1} = id$ on $X_S$ 540since the spaces are nested). 541 542Thus, given a functional $f$ which we can evaluate exactly 543for all basis functions of a grid $\tilde \tria$ and its refinements, 544we can also calculate $\langle f, \, \vphi^i \rangle$ exactly for all basis 545functions $\vphi^i$ of a grid $\tria$ obtained by 546refinement and coarsening of $\tilde \tria$ in the following way: 547First refine all elements of the grid that have to be refined; 548calculate $\langle f, \, \vphi \rangle$ for all basis functions $\vphi$ 549of this intermediate grid; in the last step coarsen all elements 550that may be coarsened and restrict this vector during each coarsening 551step as described above. 552 553In \ALBERTA the assemblage of the discrete system inside the adaptive 554method can be split into three steps: one initialization step before 555refinement, the second step between refinement and coarsening, and the 556last, and usually most important, step after coarsening, when all grid 557modifications are completed, see \secref{man:S:adapt_stat_in_ALBERTA}. 558The second assemblage step can be used for an exact computation of a 559functional $\langle f, \, \vphi \rangle$ as described above. 560 561 562\subsubsection{Interpolation during coarsening}% 563\idx{coarsening!interpolation of DOF vectors}% 564\idx{interpolation of DOF vectors} 565 566Finally, we can interpolate a finite element function 567during coarsening. The matrix for transforming the coefficient vector 568$\vec{u}_\psi = (u^1_\psi, \dots, u^k_\psi)^t$ of a fine grid function $\uh$ 569to the coefficient vector $\vec{u}_\vphi = (u^1_\vphi, \dots, u^m_\vphi)^t$ 570of the interpolant on the coarse grid is given by matrix $B$ defined 571in \mathref{E:b_matrix}: 572\begin{align*} 573\I^{S_0,S_1}_S \uh & = \I^{S_0,S_1}_S \sum_{j = 1}^k u^j_\psi\, \psi^j 574 = \sum_{j = 1}^k u^j_\psi\, \I^{S_0,S_1}_S \psi^j\\ 575& = \sum_{j = 1}^k u^j_\psi\, \sum_{i = 1}^m b_{ij}\, \vphi^i 576 = \sum_{i = 1}^m \left(\sum_{j = 1}^k b_{ij}\, u^j_\psi\right) \vphi^i. 577\end{align*} 578Thus we have the following equation for the coefficient vector of 579the interpolant of $\uh$: 580\[ 581\vec{u}_\vphi = B\, \vec{u}_\psi. 582\] 583In contrast to the interpolation during refinement and the above 584described transformation of a linear functional, information is 585lost during an interpolation to the coarser grid. 586 587\begin{example}[Lagrange elements] 588Lagrange finite elements are connected to Lagrange nodes $x^i$. For 589linear elements, these nodes are the vertices of the triangulation, 590and for quadratic elements, the vertices and the edge--midpoints. 591The Lagrange basis functions $\big\{\vphi^i\big\}$ satisfy 592\[ 593\vphi_i(x_j) = \delta_{ij} \qquad \mbox{for } i,j = 1,\dots, \dim X_h. 594\] 595Consider the situation of a simplex $S$ with children $S_0$, $S_1$. 596Let $\big\{\vphi^i\big\}_{i = 1,\dots,m}$ be the Lagrange basis functions 597of $X_S$ with Lagrange nodes $\big\{x^i_\vphi\big\}_{i = 1,\dots,m}$ on $S$ 598and $\big\{\psi^j\}_{j = 1,\dots,k}$ be the Lagrange basis functions of 599$X_{S_0 \cup S_1}$ with Lagrange nodes 600$\big\{x^j_\psi\}_{j = 1,\dots,k}$ on $S_0 \cup S_1$. The Matrix $A$ is then 601given by 602\[ 603 a_{ij} = \vphi^i(x^j_\psi), \qquad i = 1,\dots,m,\; j = 1,\dots,k 604\] 605and matrix $B$ is given by 606\[ 607 b_{ij} = \psi^j(x^i_\vphi), \qquad i = 1,\dots,m,\; j = 1,\dots,k. 608\] 609\end{example} 610\idx{interpolation and restriction of DOF vectors|)} 611 612\subsection{Discretization of 2nd order problems}% 613\label{S:Dis2ndorder}% 614\idx{finite element discretization|(} 615 616In this section we describe the assembling of the discrete 617system in detail. 618We consider the following second order differential equation in 619divergence form: 620\begin{subequations}\label{E:strong} 621\begin{alignat}{2} 622L u := -\nabla \cdot A \nabla u + b \cdot \nabla u + c\, u 623 &= f \qquad & &\mbox{in } \Omega,\label{E:strong.a}\\ 624 u &= g & &\mbox{on } \GD,\label{E:strong.b}\\ 625\nO \cdot A \nabla u & = 0 &&\mbox{on } \GN\label{E:strong.c}, 626\end{alignat} 627\end{subequations} 628where $A \in L^\infty(\Omega;\R^{n \times n})$, 629$b \in L^\infty(\Omega;\R^n)$, $c \in L^\infty(\Omega)$, and 630$f \in L^2(\Omega)$. 631By $\GD \subset \partial \Omega$ (with $|\GD| \ne 0$) we denote the 632Dirichlet boundary and we assume that the Dirichlet boundary values 633$g \colon \GD \to \R$ have an extension to some function 634$g \in H^1(\Omega)$. 635 636By $\GN = \partial \Omega \backslash \GD$ we denote the 637Neumann boundary, and by $\nO$ we denote the outer unit normal vector 638on $\partial \Omega$. The boundary condition \mathref{E:strong.c} 639is a so called \emph{natural} Neumann condition. 640 641Equations \mathref{E:strong} describe not only a simple model problem. 642The same kind of equations result from a linearization of nonlinear 643elliptic problems (for example by a Newton method) as well as from a time 644discretization scheme for \hbox{(non--)} linear parabolic problems. 645 646Setting 647\begin{equation*}\label{E:X-X0} 648X = H^{1}(\Omega) \qquad \mbox{and} \qquad 649\Xc = \left\{v \in H^{1}(\Omega);\, v = 0 \mbox{ on } \GD\right\} 650\end{equation*} 651this equation has the following weak formulation: We are looking 652for a solution $u \in X$, such that $u \in g+\Xc$ and 653\begin{equation}\label{E:weak} 654\int_\Omega (\nabla \varphi(x)) \cdot A(x) \nabla u(x) 655+ \varphi(x)\, b(x) \cdot \nabla u(x) 656+ c(x)\,\varphi(x) \, u(x) \, dx = 657\int_\Omega f(x)\, \varphi(x)\, dx 658\end{equation} 659for all $\vphi \in \Xc$ 660 661Denoting by $\Xc^*$ the dual space of $\Xc$ we 662identify the differential operator $L$ with the linear operator 663$L \in \Lin(X,\Xc^*)$ defined by 664\begin{equation*}\label{E:Lu} 665\ldual{L v}{\varphi}{\Xc^* \times \Xc} := 666\int_\Omega \nabla \varphi \cdot A \nabla v 667 + \int_\Omega \varphi \, b \cdot \nabla v 668 + \int_\Omega c\,\varphi \, v 669\qquad \mbox{for all } v, \varphi \in \Xc 670\end{equation*} 671and the right hand side $f$ with 672the linear functional $f \in \Xc^*$ defined by 673\begin{equation*} 674\ldual{F}{\varphi}{\Xc^* \times \Xc} := 675\int_\Omega f\, \varphi \qquad \mbox{for all } \varphi \in \Xc. 676\end{equation*} 677 678With these identifications we use the following reformulation of 679\mathref{E:weak}: Find $u \in X$ such that 680\begin{equation}\label{E:weak2} 681u \in g + \Xc: \qquad L\, u = f \qquad\mbox{in }\Xc^* 682\end{equation} 683holds. 684 685Suitable assumptions on the coefficients imply that $L$ is elliptic, 686i.e.\ there is a constant $C = C_{A,b,c,\Omega}$ such that 687\[ 688 \ldual{L\,\vphi}{\vphi}{\Xc^* \times \Xc} \ge C\,\lnorm{\vphi}{X}^2 689 \qquad\mbox{for all }\vphi\in\Xc. 690\] 691The existence of a unique solution $u \in X$ of \mathref{E:weak2} 692is then a direct consequence of the Lax--Milgram--Theorem. 693 694We consider a finite dimensional subspace 695$X_h \subset X$ for the discretization of \mathref{E:weak2} with 696$N = \dim\ X_h$. We set $\Xc_h = X_h \cap \Xc$ with 697$\Nc = \dim\ \Xc_h$. Let $\gh \in X_h$ be an approximation of $g \in X$. 698A discrete solution of \mathref{E:weak2} is then given by: Find 699$\uh \in X_h$ such that 700% 701\begin{equation}\label{E:discrete} 702\uh \in \gh + \Xc_h: \qquad L\, \uh = f \qquad\mbox{in } \Xc_h^*, 703\end{equation} 704i.e. 705\[ 706\uh \in \gh + \Xc_h: \qquad 707\ldual{L \uh}{\ph}{\Xc_h^* \times \Xc_h} 708 = 709\ldual{f}{\ph}{\Xc_h^* \times \Xc_h} \quad 710\mbox{for all } \ph \in \Xc_h 711\] 712holds. If $L$ is elliptic, we have a unique discrete solution 713$\uh \in X_h$ of \mathref{E:discrete}, again using the 714Lax--Milgram--Theorem. 715 716\def\vv{\vec{v}} 717\def\vu{\vec{u}} 718\def\vL{\vec{L}} 719\def\vf{\vec{f}} 720 721Choose a basis $\left\{\vphi_1,\dots,\vphi_{N\vphantom{\Nc}}\right\}$ 722of $X_h$ such that $\left\{\vphi_1,\dots,\vphi_{\Nc}\right\}$ is a 723basis of $\Xc_h$. For a function $\vh \in X_h$ we denote by $\vv = 724(v_1,\dots,v_N)$ the coefficient vector of $\vh$ with respect to the 725basis $\left\{\vphi_1,\dots,\vphi_{N\vphantom{\Nc}}\right\}$, i.e. 726\[ 727\vh = \sum_{j = 1}^N v_j \vphi_j. 728\] 729Using \mathref{E:discrete} with test functions $\vphi_i$, $i = 7301,\dots,\Nc$, we get the following $N$ equations for the coefficient 731vector $\vu = (u_1,\dots,u_N)$ of $\uh$: 732%\begin{subequations*}\label{E:discrete.system} 733\begin{alignat*}{2} 734\sum\limits_{j = 1}^N u_j \,\ldual{L \vphi_j}{\vphi_i}{\Xc_h^* \times \Xc_h} 735&= \ldual{f}{\vphi_i}{\Xc_h^* \times \Xc_h} 736\qquad & &\mbox{for } i = 1,\dots,\Nc, \\ 737u_i & = g_i & &\mbox{for } i = \Nc+1,\dots,N. 738\end{alignat*} 739%\end{subequations*} 740Defining the \emph{system matrix} $\vec{L}$ by 741\begin{equation}\label{E:matrix} 742\vL := 743\left[ 744\begin{matrix} 745\ldual{L \vphi_1}{\vphi_1}{} & \dots & \ldual{L \vphi_{\Nc}}{\vphi_1}{} & 746\ldual{L \vphi_{\Nc+1}}{\vphi_1}{} & \dots & \ldual{L \vphi_{N\vphantom{\Nc}}}{\vphi_1}{}\\ 747\vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ 748\ldual{L \vphi_1}{\vphi_{\Nc}}{} & \dots & \ldual{L \vphi_{\Nc}}{\vphi_{\Nc}}{} & 749\ldual{L \vphi_{\Nc+1}}{\vphi_{\Nc}}{} & \dots & \ldual{L \vphi_{N\vphantom{\Nc}}}{\vphi_{\Nc}}{}\\ 7500 & \dots & 0 & \multicolumn{3}{c}{1 \hfil 0 \hfil \dots \hfil 0}\\ 7510 & \dots & 0 & \multicolumn{3}{c}{0 \hfil 1 \hfil \dots \hfil 0}\\ 752\vdots & \ddots &0&\multicolumn{3}{c}{0 \hfil 0 \hfil \ddots \hfil \vdots\,}\\ 7530 & \dots & 0 & \multicolumn{3}{c}{0 \hfil 0 \hfil \dots \hfil 1}\\ 754\end{matrix} 755\right] 756\end{equation} 757and the \emph{right hand side vector} or \emph{load vector} $\vf$ by 758\begin{equation}\label{E:Fright} 759\vf := 760\left[ 761\begin{matrix} 762\ldual{f}{\vphi_1}{}\\ 763\vdots\\ 764\ldual{f}{\vphi_{\Nc}}{}\\ 765g_{\Nc+1}\\ 766\vdots\\ 767g_{N\vphantom{\Nc}} 768\end{matrix} 769\right], 770\end{equation} 771we can write the discrete equations as the linear $N \times N$ system 772\begin{equation}\label{E:LuF} 773\vL \, \vu = \vf, 774\end{equation} 775which has to be assembled and solved numerically. 776 777\subsection{Discretization of coupled vector valued problems}% 778\label{S:DisCoupled}% 779\idx{coupled vector valued problems} 780 781This section describes the discretization and assembling of coupled 782vector valued problems. Consider the following artificial coupled Poisson 783problem: 784 785Let $C\in\R^{n\times n}$ be a regular coupling matrix, $f \in 786L^2(\Omega;\R^n)$ a given right hand side, and 787$g\colon\partial\Omega\to\R^n$ suitable boundary values. Find 788$u:\Omega\to\R^n$ with 789 790\begin{subequations}\label{E:coupled_example} 791\begin{alignat}{2} 792 -\sum_{\nu=1}^n C_{\mu\nu} \Delta u_\nu &= f_\mu \qquad\text{in }\Omega\text{ 793 for }\mu=1,\dots,n\\ u &= g \qquad\text{on }\partial\Omega, 794\end{alignat} 795\end{subequations} 796 797By a left multiplication with $C^{-1}$ this problem decouples 798into a set of independent scalar Poisson problems, for which we could 799apply the same existence and uniqueness theory as above. However, we 800will refrain from doing this in order to illustrate the concepts of 801this section. Generally, the weak form of a coupled system of linear second 802order equations can be written as follows: 803 804Define vector valued spaces $X=H^1(\Omega;\R^n)$, 805$\Xc=H_0^1(\Omega;\R^n)$. Find a solution $u \in X$, such that $u \in 806g+\Xc$ and 807\begin{equation}\label{E:coupled_weak} 808 \ldual{L u}{\vphi}{\Xc^* \times \Xc} := 809 \sum_{\mu,\nu=1}^n\int_\Omega \nabla \vphi_\mu \cdot A^{\mu\nu} \nabla u_\nu 810 + \vphi_\mu\, b^{\mu\nu} \cdot \nabla u_\nu 811 + c^{\mu\nu}\,\vphi_\mu \, u_\nu \, dx = 812 \int_\Omega f \cdot \vphi\, dx 813\end{equation} 814for all $\vphi \in \Xc$. 815 816To obtain the weak form of problem \mathref{E:coupled_example} for 817example, we would set $b:=0$, $c:=0$ and 818$A_{ij}^{\mu\nu}:=\delta_{ij}C_{\mu\nu}$. The next step is to derive a 819suitable linear system for a discretization as in the last section. 820 821As mentioned in \secref{S:FES}, basis functions are always 822scalar-valued. To gain a vector valued finite element space $X_h$, we 823use vector valued coefficients. Choose a set of scalar basis functions 824$\left\{\vphi_1,\dots,\vphi_{N\vphantom{\Nc}}\right\}$ as above. For a 825function $\vh \in X_h$ we denote by $\vv = (v_1,\dots,v_N)$ the 826coefficient vector of $\vh$. Each $v_j$ is now itself a vector 827$v_j=(v_{j\mu})_{\mu=1}^n$. Thus, we have the following decomposition: 828\[ 829\vh = \sum_{j = 1}^N v_j \vphi_j. 830\] 831Define $\vphi_j^\mu:=(\delta_{\mu\nu}\vphi_j)_{\nu=1}^n$. The discrete problem 832can now be written as a set of linear equations for the coefficients 833$u_{j\mu}$: 834\begin{alignat*}{2} 835 \sum_{j = 1}^N \sum_{\mu=1}^n u_{j\mu} \,\ldual{L 836 \vphi_j^\mu}{\vphi_i^\nu}{\Xc_h^* \times \Xc_h} &= 837 \ldual{f}{\vphi_i^\nu}{\Xc_h^* \times \Xc_h} \qquad & &\text{for } i = 838 1,\dots,\Nc; \nu = 1,\dots,n,\\ u_{i\nu} & = g_{i\nu} & &\text{for } i = 839 \Nc+1,\dots,N; \nu=1,\dots,n. 840\end{alignat*} 841The corresponding system matrix $\vL$ is defined by 842\begin{equation}\label{E:coupled_matrix} 843\vL := 844\left[ 845\begin{matrix} 846 \vL^{11} & \dots & \vL^{1\Nc} & \vL^{1,\Nc+1} & \dots & \vL^{1N}\\ 847 \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ 848 \vL^{\Nc1} & \dots & \vL^{\Nc\Nc} & \vL^{\Nc,\Nc+1} & \dots & \vL^{\Nc N}\\ 849 0 & \dots & 0 & \multicolumn{3}{c}{\vI \hfil 0 \hfil \dots \hfil 0}\\ 850 0 & \dots & 0 & \multicolumn{3}{c}{0 \hfil \vI \hfil \dots \hfil 0}\\ 851 \vdots & \ddots &0&\multicolumn{3}{c}{0 \hfil 0 \hfil \ddots \hfil\vdots\,}\\ 852 0 & \dots & 0 & \multicolumn{3}{c}{0 \hfil 0 \hfil \dots \hfil \vI}\\ 853 \end{matrix} 854\right] 855\end{equation} 856with $\vI\in\R^{n\times n}$ an identity matrix and 857\[ 858\vL^{ij} := (\ldual{L\vphi_j^\mu}{\vphi_i^\nu}{\Xc_h^* \times \Xc_h})_{\mu,\nu=1}^n. 859\] 860The load vector $\vf$ is defined by 861\begin{equation}\label{E:coupled_Fright} 862\vf := 863\left[ 864\begin{matrix} 865\ldual{f}{\vphi_1^1}{}\\ 866\vdots\\ 867\ldual{f}{\vphi_1^n}{}\\ 868\vdots\\ 869\ldual{f}{\vphi_{\Nc}^1}{}\\ 870\vdots\\ 871\ldual{f}{\vphi_{\Nc}^n}{}\\ 872g_{\Nc+1,1}\\ 873\vdots\\ 874g_{Nd} 875\end{matrix} 876\right]. 877\end{equation} 878The problem can now be written as the linear $Nd \times Nd$ system 879\begin{equation}\label{E:coupled_LuF} 880 \vL \, \vu = \vf, 881\end{equation} 882which has to be assembled and solved. The organization of vectors and matrices 883using small $n$-size blocks as components was chosen with the goal of 884efficient cache usage during matrix-vector multiplication. 885 886\subsection{Numerical quadrature}% 887\label{S:quadrature}% 888\idx{numerical quadrature} 889 890For the assemblage of the system matrix and right hand side 891vector of the linear system \mathref{E:LuF}, we have to compute 892integrals, for example 893\[ 894\int_\Omega f(x) \vphi_i(x)\,dx. 895\] 896For general data $A$, $b$, $c$, and $f$, these integrals can not be 897calculated exactly. Quadrature formulas have to be used in order to 898calculate the integrals approximately. Numerical integration in finite 899element methods is done by looping over all grid elements and using a 900quadrature formula on each element. 901 902\begin{definition}[Numerical quadrature]\label{D:numer-quad} 903\idx{numerical quadrature} 904A \emph{numerical quadrature} $\hat Q$ on $\Shat$ is a set 905$\{(w_k,\lambda_k) \in \R \times \R^{d+1}; k = 0, \dots, n_Q-1\}$ 906of weights $w_k$ and quadrature points $\lambda_k \in \bar S$ 907(i.e.\ given in barycentric coordinates) such that 908\[ 909 \int_{\Shat} f(\hat x)\, d\hat{x} \approx \hat Q(f) := 910 \sum_{k = 0}^{n_Q-1} w_k f(\hat{x}(\lambda_k)). 911\] 912It is called \emph{exact of degree $p$} for some $p \in \N$ if 913\[ 914 \int_{\Shat} q(\hat x)\, d\hat{x} = \hat Q(q) 915 \qquad\mbox{for all }q \in \P_p(\Shat). 916\] 917It is called \emph{stable} if 918\[ 919 w_k > 0 \qquad\mbox{for all } k = 0, \dots, n_Q-1. 920\] 921\end{definition} 922 923\begin{remark}\label{R:numerical_int} A given numerical quadrature $\hat Q$ on 924$\Shat$ defines for each element $S$ a numerical quadrature $Q_S$. 925Using the transformation rule we define 926$Q_S$ on an element $S$ which is parameterized by 927$F_S \colon\, \Shat \to S$ and a function $f : S \to \R$: 928\begin{equation*}\label{E:numer-quad-S} 929\int_{S} f(x)\, dx \approx Q_S(f) := 930\hat Q((f \circ F_S) |\det DF_S|) = 931\sum_{k = 0}^{n_Q-1} w_k f(x(\lambda_k)) |\det DF_S(\hat x(\lambda_k)|. 932\end{equation*} 933For a simplex $S$ this results in 934\begin{equation*}\label{E:numer-quad-sim} 935Q_S(f) = d!\, |S| \sum_{k = 0}^{n_Q-1} w_k f(x(\lambda_k)). 936\end{equation*} 937\end{remark} 938 939\subsection{Finite element discretization of 2nd order problems}% 940\label{S:FE-dis-2nd}% 941\idx{assemblage of discrete system|(} 942 943Let $\Pbar$ be a finite dimensional function space on $\Sbar$ with 944basis $\{\pbar^1,\dots,\pbar^m\}$. For a conforming finite element 945discretization of \mathref{E:weak} we use the finite element space 946$X_h = X_h(\tria,\Pbar,X)$. For this space $\Xc_h$ is 947given by $X_h(\tria,\Pbar,\Xc)$. 948 949\idx{assemblage of discrete system!load vector} 950By the relation \mathref{E:global-lokal.a} for 951global and local basis functions, we obtain for the $j$th component of 952the right hand side vector $\vf$ 953\begin{align*} 954\ldual{f}{\vphi_j}{} &= 955\int_\Omega f(x)\, \vphi_j(x)\,dx 956= 957\sum_{S \in \tria} \int_S f(x)\, \vphi_j(x)\,dx 958= 959\sum_{\atop{S \in \tria}{S \subset \supp(\vphi_j)}} 960 \int_S f(x) \pbar^{i_S(j)}(\lambda(x))\,dx\\ 961&= 962\sum_{\atop{S \in \tria}{S \subset \supp(\vphi_j)}} 963 \int_{\Shat} f(F_S(\xhat)) \pbar^{i_S(j)}(\lambda(\xhat)) 964 |\det DF_S(\xhat)|\,d\xhat, 965\end{align*} 966where $S$ is parameterized by $F_S\colon\Shat\to S$. 967The above sum is reduced to a sum over all $S \subset \supp(\vphi_j)$ 968which are only few elements due to the small support of $\vphi_j$. 969 970The right hand side vector can be assembled in the following way: 971First, the right hand side vector is initialized with zeros. For each 972element $S$ of $\tria$ we calculate the \emph{element load vector} 973$\vf_S = (f_S^1,\dots,f_S^m)^t$, where 974\begin{equation}\label{E:int_S-f-phi} 975f_S^i = \int_{\Shat} f(F_S(\xhat)) \pbar^{i}(\lambda(\xhat)) 976 |\det DF_S(\xhat)|\,d\xhat, \qquad i = 1,\dots,m. 977\end{equation} 978Denoting again by $j_S: \{1,\dots,m\} \to J_S$ the function which connects 979the local DOFs with the global DOFs (defined in \mathref{E:global-lokal.b}), 980the values $f_S^i$ are then added to the $j_S(i)$th 981component of the right hand side vector $\vf$, $i = 1,\dots,m$. 982 983For general $f$, the integrals in \mathref{E:int_S-f-phi} can not be 984calculated exactly and we have to use a quadrature formula for the 985approximation of the integrals (compare Section~\ref{S:quadrature}). 986Given a numerical quadrature $\hat Q$ on $\Shat$ we approximate 987\begin{align}\label{EX:quad_right_param} 988f_S^i &\approx \hat Q\left( 989 (f \circ F_S)\,(\pbar^{i}\circ \lambda) |\det DF_S| \right) 990 = \sum_{k = 0}^{n_Q-1} w_k f(x(\lambda_k))\,\pbar^{i}(\lambda_k) 991 |\det DF_S(\hat x(\lambda_k)|. 992\end{align} 993For a simplex $S$ this is simplified to 994\begin{equation*}\label{E:quad_right} 995f_S^i \approx 996d!\, |S| \,\sum_{k = 0}^{n_Q-1} w_k f(x(\lambda_k))\,\pbar^{i}(\lambda_k). 997\end{equation*} 998 999In \ALBERTA, information about values of basis functions and its 1000derivatives as well as information about the connection of global and 1001local DOFs (i.e.\ information about $j_S$) is stored in special data 1002structures for local basis functions (compare 1003Section~\ref{man:S:basfct_impl}). By such information, the element 1004load vector can be assembled by a general routine if a function for 1005the evaluation of the right hand side is supplied. For parametric 1006elements, a function for evaluating $|\det DF_S|$ is additionally 1007needed. The assemblage into the global load vector $\vf$ can again be 1008performed automatically, using information about the connection of 1009global and local DOFs. 1010 1011\idx{assemblage of discrete system!system matrix} 1012The calculation of the system matrix is also done element--wise. 1013Additionally, we have to handle derivatives of basis functions. 1014Looking first at the second order term we obtain by the chain 1015rule \mathref{E:grad_phi} and the relation \mathref{E:global-lokal} for 1016global and local basis functions 1017\begin{align*} 1018\int_S \nabla \vphi_i(x) \cdot A(x) \nabla \vphi_j(x)\,dx &= 1019\int_S \nabla (\pbar^{i_S(i)} \circ \lambda)(x) \cdot A(x) 1020 \nabla (\pbar^{i_S(j)} \circ \lambda)(x) \,dx\\ 1021&= 1022\int_S 1023 \nablal \pbar^{i_S(i)}(\lambda(x)) \cdot (\Lambda(x)\, A(x)\, \Lambda^t(x)) 1024 \nablal \pbar^{i_S(j)}(\lambda(x))\,dx, 1025\end{align*} 1026where $\Lambda$ is the Jacobian of the barycentric coordinates $\lambda$ 1027on $S$. In the same manner we obtain for the first and zero order terms 1028\[ 1029 \int_S \vphi_i(x) \, b(x) \cdot \nabla \vphi_j(x) \,dx = 1030\int_S \pbar^{i_S(i)}(\lambda(x))\,(\Lambda(x)\, b(x)) \cdot 1031 \nablal \pbar^{i_S(j)}(\lambda(x)) \,dx 1032\] 1033and 1034\[ 1035 \int_S c(x)\,\vphi_i(x) \, \vphi_j(x)\,dx = 1036 \int_S c(x)\, \pbar^{i_S(i)}(\lambda(x)) \, \pbar^{i_S(j)}(\lambda(x))\,dx. 1037\] 1038Using on $S$ the abbreviations 1039%\begin{subequations*}\label{E:bar_A_b_c} 1040\begin{align*} 1041\bar A(\lambda) &:= 1042\left(\bar a_{kl}(\lambda)\right)_{k,l = 0,\ldots,d} 1043 := |\det DF_S(\xhat(\lambda))| \, 1044\Lambda(x(\lambda)) \, A(x(\lambda)) \, \Lambda^t(x(\lambda)),\\ 1045\bar b(\lambda) & := 1046\left(\bar b_{l}(\lambda)\right)_{l = 0,\ldots,d} 1047 := |\det DF_S(\xhat(\lambda))| \, 1048\Lambda(x(\lambda)) \, b(x(\lambda)), \quad \mbox{and}\\ 1049\bar c(\lambda) & := |\det DF_S(\xhat(\lambda))| \, c(x(\lambda)) 1050\end{align*} 1051%\end{subequations*} 1052and transforming the integrals to the reference simplex, we can write 1053the \emph{element matrix} $\vL_S = (L_S^{ij})_{i,j=1,\dots,m}$ as 1054\begin{align}\label{E:L-phii-phij} 1055L_S^{ij} &= 1056\int_{\Shat} 1057 \nablal \pbar^{i}(\lambda(\xhat)) \cdot \bar A(\lambda(\xhat))\, 1058 \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat 1059+ 1060\int_{\Shat} \pbar^{i}(\lambda(\xhat))\; 1061 \bar b(\lambda(\xhat)) \cdot 1062 \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat \notag\\ 1063&\qquad + 1064\int_{\Shat} \bar c(\lambda(\xhat))\, \pbar^{i}(\lambda(\xhat)) \, 1065 \pbar^{j}(\lambda(\xhat))\,d\xhat, 1066\end{align} 1067or writing the matrix--vector and vector--vector products explicitly 1068\begin{align*}\label{E:L-phii-phij'} 1069L_{S}^{ij} &= 1070\sum_{k,l = 0}^d \int_{\Shat} 1071\bar a_{kl}(\lambda(\xhat)) \, \pbar_{,\lambda_k}^i(\lambda(\xhat)) \, 1072 \pbar_{,\lambda_l}^j(\lambda(\xhat)) \, d\xhat 1073+ \sum_{l = 0}^d \int_{\Shat} 1074\bar b_{l}(\lambda(\xhat)) \, \pbar^i(\lambda(\xhat)) \, 1075 \pbar_{,\lambda_l}^j(\lambda(\xhat))\, d\xhat 1076%\notag 1077\\ 1078%\tag{\ref{E:L-phii-phij}'} 1079&\qquad + \int_{\Shat} 1080\bar c(\lambda(\xhat)) \, \pbar^i(\lambda(\xhat))\,\pbar^j(\lambda(\xhat))\, 1081 d\xhat, 1082\end{align*} 1083$i,j = 1,\dots,m$. Using quadrature formulas $\hat Q_2$, $\hat Q_1$, and 1084$\hat Q_0$ for the second, first and zero order term we approximate 1085the element matrix 1086\begin{equation*}\label{E:quad_L} 1087L_S^{ij} \approx 1088\hat Q_2\left(\sum_{k,l = 0}^d (\bar a_{kl} \pbar^i_{,\lambda_k} \, 1089 \pbar^j_{,\lambda_l}) \circ\lambda\right) 1090 + 1091\hat Q_1\left(\sum_{l = 0}^d ( 1092\bar b_{l} \, \pbar^i \, \pbar^j_{,\lambda_l})\circ\lambda\right) 1093 + \hat Q_0\Big( 1094 (\bar c \, \pbar^i \,\pbar^j)\circ\lambda\Big), 1095\end{equation*} 1096$i,j=1,\dots,m$. Having access to functions for the evaluation of 1097\begin{equation*}\label{E:LALt-Lb-c} 1098\bar a_{kl}(\lambda_q), \quad \bar b_{l}(\lambda_q),\quad 1099\bar c(\lambda_q) 1100\end{equation*} 1101at all quadrature points $\lambda_q$ on $S$, $\vec{L}_S$ can 1102be computed by some general routine. The assemblage into the 1103system matrix can also be done automatically (compare the assemblage 1104of the load vector). 1105 1106\begin{remark} 1107Due to the small support of the global basis 1108function, the system matrix is a sparse matrix, i.e.\ the maximal 1109number of entries in all matrix rows is much smaller than the 1110size of the matrix. Special data structures are needed for an 1111efficient storage of sparse matrices and they are described in 1112Section~\ref{man:S:DOF_MATRIX}. 1113\end{remark} 1114 1115\begin{remark} 1116The calculation of the gradient of the barycentric coordinates 1117$\Lambda(x(\lambda))$ usually involves the 1118determinant of the parameterization's Jacobian 1119$|\det DF_S(\xhat(\lambda))|$. Thus, a calculation of 1120$|\det DF_S(\xhat(\lambda))| \,\Lambda(x(\lambda)) \, A(x(\lambda)) 1121\, \Lambda^t(x(\lambda))$ may be much faster than the calculation of 1122$\Lambda(x(\lambda)) \, A(x(\lambda)) \, \Lambda^t(x(\lambda))$ only; 1123the same holds for the first order term. 1124\end{remark} 1125 1126Assuming that the coefficients $A$, $b$, and $c$ are piecewise 1127constant on a non--parametric triangulation, $\bar A(\lambda)$, 1128$\bar b(\lambda)$, and $\bar c(\lambda)$ are constant on each simplex $S$ 1129and thus simplified to 1130\begin{equation*}\label{E:bar_A_b_cS} 1131\bar A_S = \left(\bar a_{kl}\right)_{k,l = 0,\ldots,d} 1132 = d!|S| \, \Lambda \, A_{|S} \, \Lambda^t, \quad 1133\bar b_S = \left(\bar b_{l}\right)_{l = 0,\ldots,d} 1134 = d!|S| \,\Lambda \, b_{|S}, \quad 1135\bar c_S = d!|S| \, c_{|S}. 1136\end{equation*} 1137For the approximation of the element matrix by quadrature we then 1138obtain 1139\begin{equation}\label{E:quad_LS} 1140\vec{L}_S^{ij} \approx 1141\sum_{k,l = 0}^d \bar a_{kl} 1142\hat Q_2\Big((\pbar^i_{,\lambda_k} \, \pbar^j_{,\lambda_l})\circ\lambda\Big) 1143 + 1144\sum_{l = 0}^d \bar b_{l} \, 1145\hat Q_1\Big((\pbar^i \, \pbar^j_{,\lambda_l})\circ\lambda\Big) 1146 + \bar c_S\, \hat Q_0\Big( 1147 (\pbar^i \,\pbar^j)\circ\lambda\Big) 1148\end{equation} 1149$i,j=1,\dots,m$. Here, the numerical quadrature is only applied 1150for approximating integrals of the basis 1151functions on the standard simplex. Theses values can be computed 1152only once, and can then be used on each simplex $S$. This will speed 1153up the assembling of the system matrix. Additionally, for polynomial 1154basis functions we can use quadrature formulas which integrate these 1155integrals exactly. 1156\idx{assemblage of discrete system|)} 1157\idx{assemblage of discrete system!coupled case|(} 1158So far we have only considered the case of scalar problems. The 1159transition to (coupled) vector valued problems is straight forward and 1160simply involves two more indices. The entries of the element matrix 1161are now $d\times d$ matrices themselves: 1162\begin{align*} 1163L_{S,\mu\nu}^{ij} &:= 1164 \int_S \nabla \vphi_i(x) \cdot A^{\mu\nu}(x) \nabla 1165 \vphi_j(x) + \vphi_i(x)\, b^{\mu\nu}(x) \cdot \nabla \vphi_j(x) 1166 + c^{\mu\nu}(x)\,\vphi_i(x) \, \vphi_j(x) \, dx\\ 1167 &= 1168 \int_{\Shat} 1169 \nablal \pbar^{i}(\lambda(\xhat)) \cdot \bar A^{\mu\nu}(\lambda(\xhat))\, 1170 \nablal \pbar^{j}(\lambda(\xhat))\,d\xhat 1171 + 1172 \int_{\Shat} \pbar^{i}(\lambda(\xhat))\; 1173 \bar b^{\mu\nu}(\lambda(\xhat)) \cdot 1174 \nablal \pbar^{j}(\lambda(\xhat)) \,d\xhat \notag\\ 1175 &\qquad + 1176 \int_{\Shat} \bar c^{\mu\nu}(\lambda(\xhat))\, \pbar^{i}(\lambda(\xhat)) \, 1177 \pbar^{j}(\lambda(\xhat))\,d\xhat, 1178\end{align*} 1179with 1180\begin{align*} 1181 \bar A^{\mu\nu}(\lambda) &:= \left(\bar 1182 a_{kl}^{\mu\nu}(\lambda)\right)_{k,l = 0,\ldots,d} := |\det 1183 DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, 1184 A^{\mu\nu}(x(\lambda)) \, \Lambda^t(x(\lambda)),\\ 1185 \bar b^{\mu\nu}(\lambda) & := \left(\bar 1186 b_{l}^{\mu\nu}(\lambda)\right)_{l = 0,\ldots,d} := |\det 1187 DF_S(\xhat(\lambda))| \, \Lambda(x(\lambda)) \, b^{\mu\nu}(x(\lambda)), 1188 \quad\text{and}\\ 1189 \bar c^{\mu\nu}(\lambda) & := |\det 1190 DF_S(\xhat(\lambda))| \, c^{\mu\nu}(x(\lambda)) 1191\end{align*} 1192for $\mu,\nu=1,\dots,d$. The approximation of the integrals 1193using quadratures is done analogously to the scalar case. 1194 1195As a result, using information about values of basis functions and 1196their derivatives, and information about the connection of global and 1197local DOFs, the linear system can be assembled automatically by some 1198general routines. Only functions for the evaluation of given data have 1199to be provided for special applications. The general assemble 1200routines are described in Section~\ref{man:S:ass_tools}. 1201\idx{assemblage of discrete system!coupled case|)}% 1202\idx{finite element discretization|)} 1203 1204%%% Local Variables: 1205%%% mode: latex 1206%%% TeX-master: "alberta-book" 1207%%% End: 1208