Let $f=\frac{1}{|x|},x\in\mathbb{R^3}$ and $\Omega=[-b,b]^3$. How to construct a quadrature scheme to solve $$ \int_\Omega f\phi\psi dx\quad ? $$ where $\phi\psi$ is smooth function.

I know there exists a transformation called Duffy transformation which can eliminate the singularity at $x=0$ by using the jacobian. Indeed, let $$\begin{aligned}T&=\{(x_1,x_2,x_3);0\leq x_1\leq 1,0\leq x_2\leq x_1,0\leq x_3\leq x_1\}\\ \hat{T}&=\{(\hat{x}_1,\hat{x}_2,\hat{x}_3);0\leq \hat{x}_1\leq 1,0\leq\hat{x_2}\leq1,0\leq\hat{x_3}\leq 1\}\end{aligned}$$

Define $F:\hat{T}\rightarrow T$ as $$ F(\hat{x}_1,\hat{x}_2,\hat{x}_3)=(\hat{x}_1,\hat{x}_1\hat{x}_2,\hat{x}_1\hat{x}_3) $$ Then $$ \begin{aligned} \int_{T}\frac{1}{|x|}\phi\psi dx&=\int_{\hat{T}}\frac{1}{\sqrt{x_1^2+x_2^2+x_3^2}}\phi\cdot\psi(x)dx =\int_{\hat{T}}\frac{\phi\cdot\psi(\hat{x_1},\hat{x}_1\hat{x}_2,\hat{x}_1\hat{x}_3)}{\sqrt{\hat{x}_1^2+\hat{x}_1^2\hat{x}_2^2+\hat{x}_1^2\hat{x}_3^2}}\hat{x}_1^2d\hat{x}\\ &=\int_\hat{T}\frac{\hat{x}_1\phi\cdot\psi(\hat{x_1},\hat{x}_1\hat{x}_2,\hat{x}_1\hat{x}_3)}{\sqrt{1+\hat{x}_2^2+\hat{x}_3^2}}d\hat{x} \end{aligned} $$ which can be solved by usual quadrature scheme.

On the other hand, someone can also avoid the singularity by omitting a $\varepsilon$ cube which contains the origin.

However, these are not enough for me when I want to study the convergence rate of numerical integration on $\Omega$. I hope the error can be bounded by $O(h)\Vert \phi\Vert_{1,2,\Omega}\Vert \psi\Vert_{1,2,\Omega}$

Is there any other way to integrate a $L^2$ function or a function with singularity at vertex directly with the error bounded by $L^2$ norm?

Thank you very much.