## 4 Stochastic Finite Element Methods
### Deterministic Part: \(p\)-Version  
In order to accurately model and simulate the effects of changes in the domain, a method capable of accurate geometry description is needed. The \(p\)- and \(hp\)-finite element methods are ideally suited to this requirement. The paper by [PERSON] and [PERSON] [18] gives an accessible overview of the method. For a more detailed exposition, see [PERSON] [19], and for those familiar with engineering approach, the book by [PERSON] and [PERSON] [13] is recommended. The selection of shape functions is not unique. The [PERSON] and [PERSON] hierarchic shape functions are based on the integrated Legendre polynomials: For \(x\in[-1,1]\)  
\[\phi_{n}(\xi)=\sqrt{\frac{2n-1}{2}}\int_{-1}^{\xi}P_{n-1}(t)\,dt,\quad n=2,3,\ldots \tag{5}\]  
The normalising coefficients are chosen so that \(\int_{-1}^{1}\frac{d\phi(\xi)}{d\xi}\frac{d\phi(\xi)}{d\xi}\,d\xi=\delta_{ij},\quad i,j\geq 2.\) For a quadrilateral reference element over the domain \([-1,1]\times[-1,1]\), the shape functions are divided into three categories: nodal shape functions, side modes, and internal modes. There are four nodal shape functions, \(4(p-1)\) side modes, and \((p-1)(p-1)\) interior modes per quadrilateral with order \((p\geq 2)\). Taken alone, the nodal shapes define the standard four-node quadrilateral finite element. For the internal modes have the tensor product structure  
\[N_{i,j}^{0}(\xi,\eta)=\phi_{i}(\xi)\phi_{j}(\eta),\quad i=2,\ldots,p,\quad j=2,\ldots,p. \tag{6}\]  
The internal shape functions are often referred to as bubble-functions. The Legendre polynomials have the property \(P_{n}(-x)=(-1)^{n}P_{n}(x)\), which is inherited by the integrated Legendre polynomials. It follows that since in 2D all internal edges of the mesh are shared by two different elements, additional bookkeeping is required to ensure that each edge has the same global parameterisation in both elements. The option of increasing discretisation accuracy with the polynomial order makes it possible to have as large elements as possible, provided the possibly curved boundary segments are represented accurately. The linear blending function method of [PERSON] [20] is the standard choice for this purpose. In the general case, all sides of an element can be curved; however, in the examples considered in this study, only one edge is curved, as in Figure 4. If every side is parameterised:  
\[x=x_{i}(t),\;y=y_{i}(t),\quad-1\leq t\leq 1,\quad i=1,\ldots,\text{number of sides}, \tag{7}\]then using capital letters as coordinates of the corner points, \((X_{i},Y_{i})\), the mapping for the global \(x\)-coordinates of a quadrilateral can be written as  
\[ \begin{split} x&=\frac{1}{2}(1-\eta)x_{1}(\xi)+\frac{1 }{2}(1+\xi)x_{2}(\eta)+\frac{1}{2}(1+\eta)x_{3}(\xi)+\frac{1}{2}(1-\xi)x_{4}( \eta)\\ &-\frac{1}{4}(1-\xi)(1-\eta)X_{1}-\frac{1}{4}(1+\xi)(1+\eta)X_{2 }-\frac{1}{4}(1+\xi)(1+\eta)X_{3}\\ &-\frac{1}{4}(1-\xi)(1+\eta)X_{4},\end{split}  \tag{8}\]  
and symmetrically for the \(y\)-coordinate. If the side parametrisations represent straight edges, the mapping simplifies to the standard bilinear mapping of quadrilaterals. Consider the abstract problem setting with \(u\) defined on the standard piecewise polynomial finite element space on some discretisation \(\mathcal{T}\) of the computational domain \(\Omega\). Assuming that the exact solution \(u\in H^{1}_{0}(D)\) has finite energy, we arrive at the approximation problem: Find \(\hat{u}\in V\) such that  
\[a(\hat{u},v)=l(v)\ (=a(u,v))\quad(\forall v\in V), \tag{9}\]  
where \(a(\cdot,\cdot)\) and \(l(\cdot)\), are the bilinear form and the load potential, respectively. Additional degrees of freedom can be introduced by enriching the space \(V\). This is accomplished via introduction of an auxiliary subspace or \"error space\" \(W\subset H^{1}_{0}(D)\) such that \(V\cap W=\{0\}\).