Eldar,
As you already found out, interpolate_boundary_values can only be used with
a ComponentMask if the element you are using is primitive.
Therefore, the workaround you are using seems to be suitable for the moment.
Since BDM1 is div-conforming, project_boundary_values and
I managed to imply the Neumann conditions on stress only by first
projecting the boundary values to the entire mixed space (affecting
rotations), and then manually removing the constraints associated with
rotations. The code is below if someone is interested.
