Need EquationSystems to reinit constraints

This whole section of code is hideous and needs libMesh-based
refactoring, but at least it works now.
parsed_fem_dirichlet boundaries now work

So let's get the regression test set down.

This is a horribly understabilized case, but that should be fine for a
regression test.  I don't want to use a finer mesh because I'm afraid
that the Dirichlet<->Neumann transition might be borderline somewhere
and rounding error flipping the sign of u=+/-epsilon might then
trigger a larger change.