In modified form taken from

Ralf Hartmann
Adaptive Finite Element Methods for the Compressible Euler Equations
PhD thesis, University of Heidelberg, 2002.

Higher order Boundary approximation

Introduction: In many numerical applications the domain $\Omega\subset\mathbb R^d$, $d\geq 2$ is not a polygonal domain but includes curved boundaries. For these cases the boundary cannot be represented exactly by the discretised boundary. Approximating the boundary by a piecewise linear boundary interpolation, i.e. by a polygonal boundary, may in some applications not be sufficient. In these cases a higher order boundary approximation, for example by piecewise quadratic or cubic boundary interpolation, must be employed. In the finite element framework this higher order boundary approximation is realized by mapping the reference element $\hat K$ to the element $K=\sigma_K(\hat K)$ in real space, whereas on cells K at the boundary, i.e. $\partial K\cap\Gamma=\emptyset$, the mappings $\sigma_K$ are given by polynomial functions of higher degree.

Elements with general mapping functions $\sigma_K$. We begin by first introducing some notation. Let $K\in T_h$ be a cell of the triangulation Th with $K=\sigma_K(\hat K)$, where $\sigma_K$is a smooth bijective mapping of the reference element (unit square) $\hat K$ to the element K in real space, see Figure 1.

  
Figure: Mapping $\sigma $ of reference element $\hat K$ to the element K in real space.
\includegraphics[scale=0.3]{figures/m1.ps}

In the following and for the sake of simplicity we suppress the letter K in the subscript and write $\sigma $ instead of $\sigma_K$.

Mapping functions of higher polynomial degree. A mapping function $\sigma $, that maps the reference element $\hat K$ to an arbitrary quadrilateral cell K with straight boundaries, can in general be represented by a bilinear function, i.e. by a Q1-mapping. For the case that the cell K includes curved boundaries it might be necessary to employ polynomial mapping functions of higher degree.

Given a degree p>0, a cell $K\in T_h$, and (p+1)d mapping support points $p_i\in K$, $i=0,\ldots,(p+1)^d-1$, we define a Qp-mapping $\sigma\in [Q_p]^d$ as follows

 \begin{displaymath}
\sigma(\hat x)=\sum_{i=0}^{(p+1)^d-1}p_i\phi_i(\hat x).
\end{displaymath} (1)

Here, $\phi_i$, $i=0,\ldots,(p+1)^d-1$ denote the Lagrange interpolation basis functions, that satisfy

\begin{displaymath}\phi_i(\hat p_j)=\delta_{ij}, \quad i,j=0,\ldots,(p+1)^d-1,
\end{displaymath}

where $\hat p_i$, $i=0,\ldots,(p+1)^d-1$ denote the Lagrange support points on the unit cell $\hat K$. The definition of $\sigma $(1) ensures that each of the unit support points $\hat p_i$ is mapped onto the corresponding mapping support points pi, i.e.

 \begin{displaymath}
\sigma(\hat p_i)=p_i, \quad i=0,\ldots,(p+1)^d-1
\end{displaymath} (2)

Analogous to Lagrange finite elements the unit Lagrange support points $\hat p_i$ are equidistantly distributed on $\hat K$ based on a tensor product mesh. In the following we only consider the two-dimensional case, d=2. For that case, Figure 2 shows the distributions of the unit support points $\hat p_i$, $i=0,\ldots,(p+1)^2-1$ for degrees $p=1,\ldots ,4$.
  
Figure: Unit support points $\hat p_i$, $0\leq i<(p+1)^2$, for degrees $p=1,\ldots ,4$.
\includegraphics*[width=2.750cm]{../../../phd/figures/boundary_approximation/mapping/unit_points_q1_n.eps} \includegraphics*[width=2.750cm]{../../../phd/figures/boundary_approximation/mapping/unit_points_q2_n.eps} \includegraphics*[width=2.750cm]{../../../phd/figures/boundary_approximation/mapping/unit_points_q3_n.eps} \includegraphics*[width=2.750cm]{../../../phd/figures/boundary_approximation/mapping/unit_points_q4_n.eps}

Let the ordering and numbering of the unit support points be as follows: first the corners, then the points on the edges and finally the inner support points, see also Figure 2. Thus the first 4p points are placed on the boundary $\partial\hat K$ of the reference cell, i.e.

\begin{displaymath}\hat p_k\in\partial\hat K, \quad k=0,\ldots,4p-1.
\end{displaymath}

According to (2) these points are mapped to the mapping support points pk, $k=0,\ldots,4p-1$ that are chosen to be placed on the boundary of the real cell in approximatively equal distances, i.e.

\begin{displaymath}p_k\in\partial K, \quad k=0,\ldots,4p-1.
\end{displaymath}

While the support points pk, $k=0,\ldots,4p-1$ on the boundary are given by the boundary description of the real cell K, the inner mapping support points

\begin{displaymath}p_i\in K\setminus\partial K, \quad i=4p,\ldots,(p+1)^2-1
\end{displaymath}

are not uniquely determined. Numerical tests show that it is not a trivial task to define the positions of the inner mapping support points appropriately. If they are not chosen appropriately the resulting mapping $\sigma $ for a cell K may degenerate, i.e. the mapping $\sigma $ for some cell K may not be bijective.

Computation of inner support points by smooth transformation. In the following we will define the positions of the inner mapping support points so that the mapping does - in all practical cases - not degenerate. To this end, we employ an approach for the mapping of the support points, that is in the style of the smooth transformations that is used to transform structured triangulations to match complex boundary discriptions. In the following, again for notational convenience, we consider only the two-dimensional case.

The smooth transformation mentioned above is based on solutions to the Laplace equation that is solved on the reference cell $\hat K$. Discrete boundary conditions are imposed that are given by the coordinates of the mapping support points pk, $k=0,\ldots,4p-1$, on the boundary of the cell K in real space.

To be more explicite we define a Laplace problem on the unit cell $\hat K$

 \begin{displaymath}
\begin{array}{rcll}
-\hat\Delta \sigma_l(\hat x)&=&0, \quad...
...l \hat K}(\hat x), \quad &\hat x\in \partial\hat K,
\end{array}\end{displaymath} (3)

for each component $\sigma_l$, l=1,2, of the Qp mapping $\sigma $. Here, the discrete boundary function $g\in [Q_p]^2$ is given by

 \begin{displaymath}
g_l(\hat x)=\sum_{i=0}^{4p-1}(p_i)_l\phi_i(\hat x), \quad l=1,\ldots,d
\end{displaymath} (4)

where (pi)l denotes the lth component of the support point pi, and $\phi_i$ the corresponding Lagrangian interpolation basis function. We recall that the numbering of the mapping support points involves $p_k\in\partial K$ for $k=0,\ldots,4p-1$. Substituting

 \begin{displaymath}
\tilde \sigma_l:=\sigma_l-g_l, \quad l=1,2,
\end{displaymath} (5)

into the Laplace problem (3) yields the zero boundary value problem,

 \begin{displaymath}
\begin{array}{rcll}
-\hat\Delta \tilde\sigma_l(\hat x)&=&\h...
...gma_l(\hat x)&=&0, \quad &\hat x\in \partial\hat K,
\end{array}\end{displaymath} (6)

that is equivalent to the following variational formulation

\begin{displaymath}\tilde\sigma_l\in H^1_0(\hat K):\qquad
(\hat\nabla \tilde\si...
...\hat\nabla \phi)_{\hat K} \quad \forall \phi\in H^1_0(\hat K). \end{displaymath}

Discretisation of this problem

\begin{displaymath}\tilde\sigma_l\in Q_p(\hat K):\qquad
(\hat\nabla \tilde\sigm...
...bla \phi_{4p+i})_{\hat K} \quad \forall i=0,\ldots, (p-1)^2-1, \end{displaymath}

and recalling definitions (1), (5) and (4) gives

 \begin{displaymath}
\sum_{j=1}^{(p-1)^2}S_{ij}(p_{4p+j})_l=-\sum_{k=0}^{4p-1}T_{ik}(p_k)_l, \quad i=0, \ldots, (p-1)^2-1,
\end{displaymath} (7)

with the matrices $S_{ij}\in\mathbb R^{(p-1)^2\times(p-1)^2}$ and $T_{ik}\in\mathbb R^{(p-1)^2\times 4p}$ given by

\begin{displaymath}S_{ij}=(\hat\nabla \phi_{4p+i}, \hat\nabla \phi_{4p+j})_{\hat K}, \quad i,j=0,\ldots,(p-1)^2-1,
\end{displaymath}

and

\begin{displaymath}T_{ik}=(\hat\nabla \phi_{4p+i}, \hat\nabla \phi_k)_{\hat K}, \quad i=0,\ldots,(p-1)^2-1,\quad k=0,\ldots,4p-1.
\end{displaymath}

The solutions to problem (7) for l=1,2 are

\begin{displaymath}(p_{4p+j})_l=-\sum_{i=0}^{(p-1)^2-1}\sum_{k=0}^{4p-1}S^{-1}_{ji}T_{ik}(p_k)_l, \quad j=0,\ldots,(p-1)^2-1,
\end{displaymath}

that may be written in compact form:

 \begin{displaymath}
p_{4p+j}=\sum_{k=0}^{4p-1}c_{jk}p_k, \quad j=0,\ldots,(p-1)^2-1,
\end{displaymath} (8)

where cjk represents the coefficient

 \begin{displaymath}
c_{jk}=-\sum_{i=0}^{(p-1)^2-1}S^{-1}_{ji}T_{ik}.
\end{displaymath} (9)

of the linear combination (8), that represents the dependency of the jth inner mapping support point p4p+j on the support points pk, $k=0,\ldots,4p-1$, that are placed on the boundary of the cell K. For a fixed degree p, these coefficients cjk are the same for the mapping of all cells K in real space because the cjk depend only on the reference element $\hat K$. Therefore the coefficients cjk can be precomputed and result in following linear combinations:

For p=2 the linear combination turns out to be

\begin{displaymath}p_8=\tfrac{1}{16}\sum_{k=0}^3p_k+\tfrac{3}{16}\sum_{k=4}^7p_k,
\end{displaymath}

see also Figure 3, left.
  
Figure 3: Left: Coefficients c8,k for Q2 mapping. Right: Coefficients c12,k for Q3 mapping.
\includegraphics*[scale=0.4]{figures/m2.ps}

For the case that p=3, Figure 3, right, shows the coefficients c12,k of the linear combination for the inner mapping support point p12. The coefficents for the points p13, p14 and p15 can be obtain by rotation of the coefficients.

Implementation in deal.II. The coefficients cjk, see (9), are represented in the MappingQ class by the laplace_on_quad_vector as follows

\begin{displaymath}\mbox{\textit{laplace\_on\_quad\_vector}[j][k]}=c_{jk}
\end{displaymath}

These coefficients are the same for the mapping of all cells K in real space because the cjk depend only on the reference element $\hat K$. Hence for a given degree p the laplace_on_quad_vector can be filled by the constructor of the MappingQ class. This is done by calling the MappingQ::set_laplace_on_quad_vector function that includes the coefficients hardcoded for p=2 and p=3 in d=2 dimensions, and a routine for computing the coefficients according to (9) for all other cases. The mapping support points $p_i\in K$, $i=0,\ldots,(p+1)^d-1$ are computed once for each MappingQ object by the MappingQ::compute_support_points_laplace that is invoked by the virtual compute_mapping_support_points function of the base Mapping class. In MappingQ::compute_support_points_laplace, first the 4p points on the boundary of the cell are computed (by calling MappingQ::add_line_support_points), then by calling MappingQ::apply_laplace_vector the remaining (p-1)2 inner mapping supports points are computed, where MappingQ::apply_laplace_vector just performs the linear combination given in (8).

Ralf Hartmann
2001-09-03