Recall that Stokes flow in $\mathbb{R}^3$ is the solution $(\textbf{u}, p)$ of
$$\begin{align} -\nabla p + \mu \Delta \textbf{u} &= \textbf{f} \\ \nabla \cdot \textbf{u} &= 0, \end{align} \tag{1}$$
together with the boundary conditions at infinity, which are $\|\textbf{u}\|, \left|p\right| \to 0$ as $\|\textbf{x}\| \to \infty$.
If we let $\textbf{f} = -\textbf{f}_0 \delta(\textbf{x})$ (where $\textbf{f}_0$ is a constant vector and $\delta(\textbf{x})$ is the delta distribution), then $(\textbf{u}, p)$ is the fundamental solution, written as
$$\begin{align} \textbf{u}(\textbf{x}) &= \frac{1}{8 \pi \mu} \left( \frac{\textbf{f}_0}{r} + \frac{\textbf{x} \textbf{x}^T \textbf{f}_0}{r^3} \right) \\ p(\textbf{x}) &= \frac{\textbf{x}^T \textbf{f}_0}{4 \pi r^3}. \end{align} \tag{2}$$
I understand the derivation of the fundamental solution $(\textbf{u}, p)$ (e.g. in Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow (1992), or in two related questions with excellent answers by Willie Wong [1] [2]). What I would like to understand now is how to verify that $(\textbf{u}, p)$ is the fundamental solution. Therefore my question is this: With $(\textbf{u}, p)$ defined as in Equation (2), and the Stokes equation in $\mathbb{R}^3$ defined as in Equation (1), how do I "plug in $(\textbf{u}, p)$," do some calculations, and conclude that $\textbf{f} = -\textbf{f}_0 \delta(\textbf{x})$?
Note: A colleague and I have tried the calculations several times, without success. We are familiar with the identity $\Delta(\frac{1}{r}) = -4 \pi \delta(\textbf{x})$ and its derivation, and with the need to consider things in terms of distributions. We are interested in this question because it is similar to (and perhaps more tractable than) a related question for Stokes equation in $\mathbb{R}^2 \times [0, \infty)$ (where the fundamental solution is given e.g. in Blake & Chwang, "Fundamental singularities of viscous flow" (1974)).