Clarify and extend BM expansion order

[bobmyhill]

Currently, the code is written for BM 3rd order expansion, with a second or third order expansion for the shear modulus. This is not obvious in the user code (although it is in the documentation)

Would it be better to define two parameters, and then have two separate flags: self.shear_order and self.bulk_order, to clarify what is actually being done in the calculations?

[bobmyhill]

The easiest way to do this would be to rewrite the PVK expressions in BM up to fourth order, then have a flag which sets

• K'' = -1./K * (K'(K' - 7.) +143/9) for BM3 and BM2
• K' = 4 for BM2

Different functions would be more efficient though.

[bobmyhill]

Addressed more simply by renaming functions (#649)