Weak Laplacians on disks?

It would be nice to support weak laplacians. For example, it leads to a nice way of achieving Neumann condition via natural conditions: we would use Weighted(Zernike(1)) combined with the harmonic polynomials real(z^m) and imag(z^m) to pick up the rest of the polynomials. (see the new Neumann p-FEM test for an example in 1D) But it seems like there are a few options:

1. Support gradients of Zernike. But how do we incorporate the rotational invariance? And in a way so dot(∇, ∇*Weighted(Zernike(1))) still gives a diagonal operator?
2. Support tensor calculus a la Vasil et al. If I recall correctly this works on ∇_r and ∇_θ . This will successfully capture rotational invariance, but technically leaves the world of polynomials. That is, we can't directly view it as a vector orthogonal polynomial basis. Though perhaps (1) can be built from these.
3. Do weak Laplacian's directly without constructing gradients. In the case of m = 0 this seems straightforward via:

Start with the new Neumann p-FEM example in test/test_odes.jl and review the three approaches described: gradients of Zernike, tensor calculus, or direct weak Laplacians. Compare how each handles rotational invariance and the stated diagonal-operator goal. Done requires choosing and specifying an implementation path for weak Laplacians on disks.