An example of mixed elliptic problem using lowest-order Raviart Thomas element

Question

I try to solve the following mixed second order elliptic PDE in the domain $D=[0, 1]^2$

\begin{eqnarray*} v+\nabla p=&0 \quad &\text{in} \quad D,\\ \text{div}(v)=&1/2 \quad &\text{in} \quad D,\\ v\cdot n =&0 \quad &\text{on} \quad \partial D,\\ \end{eqnarray*} where $n$ is the unit outward normal vector for $D$. Divided the domain into four sub-squares: $[0,0.5]^2$, $[0.5,1] \times [0, 0.5]$, $[0,0.5] \times [0.5,1]$ and $[0.5,1]^2$, denote the partition as $\mathcal{T}^h$. We want to solve the above problem using (rectangular) lowest-order Raviart Thomas element $RT_0$.

If we take the reference square to be $[0,1]^2$, then we can compute four $RT$ basis functions as follow $$ \psi_1(x_1,x_2) = (x_1,0)^T $$ $$ \psi_2(x_1,x_2) = (0, x_2-1)^T $$ $$ \psi_3(x_1,x_2) = (x_1-1, 0)^T $$ $$ \psi_4(x_1,x_2) = (0,x_2)^T.$$

One can transform these basis functions to our sub-squares defined above. For example, in $[0,0.5]^2$, $$\tilde{\psi_1} (x_1,x_2) = (2x_1,0)^T.$$

The weak formulation of the above PDE is: find $(v,p)\in V_h \times Q$ such that \begin{eqnarray*} \int_D v \cdot w - \int_D p \text{div}(w)= & 0 \quad & \forall w \in V_h,\\ \int_D \text{div}(v) q = &\frac{1}{2}\int_D q \quad & \forall q \in Q, \end{eqnarray*} where $V_h = \{ w\in \big(L^2(D)\big)^2: w\cdot n =0\}$ and $Q = \{ q: q|_K = const, K\in \mathcal{T}^h\}.$ After we derive the weak formulation, the finite element system will look like the following $$ \left [ \begin{array}{cc} B & C \\ C^T & 0 \end{array} \right ] \left [ \begin{array}{c} \tilde{v} \\ \tilde{p} \end{array} \right ] = \left [ \begin{array}{c} f_v \\ f_p \end{array} \right ]. $$ Here is my question: In our special case, it should be $B,C \in \mathbb{R}^{4\times 4}$ and $\tilde{f_v}$ zero vector. I am not clear that which $RT$ basis functions I should consider when forming the finite element matrix. I think there is one basis function per interior edge for $RK_0$. Also, I want to know how the matrices $B$ and $C$ and the vector $\tilde{f_p}$ exactly look like. (Both the numeric and analytical expression)

Another question is: What is the compact support of these $RT$ basis functions from the corresponding sub-squares? Say $\tilde{\psi_1}$ in the sub-square $[0,0.5]^2$, what is the compact support of $\tilde{\psi_1}$? How about the other $\tilde{\psi_i}$ in the same sub-squares $[0,0.5]^2$?

Remark: To approximate $v$, we use $RT_0$ basis functions. For $p$, we use the piecewise constant function.

Answer 1

The size of the system depends on whether you eliminate the degrees of freedom from your system or not. If you eliminate them, then $B$ will be a 4x4 matrix (one unknown per interior edge), $C$ will be also 4x4 (but in general case it will be rectangular).

Assuming that you enumerate first dofs for $x$-components and then for $y$-components in a standard way, $B$ will be a block diagonal matrix with blocks corresponding to $x$ and $y$-directions. Each block is (before elimination of boundary dofs) a tridiagonal matrix with a row $$ \left\{ h_x h_y \frac{1}{6}, h_x h_y \frac{4}{6}, h_x h_y \frac{1}{6} \right\} $$

If you again use the standard enumeration for dofs for $p$, matrix $C$ before elimination will be a block matrix which basically gives you simply the two-point approximation of the $- \nabla$ operator, with entries $\pm h_x$ and $\pm h_y$.

The support of each basis function for $\mathbf{v}$ are the two rectangles which shared the edge with the corresponding degree of freedom.

Vector $\tilde{f_p}$ has the elements $\frac{1}{2}\int_D \chi_j$ where $\chi_j$ is the characteristic function of a mesh cell, that is, it will be just the area which equals $h_x h_y$.

I suggest you first write down the local $B_{loc}$ and $C_{loc}$ matrices and then look at the assembling procedure which will produce the global matrix in a straight-forward way.