Thanks a lot for detailed response! I will need to adapt this to my actual problem. Like the normal vector, I guess we can work with any solution variable?
For now, could you please have a look again at this which is running into error on my end:
# Create custom integration measure for one-sided integrals
ds = ufl.Measure("ds", domain=mesh, subdomain_data=[
(8, np.asarray(integration_entities, dtype=np.int32))])