I'd like to experiment with setting up a simulation in which the jump condition on an internal boundary is evaluated then, as a secondary step, exported to an internal boundary mesh between the two domains across which the jump occurs. I've tried including something to the tune of
domain1_submesh = SubMesh(mesh, subdomains, 1)
boundaries_domain1 = FacetFunction("size_t", domain1_submesh)
boundaries_domain1.set_all(0)
domain1_vmap = domain1_submesh.data().array('parent_vertex_indices', 0)
boundary_domain1 = boundaries.array()[domain1_vmap]
V1 = FunctionSpace(cell_submesh, 'DG', 2)
u0_domain1 = interpolate(u0, V1)
domain2_submesh = SubMesh(mesh, subdomains, 0)
boundaries_domain2 = FacetFunction("size_t", domain2_submesh)
boundaries_domain2.set_all(0)
domain2_vmap = domain2_submesh.data().array('parent_vertex_indices', 0)
boundary_domain2 = boundaries.array()[domain2_vmap]
V2 = FunctionSpace(domain2_submesh, 'DG', 2)
u0_domain2= interpolate(u0, V2)
bmesh = BoundaryMesh(domain1_submesh, 'exterior')
Vb = FunctionSpace(bmesh, 'DG', 2)
u0_b = interpolate(u0_domain1 - u0_domain2, Vb)
but I keep coming up with incorrect solutions. I'm using a DG method on to solve for u0 prior to trying to interpolate this function. I'm really after finding the scalar value of the jump between the two domains on bmesh. Any help and/or suggestions would be greatly appreciated! Thanks!
