Higher order RT elements require more delicate mechanism to mark boundary indices.
    
Working:
2D: RT0 x P0, RT1 x P1, BDM1 x P0, BDM2 x P1
3D: RT0 x P0, RT1 x P1
    
Not working:
2D: RT2 x P2 - solver does not converge.
3D: BDM1 x P0 - no interpolation provided by dune-localfunctions.