Force on a Magnetized Body in a Non-Uniform Magnetic Field In Griffith's Electrodynamics, it is shown that when a magnetic moment, $\mathbf{m}$, is subjected to a non-uniform magnetic field, then it experiences a net force given by  $\mathbf{F}=\mathbf{\nabla} (\mathbf{m} \cdot\mathbf{B})$.
Now suppose we have a magnetized body of volume, $V$, with a magnetization $\mathbf{M}\equiv \frac{\mathbf{m}}{V} $. Is it correct to say that the force the body experiences is given by  $\mathbf{F}= \int \mathbf{\nabla} (\mathbf{M} \cdot\mathbf{B}) dV$ ? If yes, would this force be uniform throughout the magnetized body?
 A: Assuming the magnetic field is linear (i.e. we can Taylor expand $\mathbf{B}$ around the origin and truncate all terms of order higher than $1$), then you are correct. The formula is also approximately true for more general $\mathbf{B}$ as long as the dimensions of the extended body are small compared to the region over which $\mathbf{B}$ is approximately linear.
The proof is almost exactly the same as the proof for the force on a magnetic moment $\mathbf{m}$, but we use the definition of magnetization $\mathbf{M}$ to rewrite it as an integral over the volume.
The relevant step for the OP is the last step in the derivation I gave in next section, which I reproduce here:
The definition of magnetic moment in terms of magnetization is $\mathbf{m} = \iiint_V\mathbf{M}(\mathbf{x}')\mathrm{d}^3x'$. From Griffiths' formula for the force on the magnetic moment we have $$\mathbf{F} = \left[\nabla(\mathbf{m}\cdot\mathbf{B}(\mathbf{x}))\right](\mathbf{0})$$
Now we can rewrite $\mathbf{m}$ in terms of $\mathbf{M}$ to obtain the desired formula by the OP:
$$\mathbf{F} = \left[\nabla\left(\iiint_V \mathbf{M}(\mathbf{x}')\mathrm{d}^3x'\cdot\mathbf{B}(\mathbf{x})\right)\right](\mathbf{0}) = \iiint_V \left[\nabla(\mathbf{M}(\mathbf{x}')\cdot\mathbf{B}(\mathbf{x}))\right](\mathbf{0})\mathrm{d}^3x'$$
So the force per unit volume is $\mathbf{f}(\mathbf{x}') = \left[\nabla(\mathbf{M}(\mathbf{x}')\cdot\mathbf{B}(\mathbf{x}))\right](\mathbf{0})$ and if the magnetization is uniform, this force density is also uniform.
Just for completeness, here is a full derivation based on Jackson, Chapter 5. I tried being explicit which variables we differentiate and evaluate against. The steps get a bit tedious, but I wanted to organize it so that all details are explained.

Write $\mathbf{B}(\mathbf{x'})$ as a Taylor expansion around $\mathbf{0}$:
$$\mathbf{B}(\mathbf{x'}) = \mathbf{B}(\mathbf{0}) + \left[\left(\mathbf{x'}\cdot\nabla\right)\mathbf{B}(\mathbf{x})\right](\mathbf{0}) + \ldots \tag{1}\label{1}$$
We assume we can neglect higher order terms (so if $\mathbf{B}$ is linear over space, this is exactly true). Also the $\nabla$ and evaluation operators do not act on $\mathbf{x}'$, but on $\mathbf{x}$.
The total force on the body due to the external magnetic field $\mathbf{B}$ is given (exactly) by
$$\mathbf{F} = \iiint_V \mathbf{J}(\mathbf{x}')\times\mathbf{B}(\mathbf{x}')\mathrm{d}^3x' \tag{2}\label{2}$$
The current $\mathbf{J}$ is assumed to be zero outside of the finite volume of the body. We may take $V$ to be a volume enclosing the body, with some slight empty space included to avoid having discontinuities on the surface of $V$ (this is important in ($\ref{15}$) below). We can include the slight empty space because it does not affect the volume integral for the force.
Combining ($\ref{1}$) and ($\ref{2}$) we obtain:
$$\mathbf{F} = -\mathbf{B}(\mathbf{0})\times\iiint_V \mathbf{J}(\mathbf{x}')\mathrm{d}^3x' + \iiint_V\mathbf{J}(\mathbf{x}')\times\left[\left(\mathbf{x'}\cdot\nabla\right)\mathbf{B}(\mathbf{x})\right](\mathbf{0})\mathrm{d}^3x' \tag{3}\label{3}$$
The first term is zero in the magnetostatic case where $\nabla'\cdot\mathbf{J}(\mathbf{x}') = 0$ (see for example here for a proof).
We have the general vector calculus identities:
$$\nabla(\mathbf{A}\cdot\mathbf{B}) = (\mathbf{A}\cdot\nabla)\mathbf{B} + (\mathbf{B}\cdot\nabla)\mathbf{A} + \mathbf{A}\times(\nabla\times\mathbf{B}) + \mathbf{B}\times(\nabla\times\mathbf{A}) \tag{4}\label{4}$$
$$\nabla\times(\psi\mathbf{A}) = \psi(\nabla\times\mathbf{A}) + (\nabla\psi)\times\mathbf{A} \tag{5}\label{5}$$
Using ($\ref{4}$) with $\mathbf{A} = \mathbf{x}'$ and $\mathbf{B} = \mathbf{B}(\mathbf{x})$ we have:
$$\mathbf{J}(\mathbf{x}')\times\nabla(\mathbf{x}'\cdot\mathbf{B}(\mathbf{x})) = \mathbf{J}(\mathbf{x}')\times\left[(\mathbf{x}'\cdot\nabla)\mathbf{B}(\mathbf{x})\right] \tag{6}\label{6}$$
because $\nabla$ acts on $\mathbf{x}$, not on $\mathbf{x}'$ and because $\nabla\times\mathbf{B} = \mathbf{0}$ (since $\mathbf{B}$ is an external magnetic field, the sources are not inside the volume $V$ and Ampere's Law implies zero curl in that case).
Using ($\ref{5}$) with $\psi = \mathbf{x}'\cdot\mathbf{B}(\mathbf{x})$ and $\mathbf{A} = \mathbf{J}(\mathbf{x}')$ we have:
$$\nabla\times\left[(\mathbf{x}'\cdot\mathbf{B}(\mathbf{x}))\mathbf{J}(\mathbf{x}')\right] = \nabla(\mathbf{x}'\cdot\mathbf{B}(\mathbf{x}))\times\mathbf{J}(\mathbf{x}') \tag{7}\label{7}$$
Combining ($\ref{6}$) and ($\ref{7}$) we obtain:
$$\mathbf{J}(\mathbf{x}')\times\left[\left(\mathbf{x'}\cdot\nabla\right)\mathbf{B}(\mathbf{x})\right] \stackrel{(\ref{6})}{=} \mathbf{J}(\mathbf{x}')\times\nabla(\mathbf{x}'\cdot\mathbf{B}(\mathbf{x})) \stackrel{(\ref{7})}{=} -\nabla\times\left[(\mathbf{x}'\cdot\mathbf{B}(\mathbf{x}))\mathbf{J}(\mathbf{x}')\right] \tag{8}\label{8}$$
Using this in ($\ref{3}$) we get:
$$\mathbf{F} = -\iiint_V \nabla\times\left[(\mathbf{x}'\cdot\mathbf{B}(\mathbf{x}))\mathbf{J}(\mathbf{x}')\right](\mathbf{0})\mathrm{d}^3x' = -\left\{\nabla\times\iiint_V\left[\mathbf{x}'\cdot\mathbf{B}(\mathbf{x})\right]\mathbf{J}(\mathbf{x}')\mathrm{d}^3x'\right\}\bigg\vert_{\mathbf{x}=\mathbf{0}} \tag{9}\label{9}$$
Now we use another vector calculus identity:
$$\mathbf{A}\times(\mathbf{B}\times\mathbf{C}) = (\mathbf{A}\cdot\mathbf{C})\mathbf{B} - (\mathbf{A}\cdot\mathbf{B})\mathbf{C} \tag{10}\label{10}$$
Use ($\ref{10}$) with $\mathbf{A} = \mathbf{B}(\mathbf{x})$ and $\mathbf{B} = \mathbf{x}'$ and $\mathbf{C} = \mathbf{J}(\mathbf{x}')$ to obtain:
$$\mathbf{B}(\mathbf{x})\times(\mathbf{x}'\times\mathbf{J}(\mathbf{x}')) = \left[\mathbf{B}(\mathbf{x})\cdot\mathbf{J}(\mathbf{x}')\right]\mathbf{x}' - \left[\mathbf{x}'\cdot\mathbf{B}(\mathbf{x})\right]\mathbf{J}(\mathbf{x}') \tag{11}\label{11}$$
We want to calculate $$\iiint_V \left[\mathbf{B}(\mathbf{x})\cdot\mathbf{J}(\mathbf{x}')\right]\mathbf{x}'\mathrm{d}^3x' = \sum\limits_{i=1}^3B_i(\mathbf{x})\iiint_V J_i(\mathbf{x}')\mathbf{x}'\mathrm{d}^3x' \tag{12}\label{12}$$
We now have vector calculus identity:
$$\nabla'\cdot(\psi\mathbf{A}) = (\nabla'\psi)\cdot\mathbf{A} + \psi(\nabla'\cdot\mathbf{A}) \tag{13}\label{13}$$
Use ($\ref{13}$) with $\psi = \mathbf{x}'\cdot\hat{\mathbf{e}}_i = x_i'$ and $\mathbf{A} = \mathbf{J}(\mathbf{x}')$:
$$\nabla'\cdot(x_i'\mathbf{J}(\mathbf{x}')) = (\nabla'x_i')\cdot\mathbf{J}(\mathbf{x}') + x_i'(\nabla'\cdot\mathbf{J}(\mathbf{x}')) = J_i(\mathbf{x}') \tag{14}\label{14}$$
Here we used $\nabla'x_i' = \hat{\mathbf{e}}_i$ and $\nabla'\cdot\mathbf{J}(\mathbf{x}') = 0$ (because we work in magnetostatic case).
So using ($\ref{14}$) and then integration by parts we have:
$$\iiint_V J_i(\mathbf{x}')\mathbf{x}'\mathrm{d}^3x' = \iiint_V\nabla'\cdot(x_i'\mathbf{J}(\mathbf{x}'))\mathbf{x}'\mathrm{d}^3x' = \sum\limits_{j=1}^3\iiint_V\nabla'\cdot(x_i'\mathbf{J}(\mathbf{x}'))x_j'\mathrm{d}^3x'\hat{\mathbf{e}}_j $$$$= \sum\limits_{j=1}^3\iiint_V\left[\nabla'\cdot(x_i'\mathbf{J}(\mathbf{x}')x_j') - x_i'\mathbf{J}(\mathbf{x}')\cdot\nabla'x_j'\right]\mathrm{d}^3x'\hat{\mathbf{e}}_j = \sum\limits_{j=1}^3\left[\mathop{\LARGE\unicode{x222f}}_{\partial V}x_i'x_j'\mathbf{J}(\mathbf{x}')\cdot\hat{\mathbf{n}}'da' - \iiint_V x_i'J_j(\mathbf{x}')\mathrm{d}^3x'\right]\hat{\mathbf{e}}_j$$$$=-\iiint_V x_i'\mathbf{J}(\mathbf{x}')\mathrm{d}^3x'\tag{15}\label{15}$$
To move from the second to the third equality, we used ($\ref{13}$) again with $\psi = x_j'$ and $\mathbf{A} = x_i'\mathbf{J}(\mathbf{x}')$. To move from the third to the fourth equality, we used the divergence theorem and the result $\nabla'x_j' = \hat{\mathbf{e}}_j$. To move from the fourth to the fifth equality, we used $\mathbf{J} = \mathbf{0}$ on the boundary of $V$ because we chose $V$ slightly larger than the region where the current  density is nonzero. We also recombined the components of $\mathbf{J}$.
After all this mess, we see that ($\ref{12}$) becomes:
$$\iiint_V \left[\mathbf{B}(\mathbf{x})\cdot\mathbf{J}(\mathbf{x}')\right]\mathbf{x}'\mathrm{d}^3x' = -\sum\limits_{i=1}^3B_i(\mathbf{x})\iiint_V x_i'\mathbf{J}(\mathbf{x}')\mathrm{d}^3x' = -\iiint_V \left[\mathbf{x}'\cdot\mathbf{B}(\mathbf{x})\right]\mathbf{J}(\mathbf{x}')\mathrm{d}^3x' \tag{16}\label{16}$$
Therefore, after taking integrals of ($\ref{11}$) we obtain:
$$\iiint_V \left[\mathbf{x}'\cdot\mathbf{B}(\mathbf{x})\right]\mathbf{J}(\mathbf{x}')\mathrm{d}^3x' = -\frac{1}{2}\mathbf{B}(\mathbf{x})\times\iiint_V\mathbf{x}'\times\mathbf{J}(\mathbf{x}')\mathrm{d}^3x' \tag{17}\label{17}$$
The magnetization is defined by $\mathbf{M}(\mathbf{x}') = \frac{1}{2}\mathbf{x}'\times\mathbf{J}(\mathbf{x}')$ and the magnetic moment is defined by $\mathbf{m} = \iiint_V \mathbf{M}(\mathbf{x}')\mathrm{d}^3x'$ so we conclude from ($\ref{9}$) and ($\ref{17}$):
$$\mathbf{F} = \frac{1}{2}\left\{\nabla\times\left[\mathbf{B}(\mathbf{x})\times\iiint_V\mathbf{x}'\times\mathbf{J}(\mathbf{x}')\mathrm{d}^3x'\right]\right\}\bigg\vert_{\mathbf{x}=\mathbf{0}} = \left\{\nabla\times\left[\mathbf{B}(\mathbf{x})\times\iiint_V\mathbf{M}(\mathbf{x}')\mathrm{d}^3x'\right]\right\}\bigg\vert_{\mathbf{x}=\mathbf{0}} = \left\{\nabla\times(\mathbf{B}(\mathbf{x})\times\mathbf{m})\right\}\big\vert_{\mathbf{x}=\mathbf{0}} \tag{18}\label{18}$$
We have yet another vector calculus identity:
$$\nabla\times(\mathbf{A}\times\mathbf{B}) = \mathbf{A}(\nabla\cdot\mathbf{B}) - \mathbf{B}(\nabla\cdot\mathbf{A}) + (\mathbf{B}\cdot\nabla)\mathbf{A} - (\mathbf{A}\cdot\nabla)\mathbf{B} \tag{19}\label{19}$$
Use ($\ref{19}$) with $\mathbf{A} = \mathbf{B}(\mathbf{x})$ and $\mathbf{B} = \mathbf{m}$. We conclude that:
$$\mathbf{F} = \left[(\mathbf{m}\cdot\nabla)\mathbf{B}(\mathbf{x})\right](\mathbf{0})$$
because $\nabla$ does not act on $\mathbf{m}$ and Gauss' law of magnetism implies $\nabla\cdot\mathbf{B} = 0$.
We now use ($\ref{4}$) with $\mathbf{A} = \mathbf{m}$ and $\mathbf{B} = \mathbf{B}(\mathbf{x})$, and again $\nabla\times\mathbf{B} = \mathbf{0}$ (external field) to obtain the formula from Griffiths for the force on magnetic moment $\mathbf{m}$:
$$\mathbf{F} = \left[\nabla(\mathbf{m}\cdot\mathbf{B}(\mathbf{x}))\right](\mathbf{0}) \tag{20}\label{20}$$
Now we can rewrite $\mathbf{m}$ in terms of $\mathbf{M}$ to obtain the desired formula by the OP:
$$\mathbf{F} = \left[\nabla\left(\iiint_V \mathbf{M}(\mathbf{x}')\mathrm{d}^3x'\cdot\mathbf{B}(\mathbf{x})\right)\right](\mathbf{0}) = \iiint_V \left[\nabla(\mathbf{M}(\mathbf{x}')\cdot\mathbf{B}(\mathbf{x}))\right](\mathbf{0})\mathrm{d}^3x' \tag{21}\label{21}$$
The force density is then given by
$$\mathbf{f}(\mathbf{x}') = \left[\nabla(\mathbf{M}(\mathbf{x}')\cdot\mathbf{B}(\mathbf{x}))\right](\mathbf{0}) \tag{22}\label{22}$$
