Add matrix-free operator application for nonlinear problems

The existing method jacobian_apply() on the GridOperator does not work correctly for nonlinear problems, as it is not possible to control the linearization point at which to evaluate the Jacobian.

This patch series adds a new method nonlinear_jacobian_apply(x,z,y) that fixes the issue by implementing the operation y += J(x) * z, where J denotes the Jacobian. This also requires new variants of the jacobian_apply_*() methods of the local operator that take the linearization point as an additional parameter. This is documented in the interface document of the local operator interface. Finally, the patch also adds a numerical implementation of the new interface and takes the opportunity to clean up the existing defaultimp.hh header and split it into separate files for the implementations of Jacobian, linear operator application and nonlinear operator application as well as the jacobian-based residual construction.

This will be backported to 2.4, but doing so first requires backporting the concept support from the current dune-common master to the 2.4 release of PDELab.