I have implemented a pretty straightforward finite element solver for the following Poisson equation. For the purposes of this question we can assume the source term and the Dirichlet data both vanish. I assume that $g\in H^{-1/2}(\partial \Omega_N) \cap L^2(\partial \Omega_N)$ (instead of just $g\in H^{-1/2}(\partial\Omega_N)$)
\begin{aligned} \mathbf{q} + \mathbf{\nabla}u &= 0 & \mathrm{in\,} \Omega,\\ \mathbf{\nabla} \cdot \mathbf{q} &= 0& \mathrm{in\ } \Omega,\\ u& = 0 &\mathrm{on\ } \partial\Omega_D,\\ \mathbf{q}\cdot \mathbf{\eta} &= g & \mathrm{on\ } \partial\Omega_N. \end{aligned}
I accounted for Neumann boundary conditions on all or some of the boundary by using Lagrange multipliers as follows.
Let the The domain $\Omega$ have a triangulation $T$ made up of simplexes $K$. We denote $$\mathbf{W}_h= \{\mathbf{p}\in \mathbf{H}^{\mathrm{div}}(\Omega)\colon \forall K\in T \quad \mathbf{p}\big|_K \in \mathrm{RT}^k(K))$$ $$V_h = \{v\in L^2(\Omega)\colon \forall K \quad v\big|_K \in P^k(K)\}$$ $$X_h = \{\mu \in L^2(\partial \Omega_N)\colon \forall K\in T\,\, \forall \mathrm{edge}\subset\partial K\cap \partial\Omega_N\quad \mu \big|_{\mathrm{edge}} \in P^k(\mathrm{edge})\}$$
Where $\mathrm{RT}^k(K)$ is the standard Raviart Thomas space and $P^k(K)$ is the space of polynomials of total degree no more than $k$.
The problem is to find the triple $(\mathbf{q}_h,u_h,\lambda_h) \in \mathbf{W}_h \times V_h \times X_h$ so that for all test functions $(\mathbf{p}, v, s) \in (W_h \times V_h \times X_h)$ we have \begin{aligned} 0 &= \int_\Omega \mathbf{p}\cdot \mathbf{q}_h\, \mathrm{d} x - \int_\Omega (\mathbf{\nabla} \cdot \mathbf{p}) u_h\, \mathrm{d} x + \int_{\partial\Omega_N} \mathbf{p} \lambda_h \cdot \mathbf{\eta}\, \mathrm{d} x \\ 0 &= \int_\Omega v (\mathbf{\nabla} \cdot \mathbf{q}_h)\, \mathrm{d} x\\ \int_{\partial\Omega_N}s g \, \mathrm{d}x&= \int_{\partial\Omega_N} s\mathbf{q}_h \cdot \mathbf{\eta}\, \mathrm{d} x \end{aligned}
The scheme works (well defined, converges to true solution, etc.) and is (I think) close to standard practice. However, in all the references I have seen in the literature (On the Mixed Finite Element Method with Lagrange Multipliers Babuska and Gatica, 2003 Numerical Methods for Partial Differential Equation and section 4.4 of A simple Introduction to the Mixed Finite Element Method by Gatica) require that the trial and test space for the Lagrange multiplier, $X_h$, have elements that are continuous across elements of the "boundary mesh".
- Is this scheme close to standard practice?
- Are there any references I could read that discuss this scheme?
- What problems am I allowing by choosing a piecewise discontinuous space for the Lagrange multiplier?
- A disadvantage of the scheme is that it requires extra regularity on $g$. Are there more disadvantages?
