Question on interface integrals at material boundaries in DG formulation


I would like to perform DG calculation on a subdomain of a two-material mesh, but I have big trouble understanding the numerical flux, in particular the .Other() operator.

Below is a minimal test to illustrate my confusion.
I have a 1D domain with 4 elements, left half of them are labelled mat1, and right half are labelled mat2.
So, I want to solve u_t + u_x = 0 only on mat1 using DG-P0 formulation.
the numerical flux is taken to be upwinding:uhat = IfPos(n, u, u.Other(bnd=1))
I am using the bilinear form element_boundary formulation

a = -uhat*v*dx(element_boundary=True, definedon=mesh.Materials("mat1"))

But the code fails to produce a constant 1 solution on mat1.
Actually, when I apply the vector u=[1,1,0,0] to the bilinear form, I would assume to get a zero-vector, but instead, I get w=[0,1,0,0].
It is clear that the bilinear form forgot to integrate the boundary term at the interface shared by mat1 and mat2: for test function v=1 at the second cell, the bilinear form returns [u0] instead of [u0-u1] as my numerical flux suggest. Do you know what’s going on?

P.S. I also tried to use skeleton-based formulation, but I don’t know how to treat the interface properly…



Hi Guosheng,

element-boundary DG terms on region boundaries are skipped, see

comments in line 4024.

As a workaround, you can multiply with the characteristic function:

mat1func = CoefficientFunction( [1 if mat=="mat1" else 0 for mat in mesh.GetMaterials() ] )
a +=  -mat1func*uhat*v*n*dx(element_boundary=True) # , definedon=mesh.Materials("mat1")) 


Hi Joachim,

Thank you for the clarification. It would be better if we can include the material interface…

Here is another issue I encounter:
I want to do HDG on a subdomain, with static condensation on.
Attached is the code, and here is the set-up:
I have two subdomains, with label “inner” and “outer”.
Then, I have three finite element spaces, the vector DG space fes_q, and the scalar DG space fes_u, both live on the whole domain, and another hybrid FESpace live only on the skeleton of subdomain “inner”. (I don’t want to use HDG on “outer” subdomain)
Then, I solve the eqn
-lap u + u = 1 on domain “inner” with zero flux bc on the material interface.
The solution shall be u = 1 on iinner" and zero on “outer”.
To solve the problem, I have to manually set those DG dofs on subdomain “outer” as not FreeDOfs. (line 40-45)
Iin the attached code, when condensation is set off, the results are correst, but when condensation is set on, the solver produce nan values for those dofs in subdomain “outer”. I noticed that line 59 on a.inner_solve create the issue. Do you know what’s going on here?



Dear Guosheng,

You simply missed to restrict your element_boundary integrals to the interior domain.


my answer above applies only to DG, not to HDG.
HDG is much nicer since we can assemble element by element, and don’t need the neighbouring elements at the same time as in DG. The neighbouring elements outside the Region were the issue in your first post.


I see…
I thought the following two lines are mathematically equivalent:



(uhat*r*n)*dx(element_boundary=True, definedon=mesh.Materials("inner"))

But the first gives a buggy result, while the second produce the correct result. Can you explain the hidden difference?

Now, I have another (minor) question on the subject:
I want to add the bilinear form <uhat, vhat> on interface edges, following the post
I can do

a += SymbolicBFI(uhat.Trace()*vhat.Trace(), definedon=mesh.Boundaries("b2"))

This works great. But I would like to change SymbolicBFI to dx/ds language, so, I use

a +=uhat*vhat*dx(skeleton=True, definedon=mesh.Boundaries("b2"))

The code runs, but the above two lines produce different results. Again, I thought these two lines are identical… So, what’s the correct equivalent of the above SymbolicBFI in dx/ds language?

Thanks in advance. Your timely response really helped me a lot.


The first version produces 0 element-matrices, and tries to assemble them. But the global matrix has no entries allocated, and thus throws an exception.

The translate to dx/ds syntax should be

a +=uhat.Trace()*vhat.Trace()*ds(definedon=mesh.Boundaries("b2"))

We want a loop over boundary elements, no need for volume elements and thus no skeleton flag.

Thanks again for the clarification.
Previously, I used ds(skeleton=True,…), it gives me an error.
Now, I tried the same thing in 1D, and it seems .Trace() is not properly implemented yet for 1D Facet Dofs. It gives me this error message:

netgen.libngpy._meshing.NgException: FacetFESpace::GetSFE: unsupported element Point

Here is the code produce the error message