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?


The Raviart-Thomas element has only 3 degrees of freedom on each triangle, namely the normal components of the vectors at the midpoints of the three edges. So the root of your misunderstanding is that $\mathbf u$ is a vector of size 3, matching the element matrix size.

  • $\begingroup$ ohhh. thanks a million! Prof. Wolfgang. But how it is that possible? The $\mathbf{u}$ at a point must have two or three components in reality. Maybe $\mathbf{u}$ is only a magnitude along the normal direction??? $\endgroup$ Nov 24 at 6:41
  • $\begingroup$ The degree of freedom is only the normal component of the vector, which is a scalar. $\endgroup$ Nov 24 at 9:01
  • $\begingroup$ understand better now. so in terms of hydraulics, it is the normal flux. Thanks! $\endgroup$ Nov 25 at 1:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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