This page was generated from jupyter-files/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 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’:

[2]:

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

[3]:

order = 3
V = HDiv(mesh, order=order, discontinuous=True)
Q = L2(mesh, order=order-1)
F = FacetFESpace(mesh, order=order, dirichlet="bottom")
X = FESpace([V,Q,F])
print ("sigmadofs:", X.Range(0))
print ("udofs:    ", X.Range(1))
print ("uhatdofs: ", X.Range(2))

sigmadofs: slice(0, 1120, 1)
udofs:     slice(1120, 1456, 1)
uhatdofs:  slice(1456, 1832, 1)

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()

condense=True

a = BilinearForm(X, eliminate_internal=condense)
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)

c = Preconditioner(a, "bddc")

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

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 condense:

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

    solvers.CG(mat=a.mat, pre=c.mat, rhs=f.vec, 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.CreateVector()
    r.data = f.vec - a.mat * gfu.vec
    inv = a.mat.Inverse(freedofs=X.FreeDofs())
    gfu.vec.data += inv * r

it =  0  err =  136.8703468856364
it =  1  err =  60.58702491012707
it =  2  err =  24.172480309852325
it =  3  err =  12.551902691915576
it =  4  err =  4.6561221507187
it =  5  err =  1.6707046939551382
it =  6  err =  0.6483905966924745
it =  7  err =  0.25370556589294635
it =  8  err =  0.10118909278378112
it =  9  err =  0.0362644711098403
it =  10  err =  0.021422651202624675
it =  11  err =  0.007698393630516215
it =  12  err =  0.003262626544392722
it =  13  err =  0.0012157696030631927
it =  14  err =  0.0005171305441789332
it =  15  err =  0.00018478710262344629
it =  16  err =  6.140422208014876e-05
it =  17  err =  2.3634731562752536e-05
it =  18  err =  8.96567905629904e-06
it =  19  err =  4.210915291606238e-06
it =  20  err =  1.09563838139071e-06
it =  21  err =  5.00341691606555e-07
it =  22  err =  1.917170519922857e-07
it =  23  err =  8.852599023217076e-08
it =  24  err =  3.7509955514489134e-08
it =  25  err =  1.2734631340019504e-08
it =  26  err =  3.9327806777442605e-09
it =  27  err =  1.5953738573307435e-09
it =  28  err =  7.259436128069845e-10
it =  29  err =  2.6189206224604935e-10

[6]:

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

[ ]: