Plotting subdomains/internal boundaries

mesh
2

There are several issues with your code. The first one is that phi does not align with the grid, thus using the following code:

will not mark any facets.
This can be verified with the following code snippet (note usage of "dS" as these are interior facets):

intBound = IntBound() 
int_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
intBound.mark(int_boundary, 1)
dS_int = Measure("dS", domain=mesh, subdomain_data=int_boundary) 

print(assemble(Constant(1)*dS_int(0)), assemble(Constant(1)*dS_int(1)))

The second issue is that the dolfin interface to matplotlib does not support plotting facet markers and the following error should be printed in your terminal:

*** Warning: Matplotlib plotting backend does not support mesh function of dim 1. Continuing without plotting...

This can be bypassed by saving to pvd and visualizing in Paraview.

File("facets.pvd") << int_boundary

Next I present a code that resolves your problem:

from dolfin import *
import numpy as np

def phi(x):
    return (x[0] - 0.5)**2 + (x[1]-0.5)**2 - 0.2**2

class Omega(SubDomain):
    def inside(self, x, on_boundary):
        return phi(x) > 0


mesh = UnitSquareMesh(30, 30)
domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

omega = Omega() 
omega.mark(domains, 1)  
dx_int = Measure("dx", domain=mesh, subdomain_data=domains)

print(assemble(Constant(1)*dx_int(0)), assemble(Constant(1)*dx_int(1)))
File("domains.pvd") << domains


# Interior facets between the two domains
int_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
tdim = mesh.topology().dim()
mesh.init(tdim-1, tdim)
facet_to_cell = mesh.topology()(tdim-1, tdim)
num_facets = facet_to_cell.size()
domain_values = domains.array()
facet_values = int_boundary.array()
for facet in range(num_facets):
    # Check if interior
    cells = facet_to_cell(facet) 
    if len(cells == 2):
        # Check if facet is on the boundary between two cells with different markers
        values = np.unique(domain_values[cells])
        if len(values) == 2:
            facet_values[facet] = 1        
    else:
        continue

dS_int = Measure("dS", domain=mesh, subdomain_data=int_boundary) 

print(assemble(Constant(1)*dS_int(0)), assemble(Constant(1)*dS_int(1)))
File("facets.pvd") << int_boundary
2 Likes