Simulate the magnetic field of magnet I'd like to numerically simulate the magnetic field due to a cylindrical rare earth magnet of known dimensions at a bunch of points $\vec{r}=(x,y,z)$ from its origin. My goal is to be as accurate as possible up to some global normalization which ultimately doesn't matter for my application.
My approach is to simply assume the magnet is made of a continuum of dipoles and do the relevant integral:
$$
\frac{\mu_0 m}{4\pi}\int\int\int \frac{3\left(\hat{m}\cdot\frac{\vec{r}-\vec{r}'}{||\vec{r}-\vec{r}'||}\right)\frac{\vec{r}-\vec{r}'}{||\vec{r}-\vec{r}'||}-\hat{m}}{||\vec{r}-\vec{r}'||^3}d\vec{r}'
$$
where the dipole moment norm $m=||\vec{m}||$ is the aforementioned global normalization.
Performing this integral is easy and it works etc., my question is about the physics: am I leaving anything out? Will a commercial solver include any physics that I haven't? In general, what will be my biggest source of error (except that I don't know $m$)?
Thanks!
 A: First of all, you need to know what $\vec{M}$ is, which is a problem! Secondly, you are doing a volume integral, but if $m$ is constant, then from the point of view of the microscopic currents that create the dipoles, only the boundary or surface currents actually contribute (if you draw a lot of directed square loops on a piece of paper, you will see that currents from neighbouring loops cancel, and only one big outer loop remains).
So you can do better by replacing the volume integral over magnetization with a Biot-savart type integral over the "surface current density". In general, this bound current density is given by $J_b = \vec\nabla \times \vec{M}$ but on the surface you get a surface current density $K_b = \vec{M} \times \vec{n}$, where $\vec{n}$ is the surface normal. Assuming your cylindrical magnet is magnetized length-wise, it will have the same field as a solenoid with the aforementioned surface current, and you can do much analytically even.
Anyway, I was wondering the same thing as you, and although I haven't gotten around the numeric part of the problem, I did a lot of digging. I'll copy what I found in my research here.
I think the answer is that you can't directly measure $B$ inside the material, but you can figure out $H$ (it depends on free currents). So you need to know how $B$ depends on $H$ for calculation purposes.  Rare earth magnets have what's called linear demagnetization curves (the demagnetization curve is the second quadrant of the B-H curve, where a negative H field is applied, and is what the manufacturers provide -- the normal part is linear). For such permanent magnets, we can linearize the $B-H$ curve in the second quadrant, and write
$H(B) = \nu B - H_c$
where $\nu = \frac{H_c}{B_r}$ is the reluctivity (inverse of permeability) and $H_c$ is the coercivity, i.e., the value of $H$ required to drive $B$ to zero. Note that if this were a linear material without hysteresis (e.g. a soft iron core), you would just have
$H(B) = \nu B$.
If you apply Ampere's law to this equation you get $\nabla \times H = J_f$. For our permanent magnet, the extra term $-\nabla \times H_c$ in the first equation is equivalent to adding an extra bound current term to the second equation. The bound current only has support on the surface, so we treat the permanent magnet like a soft iron core surrounded by a surface current of value $-H_c \hat {n}\times \hat{m}$.
Incidently, this is the way a lot of permanent magnet solvers work, and in particular, is how FEMM works (see Item 6 in their faq.)
I found the following documents useful:
http://www.permagsoft.com/english/assets/applets/DemagnetisingCurve.pdf
http://www.permagsoft.com/english/assets/applets/CompEngl.pdf
