I didn't read through the whole post on the issue. But I found some similar problem for Hex8 elements.

If you checked the attached the picture. This is a local mesh I get from libMesh AMR. There, both A and B are hanging nodes, while A depends on B. So when I do both dofmap.enforce_constraints_exactly(), A is interpolated from B and C, while B as a hanging node has value 0. Therefore, A gets a wrong value.

I am guessing that libMesh level one condition is enforced only on surface, i.e. the level difference of two elements that share surfaces can not exceed one. If this is true, it could be a potential problem for hanging nodes.

Is my conjecture right? How can I solve this problem?

Thank you very much!

-Shun