This page was generated from unit-2.7-hybrid/hybrid.ipynb.

2.7 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\):

\[\begin{split}\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}\end{split}\]

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.

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

same example as in 'mixed':

[2]:
source = sin(pi*x)
ud = CF(5)
g = mesh.BoundaryCF( {"left" : y*(1-y)}, default=0)
lam = 10

define spaces:

  • The Discontinuous FESpace-wrapper generates an element-wise discontinuous 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

[3]:
order = 3
V = Discontinuous (HDiv(mesh, order=order))
Q = L2(mesh, order=order-1)
F = FacetFESpace(mesh, order=order, dirichlet="bottom")
X = V*Q*F
print ("sigmadofs:", X.Range(0))
print ("udofs:    ", X.Range(1))
print ("uhatdofs: ", X.Range(2))
sigmadofs: [0,1120)
udofs:     [1120,1456)
uhatdofs:  [1456,1832)

Assemble forms. The jump-term is rewritten as

\[\sum_F \int_F [\sigma_n] v = \sum_T \int_{\partial T} \sigma_n v\]
[4]:
sigma,u,uhat = X.TrialFunction()
tau,v,vhat = X.TestFunction()

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

c = Preconditioner(a, "bddc")

f = LinearForm(X)
f += -source*v * dx - g*vhat.Trace() * ds

a.Assemble()
print ("A non-zero", a.mat.nze)

gfu = GridFunction(X)
A non-zero 6880

Solve system. Either we leave everything to the sparse direct solver, or use CG

[5]:
f.Assemble()
gfu.components[2].Set(ud, BND)

if a.condense:

    fmod = (f.vec + a.harmonic_extension_trans * f.vec).Evaluate()
    solvers.CG(mat=a.mat, pre=c.mat, rhs=fmod, sol=gfu.vec, initialize=False)

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

else:

    r = f.vec - a.mat * gfu.vec
    inv = a.mat.Inverse(freedofs=X.FreeDofs())
    gfu.vec.data += inv * r
CG iteration 1, residual = 136.14101710171593
CG iteration 2, residual = 59.97199896705699
CG iteration 3, residual = 24.09708229776066
CG iteration 4, residual = 12.365259367332891
CG iteration 5, residual = 4.500227573740078
CG iteration 6, residual = 1.5436222690215948
CG iteration 7, residual = 0.6253606118634962
CG iteration 8, residual = 0.23890072286213265
CG iteration 9, residual = 0.09410401770377781
CG iteration 10, residual = 0.03045052624547939
CG iteration 11, residual = 0.01423731192075639
CG iteration 12, residual = 0.004932886745754243
CG iteration 13, residual = 0.002146658066527776
CG iteration 14, residual = 0.0008801579888635546
CG iteration 15, residual = 0.0004053801141666873
CG iteration 16, residual = 0.0001073184650755895
CG iteration 17, residual = 4.241130888841856e-05
CG iteration 18, residual = 1.6149940767706976e-05
CG iteration 19, residual = 7.046696006269133e-06
CG iteration 20, residual = 2.5136423160207186e-06
CG iteration 21, residual = 6.758564102270819e-07
CG iteration 22, residual = 3.18214455629111e-07
CG iteration 23, residual = 1.1051209305651063e-07
CG iteration 24, residual = 5.400277686259754e-08
CG iteration 25, residual = 1.8901134514362357e-08
CG iteration 26, residual = 7.3118615775741006e-09
CG iteration 27, residual = 2.6927061522604897e-09
CG iteration 28, residual = 1.0284857201216218e-09
CG iteration 29, residual = 4.1485089070411144e-10
CG iteration 30, residual = 1.4204992823376103e-10
CG iteration 31, residual = 5.2792155687807855e-11
[6]:
Draw (gfu.components[0], mesh, "sigma")
Draw (gfu.components[1], mesh, "u");
[ ]:

[ ]: