I guess a more robust approach would be to push the facet normal forward to the physical domain with Piola map, as shown in: Setting boundary conditions on H(div) spaces defined on general meshes - #2 by dokken
Let’s see if I have further time to look at this soon.