I'm trying to understand how permanent magnets can be modeled in finite element methods. FEMM's manual, here, states that it treats a PM as a solenoid with the same surface current as the PM, and then solves the field from there, which I understand. However, they use the coercivity $H_c$ in place of the magnetization $M$ (I might be dropping factors of $\mu$ here.) to calculate this surface current, i.e., $H_c (\hat{m} \times \hat{n})$.
The argument given is that this is the amount of current you would need to surround the magnet in (in opposite direction) to drive the field to zero (this is the definition of coercivity), which is the same field you get with no current at all. Therefore, the surface current of the original PM can be physically cancelled out with $-H_c (\hat{m} \times \hat{n})$, so this must be it's actual value.
My problem is that when we model a permanent magnet in absence of any external currents, we are at the $H=0$ point of the $B-H$ diagram (i.e., $B = B_r$), and this is not the same physical point. Intuitively I would want to use $B_r$ in place of $M$, modulo some factors of $\mu$ again.
I know that this makes assumptions about linearity of the magnet, but I'm not exactly sure where they creep in. Thanks in advance for any input.