How to integrate on the surface of a subdomain (without making the surface its own subdomain)

[dokken]

You issue is in the definition of the MeshFunction as you define it as a function over cells.

@jefromm

sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim())

You should change it to

sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
Right().mark(sub_domains,1)
File("mf.pvd") << sub_domains

and you can visualize the markers on facets in for instance Paraview.
Using this marker will get you something close to 1 for

@jefromm

R_subdomain_s =f*ds_custom
print(assemble(R_subdomain_s))

However, I disagree with the statement that

@jefromm

R_subdomain_S =f*dS_custom
print(assemble(R_subdomain_S))

should be 1, as dS is the integration measure over internal facets (those that are not on the boundary).
As no internal markers has the value 1, the integral will be 0