Mixed FEM for diffusion problem with discontinuity and unstructured mesh

1

Following the example for MFEM and Poisson’s equation here, I am trying to solve the following system

58%20AM
where K is a spatially varying constant. I’m testing my code on a simple problem: square domain with the boundary conditions given as

09%20AM
and a discontinuous K field

07%20PM .
When I run this problem with a structured mesh, the solution for q is smooth and well behaved at the K discontinuity, so long as a string of edges matches the location of the discontinuity. Here is a plot of q overlaid with streamlines zoomed in on the K interface -

56%20AM
However, when I run the same problem with an unstructured mesh, the solution for q shows unphysical, sporadic oscillations. Of course, I have ensured that a string of edges in the unstructured mesh runs along the discontinuity in K, just as I did for the structure mesh. Here is the same plot, but for the unstructured mesh -

20%20AM
The solutions for h are essentially identical in both cases. They are smooth and well behaved.

My question may be more about FEM in general rather than FEniCS specific: Does anyone know why MFEM would have difficulties solving this system with an unstructured mesh (even with a string of mesh edges running along the discontinuity)? Can MFEM handle a heterogenous parameter space such as K in this formulation? Any comments are welcomed! Thank you in advance!

Here is a basic working code:

from fenics import *
import mshr

# Sturctured mesh
mesh = RectangleMesh(Point(0.0,0.0),Point(10.0,10.0),100,100)

# Unstructured mesh
# domain = mshr.Polygon([Point(0.0,0.0),Point(10.0,0.0),Point(10.0,10.0),Point(0.0,10.0)])
# domain.set_subdomain(1,mshr.Polygon([Point(0.0,0.0),Point(10.0,0.0),Point(10.0,5.0),Point(0.0,5.0)]))
# res = 200
# mesh = mshr.generate_mesh(domain,res)

ele_h  = FiniteElement("DG",  mesh.ufl_cell(), 0)
ele_q = FiniteElement("BDM", mesh.ufl_cell(), 1)
W = MixedElement([ele_h, ele_q])
W = FunctionSpace(mesh, W)

def K1():
    K_1_val = Constant(1.0)
    return K_1_val

def K2():
    K_2_val = Constant(100.0)
    return K_2_val

class h_boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
h_Boundary = h_boundary()

boundaries = MeshFunction('size_t',mesh,1)
boundaries.set_all(0)
h_Boundary.mark(boundaries,1)

Boundary_head = Expression('sin(x[0]) + sin(x[1])',degree = 5)

class Upper(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 5.0
upper = Upper()

domains = MeshFunction('size_t',mesh,2)
domains.set_all(2)
upper.mark(domains,1)

n = FacetNormal(mesh)
norm = as_vector([n[0], n[1]])

Xf = Function(W)

dx = Measure('dx')(subdomain_data = domains)
ds = Measure('ds')(subdomain_data = boundaries)

U = TrialFunction(W)
V = TestFunction(W)

h, q = split(U)
Ht, Qt = split(V)

a = div(q)*Ht*dx(1) + div(q)*Ht*dx(2) + inner(q,Qt)/K1()*dx(1) + inner(q,Qt)/K2()*dx(2) - h*div(Qt)*dx(1) - h*div(Qt)*dx(2)
L = - Boundary_head*inner(Qt,norm)*ds(1)

A, b = assemble_system(a,L)

solve(A, Xf.vector(), b, 'umfpack')

h, q = Xf.split(True)

h_file = XDMFFile('h.xdmf')
q_file = XDMFFile('q.xdmf')
h_file.write(h)
q_file.write(q)