This question is not entirely similar to the question here. Please read this question and the reader will see it is obviously not the same.
$\mathbf{M'}$ is a continuous vector field in volume $V'$ and $P$ be any point on the surface of $V'$ with position vector $\mathbf {r}$
PART I:
Consider the expression:
$\displaystyle \iiint_{V'} \mathbf{M(r')}.\nabla' \left( \dfrac{1}{\left| \mathbf{r}-\mathbf{r'} \right|} \right) dV'$
Take the origin of our coordinate system at $P$ (see the diagram) and write $dV'$ as spherical volume element. Then the above expression can be written as:
$\displaystyle \iiint_{V'} \mathbf{M(r')}.\dfrac{\hat{r'}}{{r'}^2} {r'}^2\ \sin\theta\ d\theta\ d\phi\ dr' = \iiint_{V'} \left[ \mathbf{M(r')}.{\hat{r'}} \right] \sin\theta\ d\theta\ d\phi\ dr' \tag{1}$
The integrand is defined everywhere except at point $P$ where $r'=0$ and the integrand is $\frac{0}{0}$.
Since $\mathbf{M(r')}$ is finite everywhere, there is no blowing up of the integrand at any point.
Therefore we can directly integrate equation $(1)$ just like ordinary integrals.
PART II:
Using the vector identity $\nabla.(\psi \mathbf{A})=\mathbf{A}.(\nabla \psi)+\psi (\nabla.\mathbf{A})$:
$\displaystyle \iiint_{V'} \left[ \nabla' . \left( \dfrac{\mathbf{M'}}{\left| \mathbf{r}-\mathbf{r'} \right|} \right) \right] dV' = \iiint_{V'} \mathbf{M'}.\nabla' \left( \dfrac{1}{\left| \mathbf{r}-\mathbf{r'} \right|} \right) dV' + \iiint_{V'} \dfrac{\nabla' . \mathbf{M'}}{\left| \mathbf{r}-\mathbf{r'} \right|} dV'\tag{2}$
Now for simplicity, let's take the origin of our coordinate system at $P$ (see the diagram). Thus equation $(2)$ becomes:
$\displaystyle \iiint_{V'} \left[ \nabla' . \left( \dfrac{\mathbf{M'}}{r'} \right) \right] dV' = \iiint_{V'} \mathbf{M'}.\nabla' \left( \dfrac{1}{r'} \right) dV' + \iiint_{V'} \dfrac{\nabla' . \mathbf{M'}}{r'} dV' \tag{3}$
Now by writing $dV'$ as spherical volume element, equation $(3)$ becomes:
\begin{align} \iiint_{V'} \left[ \nabla' . \left( \dfrac{\mathbf{M'}}{r'} \right) \right] {r'}^2\ \sin\theta\ d\theta\ d\phi\ dr &=\iiint_{V'} \mathbf{M'}.\nabla' \left( \dfrac{1}{r'} \right) {r'}^2\ \sin\theta\ d\theta\ d\phi\ dr'\\ &+\iiint_{V'} \dfrac{\nabla' . \mathbf{M'}}{r'} {r'}^2\ \sin\theta\ d\theta\ d\phi\ dr'\\ &= \iiint_{V'} (\mathbf{M'}. \hat{r'}) \sin\theta\ d\theta\ d\phi\ dr' \\ &+ \iiint_{V'} (\nabla' . \mathbf{M'})\ r'\ \sin\theta\ d\theta\ d\phi\ dr' \tag{4} \end{align}
In both terms:
The integrands are defined everywhere except at point $P$ where $r'=0$ and the integrand $\frac{0}{0}$.
Since $\mathbf{M(r')}$ is finite everywhere, there is no blowing up of the integrands at any point.
Therefore we can directly integrate both terms in the $RHS$ of equation $(4)$ just like ordinary integrals.
Therefore we can directly integrate $LHS$ of equation $(4)$ just like ordinary integrals.
Now my question is: Is Gauss divergence theorem applicable to $LHS$ of equation $(4)$ ?
If it is applicable, then, since point $P$ lies on the surface, there would be a singularity in the equation:
$\unicode{x222F}_{S'} \left[ \left(\dfrac{\mathbf{M'}}{\left| \mathbf{r}-\mathbf{r'} \right|} \right) . \hat{n} \right] dS'$
How to deal with it?