I have a 3D unit vector field $v(x,y,z)$ in a domain $\Omega$ with some boundary constraints. I would like to numerically integrate the p-norm of its "helicity", where helicity is defined as $$h(x,y,z)=v \cdot (\nabla \times v)$$ and its p-norm is $$E_p[v]=(\int_{\Omega} h(x,y,z)^p dV)^{\frac{1}{p}}$$ Naturally, the unit vector field with boundary constraints will contain singularities, around which helicity might diverge. I am looking for a numerical method of evaluating $E_{\infty}[v]$ when I don't explicitly know the location of the singularity.
$\textbf{For example}$, let rotation matrix $R(a)$ be a rotation along the x-axis by an angle 'a'. Let the unit vector field $v$ on a unit sphere $\Omega$ be $$v(x,y,z)=R(a) \frac{[x,y,0]^T}{x^2+y^2}$$ The helicity of this field analytically is $h(x,y,z)=\frac{x \; sin(a)}{x^2+y^2}$. Thus if $a= n\pi$, for integer n, then $E_{\infty}[v] = 0$. Otherwise, $E_{\infty}[v] = \infty$.
How can I numerically integrate $E_{\infty}[v]$ such that under mesh refinement, the numerical result approaches the analytic results? I have already tried naive FEM with no success which I think is due to inaccuracy in the neighborhood of the singularity. What numerical integration methods can I use to evaluate an integral that may or may not diverge, over a singularity whose location is not known?