I've been meaning to ask about this API - is there a reason for expressing Neumann BCs as a higher order tensor?  It's only the normal component of the specified input that actually gets matched, right?

Something else that might be worth discussing here is fourth-order problems.  You have a *pair* of BCs at each boundary: either a Dirichlet condition on u or a Neumann condition on something-like-Laplacian(u), plus either a Dirichlet condition on dn(u) or a Neumann condition on dn(something-like-Laplacian(u)).  If we're having to rip up the BC handling anyway it might be worth trying to keep the design compatible with higher order PDEs.