Ok, I'm with you and agree that function is probably the way to go. I'd think a function call like extract_... is what you need.
The full matrix method embeds (or used to) the constraint in the row (at the cost of introducing asymmetry), so even without enforce_constraints_exactly() you should recover the constrained value to some accuracy commensurate with your linear solver tolerance.
The constrain_element_residual() should do the trick if it replaces whatever residual value you have with
Re(i) = u(i) - u_constrained
For either dirichlet or hanging-node type constraints... Right?