There are plenty of posts on how to read MeshValueCollection objects, see for instance:How to read a physical object from a xdmf file - #7 by dokken
It is not very different in C++.
I guess n is defined by the surface of the domain, i.e. on the facets where the interface between marked subdomain and other domains are? Then you can use FacetNormal(mesh) + the dS measure, see: Integrating over an interior surface - #4 by MiroK
Do you want to solve a PDE where u is the unknown?
If so, what kind of PDE are you trying to solve?