If you want to avoid issues with inacurate normals, have a look at: Setting scalar function for BDM boundary condition with controlled flux [mixed Poisson] - #9 by dokken which can evaluate the FacetNormals exactly for every cell.
As you have not provided a new minimal example to illustrate how you currently want to use the fluxes with another program, I cannot give any further guidance.