Implement hp-adaptivity based on local element reference solutions

[solin]

This is very important as for elliptic (and perhaps other) problems
this would work almost exactly as hp-adaptivity based on global
reference solutions, but it would be much faster.