HDG mixed Poisson on hexa/Tet or Quad/Tri do not have the same results

Hello everyone

That’s the first time I post here, thanks a lot for all the developpement of this code, that’s incredible.
Nowadays, I am intersted in using the HDG versions of the mixed Poisson equation (to compute temperature flux or Darcy flow in porous media).

I have successfully implemented that, so I am happy :slight_smile:

Therefore, my next steps should be to add holes in the mesh but during my investigations, I discover that the results of the flux or the scalar field depends on the type of elements of the mesh (hexa, prism or tet in 3D) and (tri/quad in 2D). The results are perfect for tri/tet elements. Do you have any hint about that?

I did a small modification of the tutorial to build an example with an analytic function.

Elements for a comparison: In this case, the shape of the functions are perfect for all cases except that the magnitude of u is not good.

from ngsolve import *
from ngsolve.webgui import Draw

from meshes import MakeStructured3DMesh
from meshes import MakeStructured2DMesh

test3D = False
testTuto = False

if test3D:
    mesh = MakeStructured3DMesh(nx=11,ny=11,nz=11,prism=False,hexes = False)
    source = -12 * pi**2 * sin(2 *pi*x)* sin(2 *pi*y)* sin(2 *pi*z)
    mesh = MakeStructured2DMesh(nx=51,ny=51,quads = True)
    source = -8 * pi**2 * sin(2 *pi*x)* sin(2 *pi*y)

if testTuto:
    source = sin(x * pi)
    ud = CF(5)
    ud = CF(0)
lam = 1

order = 1
V = Discontinuous (HDiv(mesh, order=order))
Q = L2(mesh, order=order-1)
F = FacetFESpace(mesh, order=order, dirichlet="bottom|right|top|left|front|back")

X = V*Q*F
print ("sigmadofs:", X.Range(0))
print ("udofs:    ", X.Range(1))
print ("uhatdofs: ", X.Range(2))

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  

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

gfu = GridFunction(X)
gfu.components[2].Set(ud, BND)

if a.condense:

    fmod = (f.vec + a.harmonic_extension_trans * f.vec).Evaluate()

    import numpy as np
    inv = a.mat.Inverse(X.FreeDofs(True))
    gfu.vec.data +=  inv * (fmod- a.mat * gfu.vec)
    gfu.vec.data += a.harmonic_extension * gfu.vec
    gfu.vec.data += a.inner_solve * f.vec


    r = f.vec - a.mat * gfu.vec
    inv = a.mat.Inverse(freedofs=X.FreeDofs())
    gfu.vec.data += inv * r

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

vtk = VTKOutput(mesh,coefs=[gfu.components[0],gfu.components[1],source],names=["sigma","u","source"],filename=f"HDG")

Here the output to compare solution from quad and triangulation.

u for Quad

u for Tri

sigma for Quad

sigma for Tri


For quads you need Q of same order as V.


1 Like

Thanks a lot Christopher :slight_smile:

As expected it works perfectly.