We can use the $RT0$ to solve the Darcy equation, i.e. $$k^{-1}\mathbf{u}+\nabla p = 0, \text{ in } \Omega,$$ $$-\nabla \cdot \mathbf{u} = 0, \text{ in } \Omega,$$ $$p = p_D \text{ on } \partial\Omega,$$ $$\mathbf{u \cdot n} = 0 \text{ on } \partial\Omega, \tag{1}$$ where $\mathbf{u}$ is the velocity, $k$ is conductivity coefficient, and $p$ is the pressure.
With the incorporation of test functions, we have: $$\int_{T}k^{-1}\mathbf{u}\cdot\mathbf{v}\mathrm{d}T - \int_{T}p\nabla\cdot \mathbf{v}\mathrm{d}T=-\int_{\partial T}p_D(\mathbf{v\cdot n})\mathrm{d}(\partial T),$$ $$-\int_{T}\nabla\cdot \mathbf{u}q\mathrm{d}T=0,\tag{2}$$ where $\mathbf{v}$ and $q$ are basis functions (they can be replaced with specific basis functions, like $RT0$), and $T$ is an elementary triangle in a mesh. The Eq. (2) can be rewritten as $$\left[\begin{matrix}\mathbf{A}&\mathbf{B}\\\mathbf{B}^{TPS}&0\end{matrix}\right]\left[\begin{matrix}\mathbf{u}\\p\end{matrix}\right] = \left[\begin{matrix}\mathbf{f}\\0\end{matrix}\right], \tag{3}$$ which is a local (elementary) matrix.
The basis function $\mathbf{\hat{\phi}}$ of $RT0$ is written as $$\hat{\phi_1} = \sqrt{2}\left(\hat{x_1} \quad \hat{x_2}\right)^{T}, \hat{\phi_2} = \left(-1+\hat{x_1} \quad \hat{x_2}\right)^{T}, \text{ and } \hat{\phi_3} = \left(\hat{x_1} \quad -1+\hat{x_2}\right)^{T},\tag{4}$$ where $\hat{}$ is a reminder of reference triangle, $\hat{x_i}$ is coordinates of corners of the reference triangle, i.e. 0 or 1.
The transformation from $T$ to $\hat{T}$ is $$\begin{equation}\left[\begin{matrix}x_1\\x_2\end{matrix}\right]=B_k\left[\begin{matrix}\hat{x_1}\\ \hat{x_2}\end{matrix}\right]+b=\left[\begin{matrix}x_1^2-x_1^1&x_1^3-x_1^1\\x_2^2-x_2^1&x_2^3-x_2^1\end{matrix}\right]\left[\begin{matrix}\hat{x_1}\\ \hat{x_2}\end{matrix}\right] + \left[\begin{matrix}{x_1}\\ {x_2}\end{matrix}\right]\end{equation},\tag{5}$$ where $x_i$ is coordinates of nodes of $T$.
Now let's approximate $\mathbf{A}$ in Eq. (3) with $\mathbf{\hat{\phi}}$: $$\mathbf{A} = \alpha_{ij}^T = \frac{1}{|\text{det}B_k|}\int_{\hat{T}}k^{-1}[s_i][s_j](B_k\hat{\phi_i})\cdot(B_k\hat{\phi_j})\mathbf{d}\hat{T},$$ where $[s_i]$ is +1 or -1 if the normal direction of on edge $i$ of reference triangle is identical or not identical to the prescribed global normal direction of edge $i$. And we may use Gauss Quadrature Rule to approximate the integral, and the coordinates of Gauss points are $\hat{x_i}$ which can be input to $\hat{\phi_i}$ ($i = 1, 2, 3$ means the three edges).
My question is: since $B_k$ is $2 \times 2$ and $\hat{\phi_i}$ is $2\times 1$, the $B_k\hat{\phi_i}$ is $2\times 1$. Then the result of the dot product,i.e. $(B_k\hat{\phi_i})\cdot(B_k\hat{\phi_j})$, is a scalar. Thus, $\mathbf{A}$ is $3 \times 3$. However, $\mathbf{u}$ is $6 \times 1$, i.e. $[u_{x1}, u_{x2}, u_{x3}, u_{y1}, u_{y2}, u_{y3}]^{TPS}$ (if $\Omega$ is in 2D). They are not consistent in the number of rows. Why? Is my understanding incorrect?
