Advice for meshing an internal boundary

5

Thank you so much for your help. The code that you provided works as desired. I’m now able to read the mesh and run my original fenics code with the new mesh. It isn’t quite working yet, however, because I still have my original issue about locating the interface between the two rectangles and apply a boundary condition there. I attempted the following.

I attempt to label the interface by adding two lines to the with-block:

with pygmsh.occ.Geometry() as geom:
    ...
    interface = geom.add_surface_loop([domain[1]])
    geom.add_physical(interface,'3')
    ...

However, I don’t think this is correct because the following code does not appear the correctly identify domains/boundaries (solution is zero when it shouldn’t be).

'''
Read mesh and mesh function
'''
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)


'''
Mark domains and boundaries
'''
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1,0)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dS = Measure('dS', domain=mesh, subdomain_data=mf)
dx = Measure('dx', domain=mesh, subdomain_data=mf)
ds = Measure('ds', domain=mesh, subdomain_data=mf)

I’ve tried incorporating information from this thread: How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh) - #5 by dokken but I can’t get things quite right.