I don't understand how it could work on the small grids...

Why do problems have this method anyway? It was obviously never used...

Needs refactoring!

The previous version was too complicated and even wrong as it did not
allow to configure e.g. solvers with doubles and meshes with floats.

Implementation of TimeIndependentPDESolver.

- separate data types for input arrays of binary operations
  (necessary for the equality comparison of two arrays of different
  type)
- added separate template parameter for ResultType, which defaults to
  the RealType
  (the goal is to allow sum of bool vectors into int instead of bool)

This is particularly useful for debugging long simulations which don't
fail soon.