Neumann boundary conditions for Laplace equation with Raviart-Thomas elements

Question

I am working on creating a finite element model for the Darcy equation using Raviart-Thomas elements and the mixed hybrid formulation. The problem in mixed form is this:

$\mathbb{K}\nabla p = \vec{q}$

$\nabla\cdot\vec{q} = 0$

Where $\mathbb{K}$ is the diagonal matrix $\begin{bmatrix} k_0 & 0 \\ 0 & k_1\end{bmatrix}$

Now the weak formulation with Lagrange multipliers is: $\begin{align} \int_\Omega \mathbb{K}^{-1} \bar{\vec{q}}\cdot\vec{q} d \Omega + \int_\Omega p \nabla\bar{\vec{q}} d \Omega- \sum_E \iint_E \lambda\cdot(\bar{\vec{q}}\cdot\vec{n}_E) dE & = \int_\Gamma \bar{\vec{q}}\cdot\vec{n}p d \Gamma \\ \int_\Omega \bar{p}\nabla\cdot\vec{q} &= 0 \\ \sum_E \iint_E \bar{\lambda}\cdot(\vec{q}\cdot\vec{n}_E) dE &= 0 \end{align}$

Where a bar over the variable indicates it is a test function for that variable. Now I understand where the Dirchelet boundary conditions can be included, but for my case I would like to be able to also include Neumann boundary conditions. I found an article which does what I want to be included, but the mathematic are unfortunately too complicated for me (I'm an undergrad student Biomedical Engineering).

I am hoping that someone here can help me with this problem.

Answer 1

Are the sets $E$ edges of the triangulation? if so I don't see why you are applying your Lagrange multipliers on the edges (unless this is just a method for determining a constant that I haven't seen before), you want them on the whole domain so that you can impose the integral of a variable to be zero. Also $q$ is fully determined by the system, it is $p$ that may differ by a constant, so you want Lagrange multipliers for that one, taking a bar overhead to be a test function you want to add $$\int_\Omega\lambda\cdot\overline p+\overline\lambda\cdot p=0$$ to the system. Imagine you had $q=0$ on $\partial\Omega$, then your sum over edges would just be zero. If you want to impose Neumann boundary conditions for $q$, then we impose homogeneous BC's for $q$, and then $q$ may differ by a constant so we need a second lagrange multiplier $$\int_\Omega\lambda'\overline q+\overline\lambda' q=0$$ to the system. If you want nonhomogeneous Neumann BC's you do the same thing, and then enforce the nonhomogeneity as a boundary integral on the right hand side.

Answer 2

Neumann boundary conditions for the mixed formulation are imposed as "essential" boundaries. What that means is, you enforce it by directly setting $\vec{q}$ on the boundary. Remember that $RT_0$ gives you a single velocity degree of freedom on every edge. Since you have a Neumann condition, you want to directly restrict the velocity unknowns on the Neumann edges to the value of the boundary condition.

One way of accomplishing this is by first building the entire linear system of equation as if you have Dirichlet (pressure) boundaries. Locate in the linear system the rows that correspond to the Neumann boundaries. You can then force those unknowns to the value you need it to be by setting the row to zero and placing 1 on the diagonal. On the corresponding right hand side, place the actual flux value you want to set it to.

For example, let's say your linear system looks like this:

$$ \left(\begin{array}{cccc} a_{00} &a_{01} & a_{02} &a_{03}\\ a_{10} &a_{11} & a_{12} &a_{13}\\ a_{20} &a_{21} & a_{22} &a_{23}\\ a_{30} &a_{31} & a_{32} &a_{33} \end{array}\right) \left(\begin{array}{c} q_0\\ q_1\\ q_2\\ p_0 \end{array}\right)= \left(\begin{array}{c} 0\\ 0\\ 0\\ 0 \end{array}\right) $$

Let's assume that $q_0, q_1, q_2$ are the velocity unknowns corresponding to edges of a triangle. Say you want to set the boundary on $q_0 = 10$. The linear system would then become:

$$ \left(\begin{array}{cccc} 1 & 0 & 0 & 0\\ a_{10} &a_{11} & a_{12} &a_{13}\\ a_{20} &a_{21} & a_{22} &a_{23}\\ a_{30} &a_{31} & a_{32} &a_{33} \end{array}\right) \left(\begin{array}{c} q_0\\ q_1\\ q_2\\ p_0 \end{array}\right)= \left(\begin{array}{c} 10\\ 0\\ 0\\ 0 \end{array}\right) $$