Raviart-Thomas elements on reference square

Question

I'd like to learn how the Raviart-Thomas (RT) element works. To that end I'd like to analytically describe how the basis functions look on the reference square. The goal here is not to implement it myself, but rather just to get an intuitive understanding of the element.

I'm largely basing this work off of the triangular elements discussed here, perhaps extending it to quadrilaterals is a mistake in itself.

That said, I can define the basis functions for first RK element RK0:

$$\mathbf{\phi}_i(\mathbf{x}) = \mathbf{a} + \mathbf{b}\mathbf{x} = \begin{pmatrix} a_1 + b_1 x\\a_2 + b_2 y\end{pmatrix}$$ for $i = 1,\dots,4.$

The conditions on $\mathbf{\phi}_i$ are that:

$$\mathbf{\phi}_i(\mathbf{x}_j)\cdot\mathbf{n}_j = \delta_{ij}$$

where $\mathbf{n}_j$ is the unit normal shown below, and $\mathbf{x}_j$ is its coordinate.

This is the reference square $[-1,1]\times[1,1]$, so this leads to a system of equations for each basis function. For $\mathbf{\phi}_1$ this is:

$$\begin{pmatrix} 1 & 0 & 1 & 0\\ 0 & -1 & 0 & 1\\ -1 & 0 & 1 & 0\\ 0 & 1 & 0 & 1\end{pmatrix}\begin{pmatrix}a_1 \\ a_2\\ b_1\\ b_3\end{pmatrix} = \begin{pmatrix} 1\\0\\0\\0\end{pmatrix}$$

which can be solved to give:

$$\mathbf{\phi}_1(\mathbf{x}) = \frac{1}{2}\begin{pmatrix} 1 + x\\ 0\end{pmatrix}$$

The other basis functions can be found similarly.

Assuming this is correct, the next step is to find the basis functions for RK1. This is where I'm getting a little unsure of myself. According to the link above, the space we are interested in is:

$$P_1(K) + \mathbf{x}P_1(K)$$

A basis for $P_1$ would be $\{ 1, x, y\}$

I think this means the RK1 basis functions should take the form:

$$\mathbf{\phi}_i(\mathbf{x}) = \begin{pmatrix} a_1 + b_1 x + c_1 y + d_1 x^2 + e_1 xy\\ a_2 + b_2 x + c_2 y + d_2 xy + e_2 y^2\end{pmatrix}$$

This leaves 10 unknowns for each basis function. If we apply the same conditions as in the RK0 case, namely:

$$\mathbf{\phi}_i(\mathbf{x}_j)\cdot\mathbf{n}_j = \delta_{ij}$$, where $\mathbf{n}_j$ is the unit normal as shown below:

this gives us 8 equations. The other 2 I think can be found from some moments. I'm not really sure how exactly. The link above talks about integrating against a basis for $[P_1]^2$, but I'm having trouble figuring out what that means. Am I on the right track, or have I completely missed something here?

Answer 1

In general, you cannot just transfer the same polynomial basis from tetrahedral to quadrilateral elements.¹ In particular, the whole point of quadrilateral elements is to work with tensor products of one-dimensional polynomials, which is not possible for tetrahedral elements.

There are in fact quadrilateral Raviart-Thomas elements, but their definition is different. In two dimensions, the polynomial space for $RT_k$ is given by $$P_{k+1,k} \times P_{k,k+1},$$ where $$P_{k,l} = \left\{\sum_{i=0}^k\sum_{j=0}^l a_{ij} x^iy^j: a_{ij}\in\mathbb{R}\right\}.$$ A typical polynomial for $k=0$ would thus be as you wrote, but for $k=1$ it would be $$\begin{pmatrix} a_1 + b_1 x + c_1 x^2 + d_1 y + e_1 xy + f_2 x^2y\\ a_2 + b_2 y + c_2 y^2 + d_2 x + e_2 xy + f_2 xy^2 \end{pmatrix}.$$ Hence, $\dim RT_1 = 12$, and in general, $\dim RT_k = 2(k+1)(k+2)$. This means you need two additional degrees of freedom, which should be located on the interior of the element. (In general, for $RT_k$ you take $k+1$ normal derivatives on each facet, and the remaining degrees of freedom from the interior.)

To answer your actual question: For Raviart-Thomas elements, you usually take moments rather than point evaluations, i.e., the remaining conditions come from the conditions $$ \int_{-1}^1 \int_{-1}^1 \phi_i(x,y) q_j(x,y) dx\,dx=\delta_{ij},$$ where the $\{q_j\}$ are a basis of $P_{k-1,k}\times P_{k,k-1}$ (e.g., $\{1,x,y\}$ for $k=1$). To make it easier to get a full nodal basis, the facet degrees of freedom are usually not taken as point evaluation, but also as moment conditions: $$ \int_{e_m} \phi_i(s)^T\nu_{e_m} \, q_{m,j}(s)\,ds,$$ where $e_m$ is one of the four edges, $\nu_{e_m}$ is the corresponding outer normal, and for each $m$, the $q_{m,j}$ form a basis of $P_k(e_m)$ (e.g., $\{1,x\}$ or $\{1,y\}$ for $k=1$ depending on the edge orientation). Together, these degrees of freedom are unisolvent (i.e., the corresponding system of basis functions is always invertible).

You can find a discussion on quadrilateral Raviart-Thomas-Elements in Chapter 2.4.1 of Boffi, Brezzi, Fortin: Mixed Finite Element Methods and Applications, Springer 2013, Arnold, Boffi, Falk: Quadrilateral $H(div)$ finite elements, SINUM 42(5), 2005, pp. 2429-2451, and Chapter 3.2.3 of Ronald Hoppe's lecture notes.

---

_{1. As a rule of thumb, a polynomial space of order $k$ on tetrahedral elements contains monomials whose powers sum to $k$, while a space of order $k$ on quadrilateral elements contains monomials whose maximal power is $k$. For example, $x^2y$ would be of order $3$ on tetrahedrals, but only of order $2$ on quadrilaterals.}