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.