This commit add the following feature:
When evaluating the residual value, the rate terms of xi is now integrated through each Gauss point.
This gives more accurate residual values.

The crash of linear solver is not improved by this change.