Facet spaces and hybrid methods

Mixed methods for second order problems lead to saddle point problems, and indefinite matrices. By hybridization one obtains a positive definite system again. It's structure is similar to the non-conforming $P^1$ method, but hybridization works for any order. See text-book by Brezzi and Fortin.

One skips the normal-continuity of the $H(div)$ variable, and reinforces it by a Lagrange parameter. This leads to the following discrete system:

Find $\sigma, u, \widehat u \in \Sigma_h \times V_h \times F_h$:

$$\DeclareMathOperator{\Div}{div} \begin{array}{ccccccll} \int \lambda^{-1} \sigma \tau & + & \sum_T \int_T \Div \tau \, u & + & \sum_F \int_F [\tau_n] \widehat u & = & 0 & \forall \, \tau \in \Sigma \\ \int \Div \sigma \, v &&&&& = & \int f v & \forall \, v \in V_h \\ \int [ \sigma_n ] \, \widehat v &&&&& = & \int_{\Gamma_n} g \widehat v & \forall \, \widehat v \in F_h \end{array}$$

where $\Sigma_h$ is an discontinuous $H(div)$ finite element space, $V_h$ a sub-space of $L_2$, and $F_h$ consists of polynomials on every edge.

from netgen.geom2d import unit_square
from ngsolve import *
import netgen.gui
%gui tk
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

same example as in 'mixed':

source = sin(3.14*x)
ud = CoefficientFunction(5)
g = CoefficientFunction([y*(1-y) if bc=="left" else 0 for bc in mesh.GetBoundaries()])
lam = 10

define spaces:

• The discontinuous flag generates an element-wise $H(Div)$-space
• FacetFESpace lives only on facets (i.e. faces in 3D, edges in 2D, points in 1D)

Boundary conditions are now posed for the facet-space

order = 3
V = HDiv(mesh, order=order, flags = { "discontinuous" : True })
Q = L2(mesh, order=order-1)
F = FacetFESpace(mesh, order=order, dirichlet="bottom")
X = FESpace([V,Q,F])

Assemble forms. The jump-term is rewritten as $$\sum_F \int_F [\sigma_n] v = \sum_T \int_{\partial T} \sigma_n v$$

sigma,u,uhat = X.TrialFunction()
tau,v,vhat = X.TestFunction()

a = BilinearForm(X)
a += SymbolicBFI(1/lam * sigma*tau + div(sigma)*v + div(tau)*u)
n = specialcf.normal(mesh.dim)
a += SymbolicBFI(sigma*n*vhat+tau*n*uhat, element_boundary=True)
a.Assemble()

f = LinearForm(X)
f += SymbolicLFI(-source*v)
f += SymbolicLFI(g*vhat, BND)
f.Assemble()

Solve system. Leave everything to the sparse direct solver.

gfu = GridFunction(X)
inv = a.mat.Inverse(freedofs=X.FreeDofs())
gfu.vec.data = inv * f.vec

Draw (gfu.components[0], mesh, "sigma")
Draw (gfu.components[1], mesh, "u")