New Features

Let us consider the Poisson problem in $\mathrm{\Omega }\subset {\mathbb{R}}^d$, that is,

$$\begin{array}{rl}-\mathrm{\Delta }u& =f\text{ in }\mathrm{\Omega },\\ u& =0\text{ on }\partial \mathrm{\Omega }.\end{array}$$

from ngsolve import *
from netgen.geom2d import unit_square
from math import pi

mesh = Mesh(unit_square.GenerateMesh(maxh=1/4))

For such a problem, a meaningful choice of finite element spaces for an HDG method with $L^2$-lifting is the following:

$$\begin{array}{rl}{V}_h& =\{v_h\in L^2({\mathcal{T}}_h):v_h{|}_K\in P^k(K),\mathrm{\forall }K\in {\mathcal{T}}_h\}\\ {\hat{V}}_h& =\{{\hat{v}}_h\in L^2({\mathcal{F}}_h):{\hat{v}}_h{|}_F\in P^k(F),\mathrm{\forall }F\in {\mathcal{F}}_h;{\hat{v}}_h=0\text{on}\partial \mathrm{\Omega }\}\end{array}$$

The discrete weak formulation is then given by: With $[[v_h]]=v_h-{\hat{v}}_h$, find $(u_h,{\hat{u}}_h)\in V_h\times {\hat{V}}_h$ such that, for all $(v_h,{\hat{v}}_h)\in V_h\times {\hat{V}}_h$, there holds

$$\begin{array}{r}\sum _{K\in {\mathcal{T}}_h}({\int }_K\mathrm{\nabla }{u}_h\cdot \mathrm{\nabla }{v}_{h}dx-{\oint }_{\partial K}(\mathrm{\nabla }{u}_h\cdot {\mathbf{n}}_K)[[v_h]]ds-{\oint }_{\partial K}[[u_h]](\mathrm{\nabla }{v}_h\cdot {\mathbf{n}}_K)ds+{\oint }_{\partial K}\frac{k}{h}[[u_h]][[v_h]]ds+{\sigma }_{\ell }{\int }_K\mathcal{L}([[u_h]])\cdot \mathcal{L}([[v_h]])dx)=\sum _{K\in {\mathcal{T}}_h}{\int }_{K}f{v}_{h}dx.\end{array}$$

Here, the integral multiplied by ${\sigma }_{\ell }>1$ represents a lifting term which ensures stability for independent of $h$, $k$ and shape regularity. The lifting is defined element-by-element and depends only on element-local (including the element boundary) values. We define it using another FE space as follows:

$$\begin{array}{rl}& \text{For all }\phi \in L^2({\mathcal{F}}_h)\text{, the lifting }\mathcal{L}(\phi )\in {\mathbf{S}}_h\text{ is defined as the solution of }\\ & {\int }_K\mathcal{L}(\phi )\cdot {\mathbf{s}}_{h}dx={\oint }_{\partial K}\phi ({\mathbf{s}}_h\cdot {\mathbf{n}}_K)ds,\mathrm{\forall }{\mathbf{s}}_h\in {\mathbf{S}}_h,\mathrm{\forall }K\in {\mathcal{T}}_h,\\ & {\mathbf{S}}_h=\mathrm{\nabla }{V}_h=\{{\mathbf{s}}_h\in (L^2({\mathcal{T}}_h){)}^d:{\mathbf{s}}_h{|}_K\in (P^{k-1}(K){)}^d,\mathrm{\forall }K\in {\mathcal{T}}_h\}.\end{array}$$

However, the actual implementation of such a method looks slightly different. We introduce the compound FE space:

$$\begin{array}{rl}{X}_h& =\{(v_h,{\hat{v}}_h,{\mathbf{s}}_h):v_h\in V_h,{\hat{v}}_h\in {\hat{V}}_h,{\mathbf{s}}_h\in {\mathbf{S}}_h\}\end{array}$$

order = 8

V = L2(mesh, order=order)
Vhat = FacetFESpace(mesh, order=order, dirichlet=[1,2,3,4])
S = VectorL2(mesh, order=order-1, hide_all_dofs=True)
X = FESpace([V,Vhat,S])

In the VectorL2 FE space the "hide_all_dofs" flag marks all corresponding dofs as HIDDEN_DOFs (instead of LOCAL_DOF). Here are a few important properties of hidden d.o.f.s

• Hidden d.o.f.s will be removed (and will not be restored!) during static condensation. This has the following implications:
  • A hidden d.o.f. is a LOCAL_DOF in the sense that it only couples with other d.o.f.s on the same element
  • The entries of a r.h.s. vector of a linear system corresponding to a hidden d.o.f. are always zero
  • During static condensation the (element-local) matrix blocks for hidden d.o.f.s have to be invertible
  • No non-zero entries are reserved for hidden d.o.f.s in sparse matrices. Consequently, a hidden d.o.f. that is involved in a bilinear form will be ignored if not eliminated through static condensation.
• Hidden d.o.f.s (as other coupling types) only behave different from other coupling types when it comes to solving linear systems and static condensation is applied. This can be done with the "eliminate_internal" flag of a bilinear form or with the "eliminate_hidden" flag where the latter only applies static condensation of HIDDEN_DOFs and not for LOCAL_DOFs.

Then, the weak formulation of the implementation variant reads as follows. Find $(u_h,{\hat{u}}_h,{\mathbf{r}}_h)\in X_h$ such that, for all $(v_h,{\hat{v}}_h,{\mathbf{s}}_h)\in X_h$, there holds

$$\begin{array}{rl}\sum _{K\in {\mathcal{T}}_h}({\int }_K\mathrm{\nabla }{u}_h\cdot \mathrm{\nabla }{v}_{h}dx& -{\oint }_{\partial K}(\mathrm{\nabla }{u}_h\cdot {\mathbf{n}}_K)[[v_h]]ds-{\oint }_{\partial K}[[u_h]](\mathrm{\nabla }{v}_h\cdot {\mathbf{n}}_K)ds+{\oint }_{\partial K}\frac{k}{h}[[u_h]][[v_h]]ds)\\ & +{\sigma }_{\ell }\sum _{K\in {\mathcal{T}}_h}(-{\int }_K{\mathbf{r}}_h\cdot {\mathbf{s}}_{h}dx+{\oint }_{\partial K}[[u_h]]({\mathbf{s}}_h\cdot {\mathbf{n}}_K)ds+{\oint }_{\partial K}[[v_h]]({\mathbf{r}}_h\cdot {\mathbf{n}}_K)ds)=\sum _{K\in {\mathcal{T}}_h}{\int }_{K}f{v}_{h}dx.\end{array}$$

sigma = 2 
sigma_lift = 2 

(u,uhat,r),(v,vhat,s) = X.TnT()

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

jmp_u = u-uhat
jmp_v = v-vhat

a = BilinearForm( X, symmetric=True, eliminate_hidden=True)
a += SymbolicBFI( grad(u)*grad(v) )
a += SymbolicBFI( - (grad(u)*n)*jmp_v - jmp_u*(grad(v)*n) + (sigma/h)*jmp_u*jmp_v, element_boundary=True )
a += SymbolicBFI( - sigma_lift*(r*s) )
a += SymbolicBFI( sigma_lift*( jmp_u*(s*n) + jmp_v*(r*n) ), element_boundary=True )

f = LinearForm(X)
f += SymbolicLFI( 2*pi**2*sin(pi*x)*sin(pi*y)*v )

Here, a right-hand side term $f$ is chosen such that in the end, the exact solution is known and an error can be computed. Note that use "eliminate_hidden=True" here only (in contrast to a full "eliminate_internal=True"). This flag is not needed whenever "eliminate_internal=True" is used. However, whenever hidden DOFs are used, they have to be eliminated in the corresponding bilinear form, regardless of static condensation for higher order bubbles is used or not.

Now, one can assemble the linear system and take a look at the number of non-zero entries in the matrix A and its factors which arise due to static condensation:

a.Assemble()
f.Assemble()

print("a.mat has {} non-zero entries.".format(a.mat.nze))

a = BilinearForm( X, symmetric=True, eliminate_internal=True)
...
print("a.mat has {} non-zero entries.".format(a.mat.nze)) 
print("a.inner_solve has {} non-zero entries.".format(a.inner_solve.nze)) 
print("a.harmonic_extension has {} non-zero entries.".format(a.harmonic_extension.nze))

Comparison of non-zero entries

The numbers for the method using hidden DOFs and compression of the FE spaces with static condensation are:

a.mat has 22680 non-zero entries.
a.inner_solve has 105300 non-zero entries.
a.harmonic_extension has 63180 non-zero entries.

The output without hidden DOFs and with static condensation is as follows:

a.mat has 22680 non-zero entries.
a.inner_solve has 711828 non-zero entries.
a.harmonic_extension has 164268 non-zero entries.

Compress d.o.f.s

Hidden d.o.f.s don't play a role in the solution, but only a local role in defining some projection-type operators. Hence, we don't need to reserve a global d.o.f. number for a hidden d.o.f.. To prevent hidden d.o.f.s from obtaining a global d.o.f. number, we can use a wrapper which is called a CompressedFESpace that removes HIDDEN_DOFs and UNUSED_DOFs from the global numbering. For an arbitrary FESpace V you obtain a compressed spaces as follows:

print(V.ndof)
V = Compress(V)
print(V.ndof)

When applied to a CompoundFESpace we use 

print("dofs before compression:", X.ndof)
X = CompressCompound(X)
print("dofs after compression:", X.ndof)

resulting in a reduction of global d.o.f.s:

dofs before compression: 6876
dofs after compression: 3132

However, the non-zero entries in the matrices are essentially unaffected by the compression.

Projected jumps

Another example where hidden d.o.f.s can be used are "projected jumps" where the highest order functions in the FacetFESpace are made element-local to realize a facet-wise $L^2$ projection into the space of polynomials up to degree $k-1$. With the concept of hidden d.o.f.s this is easily achieved by changing only the one line in the definition of the FacetFESpace:

Vhat = FacetFESpace(mesh, order=order, dirichlet=[1,2,3,4], highest_order_dc=True, hide_highest_order_dc=True)

A fully functional python script is attached to this article. In this script, more options to play with are included.