Missing the marking for interface elements between subdomains

Hi everyone! I am trying to mark my domain in a mesh function with half left being marked as 1 and half right marked as 2. The marking is totally fine but it misses out the elements to mark which are on the subdomain interfaces.
What I want is to include those elements whose volume is more than half already inside the left subdomain (x[0]<=0.5) to be marked as 1 and the one with more than half of there volume inside right domain (x[0]>0.5) to be marked as 2. In this way I want to avoid any of the element to stay marked as 0. Any help would be appreciated. Thanks in advance!

Link to mesh file: https://drive.google.com/file/d/1fB4146_klPikr1RsDE1XWzsIP-pDg8tu/view?usp=sharing

#creating mesh and its function

mesh = df.Mesh('input_mesh.xml')
mf = df.MeshFunction("size_t",mesh, mesh.topology().dim(),0)

LeftDomain = df.CompiledSubDomain(' x[0] <= 0.5 ')
RightDomain = df.CompiledSubDomain('x[0] > 0.5 ')

# marking the domains
LeftDomain.mark(mf, 1)
RightDomain.mark(mf, 2)

#output for visualizsing the marks 
mxf = df.XDMFFile("labels_output/mesh_function_labels.xdmf" )

For an element to be marked as belonging to subdomain, all of its vertices must meet the geometric criterion. Thus, elements with vertices on both sides of x_0 = 0.5 are not in either subdomain. You can flag all elements by making the subdomains overlap near the interface, e.g.,

h_max = mesh.hmax()
LeftDomain = df.CompiledSubDomain(' x[0] < 0.5 + DOLFIN_EPS + '+str(h_max/2))
RightDomain = df.CompiledSubDomain('x[0] > 0.5 - DOLFIN_EPS - '+str(h_max/2))

for a general unstructured mesh (where the size of overlap must depend on element size), or just

LeftDomain = df.CompiledSubDomain(' x[0] < 0.5 + DOLFIN_EPS')
RightDomain = df.CompiledSubDomain('x[0] > 0.5 - DOLFIN_EPS')

if the interface is built directly into the mesh.


Thanks and so nice of you @kamensky ! :slight_smile: