Previously, we computed the difference ourselves and then called
the norm object to get the norm of that difference.  The change
hence doesn't make any difference in almost all cases except...
In the implementation of the 2-Lagrange-Multiplier method for the
Richards equation, the iteration variables are Robin boundary values
on the skeleton of the decomposition.  However for computing the
errors and convergence rates I want to use the energy norm of the
subdomain solutions.  Hence given Robin traces, my norm object
(MultiDomainPhysicalEnergyNorm) solves the subdomain problems
and computes the norm of that.  But since it is a nonlinear problem
calling diff is not the same as calling operator() for a difference.

[[Imported from SVN: r6330]]