Based on @MiroK’s answer here, the "+" cell for a facet is the one for which the subdomain MeshFunction has a higher value, so you want to use "+" instead of "-" to integrate flux using the function restricted to subdomain 1. Trying this, int_boundary_val is a small nonzero value that converges to zero as the mesh is refined. (You can only expect exact conservation over \Gamma\subset\partial\Omega if there exists a test function whose trace on \partial\Omega is the characteristic function on \Gamma. With continuous basis functions, this is only possible if \Gamma = \partial\Omega.)
2 Likes