Axis treatment for axisymmetric

[trevilo]

As noted in #98, the initial treatment of the axis for axisymmetric flow is not optimal. Specifically, in RHSoperator any source terms (implemented by deriving from the ForcingTerms class) are added after applying the inverse mass matrix. This leads to a requirement that we evaluate 1/radius in the axisymmetric source terms, which is clearly a problem on the axis.

This can be avoided by evaluating the source terms (at least for axisymmetric) prior to applying the inverse mass matrix. This ticket is to do that.