Zero internal Neumann boundary

Hello all, I have a similar problem with dolfinx 0.5.0 and I can’t seem to apply the method above without throwing an AssertionError. I can’t get my form to compile.

The idea is to have 2D Navier-Stockes with a 1D wall inside the domain. I’d like the pressure to verify \partial_yp^+=\partial_yp^-=0 at the wall but allowing for discontinuities through, as in the picture below.

The relevant part of my current code is largely taken from dokken’s tutorial and looks like this :

boundaries = [(1, lambda x: np.logical_or(top(x),outlet(x))), (2, symmetry), (3, nozzle)]

facet_indices, facet_markers = [], []
fdim = spy.mesh.topology.dim - 1
for (marker, locator) in boundaries:
	facets = dfx.mesh.locate_entities(mesh, fdim, locator)
	facet_indices.append(facets)
	facet_markers.append(np.full_like(facets, marker))
	
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)

face_tag = dfx.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=mesh, subdomain_data=face_tag)
n = ufl.FacetNormal(mesh)
weak=ufl.inner(p.dx(1)('+')*spy.r,s('+')*n[1])*ds(3)

p is a DG type element, which seemed the simplest way to allow for such sharp jumps.

I have looked around this forum but not yet managed to find anything related to this specific problem - i.e. a Neumann or no-flux condition inside the domain.