Hello,
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:
Code:
uhat = IfPos(n, u, u.Other(bnd=1))
I am using the bilinear form element_boundary formulation
Code:
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...
Best,
Guosheng