Assume we want to solve the Poisson equation $$ \Delta u = f $$ with free (Neumann) boundary conditions. So, the right hand side function $f$ must fulfill the compatibility condition to integrate to zero. When using linear finite elements and a uniform refined mesh, the rhs vector has also nodal average of zero. So an iterative solver will find the solution. But this does not hold anymore for adaptively refined mesh. Thus, on the level of linear algebra the rhs vector is not orthogonal to the null space of the operator and we cannot assume the solver to find a solution in this case. Is it still valid to project out the null space from the rhs vector? For the constant null space this is the same as the substruction of the rhs nodal average from each vector entry. In this case, is there any interpretation on the finite element space? For me, this does not seem to be correct.
Remove the discrete null space, solve the discrete problem, then postprocess (shift by a constant) to correct for the difference between the continuum null space and discrete null space. You can compute the difference by integrating the discrete null space using quadrature.
The correct compatibility condition is that $\int_\Omega f(x) \; dx = 0$. If you are on a uniform mesh, this means that the average of the right hand side vector entries must be zero, but on a non-uniform mesh (such as an adaptively refined one) that is no longer true. You need to ensure the analogue to this for your right hand side vector which means that you'll need to weigh the entries of your right hand side vector appropriately when taking their average.