Hybridization

2. Hybridization#

2.1. Poisson equation#

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

order = 2

Sigma = HDiv(mesh, order=order)
V = L2(mesh, order=order-1)
Vhat = FacetFESpace(mesh, order=order, dirichlet=".*")

X = PrivateSpace(Sigma) * V * Vhat
(sigma, u, uhat), (tau, v, vhat) = X.TnT()

print ("ndof.Sigma=", Sigma.ndof)
print ("Discont(Sigma).ndof=", Discontinuous(Sigma).ndof)
print ("Privat(Sigma).ndof=", PrivateSpace(Sigma).ndof)
print ("V.ndof=", V.ndof, "Vhat.ndof=", Vhat.ndof)
print ("X.ndof = ", X.ndof)

ndof.Sigma= 435
Discont(Sigma).ndof= 648
Privat(Sigma).ndof= 0
V.ndof= 162 Vhat.ndof= 273
X.ndof =  435

PrivateSpace .. have elements, but don’t assign global dofs
will be condensed out in assembling

Find \(\sigma \in BDM_k^{dc}\), \(u \in P^{k-1}\), \(\hat u \in P^k({\mathcal E})\):

\[\begin{split} \begin{array}{ccccll} (\sigma, \tau) & + & (\operatorname{div} \tau, u) + ([\sigma n], \hat u)_{\mathcal E} & = & 0 & \forall \, \tau \\ (\operatorname{div} \sigma, v) + ([\tau n], \hat v)_{\mathcal E} & & & = & (f,v) & \forall \, (v, \hat v) \end{array} \end{split}\]
n = specialcf.normal(mesh.dim)
def Div(sigma, v, vhat): return div(sigma)*v*dx - sigma*n * vhat*dx(element_boundary=True)
source = x

lhs = sigma*tau*dx + Div(sigma,v,vhat) + Div(tau, u,uhat)
rhs = source*v*dx
gf = Solve(lhs==rhs)
_, gfu, gfuhat = gf.components
print ("ndof gf =", len(gf.vec))
Draw (gfu);

ndof gf = 435

2.2. Reconstruct the flux#

Xp1 = Discontinuous(Sigma)*V
(sigmap1, up1), (taup1, vp1) = Xp1.TnT()

lhs = sigmap1*taup1*dx + div(sigmap1)*vp1*dx + div(taup1)*up1*dx
rhs = source*vp1*dx + taup1*n * gfuhat*dx(element_boundary=True)

gfpsigma,_ = Solve(lhs==rhs).components

Draw (gfpsigma[0], mesh);

Setup and solve a new problem for the flux, only if we need the flux.
Alternative:

  • Keep the flux in the product space, local condensation, and reconstruct algebraically by harmonic extension

  • Setup reconstruction operator

2.3. Postprocessing for the scalar#

obtain better scalar by solving

\[ \min_{u_p \atop \int_T u_p = \int_T u} \| \nabla u_p - \sigma \|^2 \]
Xp = L2(mesh,order=order+1) * L2(mesh,order=0)
(up,umean),(vp,vmean) = Xp.TnT()

lhs = grad(up)*grad(vp)*dx + up*vmean*dx + vp*umean*dx
rhs = gfpsigma*grad(vp)*dx + gfu*vmean*dx
gfp,_ = Solve(lhs==rhs).components
Draw (gfp);
Contents