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.

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

 Generate Mesh from spline geometry
 Boundary mesh done, np = 20
 CalcLocalH: 20 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 load internal triangle rules
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters

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, discontinuous=True)
# V = Discontinuous(HDiv(mesh, order=order))
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: [0,1080)
udofs:     [1080,1404)
uhatdofs:  [1404,1768)

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, condense=False)
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)

assemble VOL element 54/54
call wirebasket inverse ( with 86 free dofs out of 1768 )
A non-zerohas inverse
 76840

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

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

if a.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

assemble VOL element 54/54
assemble BND element 20/20
setvalues element 20/20
call pardiso ... done

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

BaseWebGuiScene