Solve poisson problem with pure Neumann B.C.

Hey,

I was wondering how to solve a poisson problem with pure Neumann conditions?

In Fenics you would use a variational formulation utilizing 2 finite element spaces, one being 0th order: https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/neumann-poisson/demo_neumann-poisson.py.rst

```
# Build function space with Lagrange multiplier
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, P1 * R)
```

But this wouldn't work in jax-fem since there are only order 1 and 2 elements implemented, right?
https://deepmodeling.github.io/jax-fem/_modules/jax_fem/basis.html#get_elements

In firedrake, you can pass a nullspace to the PetSc solver https://www.firedrakeproject.org/solving-interface.html#solving-singular-systems

Since PetSc solver are available this should be possible but I don't really understand how to do this in jax-fem

Would be happy for any tips/recommendations