\(\DeclareMathOperator{\opdiv}{div}\)

51. Hybridization Techniques

Discretizing an elliptic equation by a primal method leads to a positive definite matrix. Using a mixed method leads to a saddle point problem. This disadvantage can be overcome by the so called hybridization technique:

The idea is as follows: Break the normal-continuity of the RT - space, and reinforce it by another Lagrange parameter. This new Lagrange parameter is a polynomial on every mesh edge (or face, in 3D).

The variational formulation is as follows:

Find: \(\sigma_h \in RT_k^{dc}, u_h \in P^k, \widehat u_h \in P^k(\cup E)\) such that

\[\begin{split} \begin{array}{ccccccl} \int_\Omega \sigma_h \tau_h & + & \sum_T \int_T \opdiv \tau_h \, u_h & - & \sum_E \int_E [\![\tau_{h,n}]\!] \widehat u_h & = & 0 & \forall \, \tau_h \\ \sum_T \int_T \opdiv \sigma_h \, v_h & & & & & = & \int f v_h & \forall \, v_h \\ -\sum_E \int_E [\![\sigma_{h,n}]\!] \widehat v_h & & & & & = & 0 & \forall \, \widehat v_h \end{array} \end{split}\]

We used the definition of the jump term \([\![ \tau_n ]\!] = \tau_1 n_1 + \tau_2 n_2\)

This formulation gives the same solution as the mixed formulation, so we don’t need an extra error analysis

The physical meaning of the Lagrange paramter \(\widehat u_h\) is the primal variable, what can be seen by integration by parts of the first equation.

Now, Dirichlet boundary conditions are set by constraining the \(\widehat u\) on the Dirichlet boundary, and Neumann boundary conditions are formulated by \(\int_{\Gamma_N} g \widehat v\). The hybridized formulation is thus similar to a primal method.

The discretization system has the form

\[\begin{split} \left( \begin{array}{ccc} A & B_1^T & B_2^T \\ B_1 & & \\ B_2 & & \end{array} \right) \left( \begin{array}{c} \sigma \\ u_1 \\ u_2 \end{array} \right) = \left( \begin{array}{c} 0 \\ f \\ 0 \end{array} \right) \end{split}\]

The submatrix

\[\begin{split} \left( \begin{array}{cc} A & B_1^T \\ B_1 & 0 \end{array} \right) \end{split}\]

is regular and block-diagonal, each block corresponds to an element. Thus it can be cheaply inverted.

If \(u_2\) were known, we could compute the first two variables from \(u_2\):

\[\begin{split} \left( \begin{array}{c} \sigma \\ u_1 \end{array} \right) = \left( \begin{array}{cc} A & B_1^T \\ B_1 \end{array} \right)^{-1} \left[ \left( \begin{array}{c} 0 \\ f \end{array} \right) - \left( \begin{array}{c} B_2^T \\ 0 \end{array} \right) u_2 \right] \end{split}\]

Plugging this term into the third equation \(B_2 \sigma = 0\), we obtain the system

\[\begin{split} \left( B_2 \; \; 0 \right) \left( \begin{array}{cc} A & B_1^T \\ B_1 \end{array} \right)^{-1} \left( \begin{array}{c} B_2^T \\ 0 \end{array} \right) \; \; u_2 = \left( B_2 \; \; 0 \right) \left( \begin{array}{cc} A & B_1^T \\ B_1 \end{array} \right)^{-1} \left( \begin{array}{c} 0 \\ f \end{array} \right) \end{split}\]

The matrix on the left hand side is symmetric positive definite. It behaves like a standard system matrix (e.g. for condition number), and standard iterative methods and preconditioners can be used for solution.

The lowest order hybrid method has the same degrees of freedom as the non-conforming \(P^1\) element. Here, the extension to higher order is straight forward.

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

order=2
Sigma = Discontinuous (HDiv(mesh, order=order, RT=True))
V = L2(mesh, order=order)
Vhat = FacetFESpace(mesh, order=order, dirichlet="left|bottom")

X = Sigma*V*Vhat

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

n = specialcf.normal(mesh.dim)

a = BilinearForm(X, eliminate_internal = True)
a += (sigma*tau + div(sigma)*v + div(tau)*u)*dx
a += (-sigma*n*vhat-tau*n*uhat) * dx(element_boundary=True)

c = Preconditioner(a, "bddc")
a.Assemble()

f = LinearForm(X)
f += -v*dx
f.Assemble()

gfu = GridFunction(X)

f.vec.data += a.harmonic_extension_trans * f.vec

# gfu.vec.data = a.mat.Inverse(X.FreeDofs(True)) * f.vec
from ngsolve.solvers import CG
CG (mat=a.mat, pre=c.mat, rhs=f.vec, sol=gfu.vec, 
    printrates=True, maxsteps=200)

gfu.vec.data += a.harmonic_extension * gfu.vec
gfu.vec.data += a.inner_solve * f.vec

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

CG iteration 1, residual = 0.37594503814550667     
CG iteration 2, residual = 0.07563602425424365     
CG iteration 3, residual = 0.013738193462247741     
CG iteration 4, residual = 0.004300349628057608     
CG iteration 5, residual = 0.0015837080803695188     
CG iteration 6, residual = 0.0006602691617452845     
CG iteration 7, residual = 0.00018078405517481033     
CG iteration 8, residual = 6.109953797594372e-05     
CG iteration 9, residual = 2.234404212444103e-05     
CG iteration 10, residual = 6.648411830187808e-06     
CG iteration 11, residual = 2.302264315279244e-06     
CG iteration 12, residual = 6.302714596286466e-07     
CG iteration 13, residual = 2.179130114255554e-07     
CG iteration 14, residual = 6.763102225221426e-08     
CG iteration 15, residual = 2.1915768422357135e-08     
CG iteration 16, residual = 7.95474670424036e-09     
CG iteration 17, residual = 2.2189016160123873e-09     
CG iteration 18, residual = 8.744737564650945e-10     
CG iteration 19, residual = 2.4650109967698113e-10     
CG iteration 20, residual = 8.6487932956182e-11     
CG iteration 21, residual = 2.7909557666719956e-11     
CG iteration 22, residual = 7.290001584068349e-12     
CG iteration 23, residual = 2.922852287859793e-12     
CG iteration 24, residual = 6.323104671078037e-13     
CG iteration 25, residual = 2.3168473544001734e-13