@dokken
Thanks for the suggestion, I somehow did not think to try the dolfinx.fem.Expression, I tried something like this
WW = FunctionSpace(mesh, ("DG", 2))
divuu = Function(WW)
u_expr = dolfinx.fem.Expression(ufl.div(u), WW.element.interpolation_points())
divuu.interpolate(u_expr)
divu=divuu.x.array
Plotting shows some points in the boundary where the divu is not zero though.
@nate Thanks for bringing those to my attention! To be honest I did not know how far the rabbit hole went.
Now my problem is slightly different. I want to calculate the Magnetic Field in some 3D geometry. Now there is some experimental input for the boundary. Long story short the boundaries have some sort of uncertainty in their values and by solving the 3D Poisson equation -\nabla ^2\mathbf{B}=0 I get a family of solutions and I need to keep the physical one, i.e. the one for which \nabla \cdot \mathbf{B}=0. Introducing a Lagrange multiplier to impose the constraint you are left with a system of PDEs that are similar to the Stokes ones, .i.e. the magnetic field is the velocity and the Lagrange multiplier is the pressure. The BC are different though, so since I used fenics a lot in the past I tried to solve my problem with fenicsx using the Poisson and stokes tutorial as a basis.
Seems like I stumbled in a non trivial problem, so thanks again for the help!