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.

  • $\begingroup$ Don't we need another boundary condition? Can you write down the weak form? $\endgroup$
    – nicoguaro
    Apr 5, 2017 at 15:38
  • $\begingroup$ @nicoguaro The weak form is now added. I think it is enough to consider $v \cdot n =0$ on the boundary $\partial D$. $\endgroup$
    – Simon Pun
    Apr 5, 2017 at 16:26
  • 1
    $\begingroup$ The condition $v \cdot n=0$ is incompatible with your PDE since $\int_D div(v) = \int_D 1/2$ becomes $0 = \int_D 1/2$. $\endgroup$
    – cfdlab
    May 6, 2017 at 3:02

1 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.

  • $\begingroup$ What do you mean by "eliminate the degrees of freedom from the system"? $\endgroup$
    – Simon Pun
    Apr 6, 2017 at 2:01
  • $\begingroup$ There are two approaches, and what I like is to first assemble a global matrix using natural boundary conditions (I e. Including boundary dofs) and only then eliminate them by either reducing a matrix so that the vector of unknowns consist only of the actual dofs, or to replace the corresponding rows by something like $(0,...,0,1,0,...,0) $. And one should change the righthand side appropriately. $\endgroup$
    – VorKir
    Apr 6, 2017 at 2:06
  • $\begingroup$ Eliminate boundary dofs, mobile version makes me suffer, sorry $\endgroup$
    – VorKir
    Apr 6, 2017 at 2:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.