Currently, mixed discretizations, at least with the scalar primal and dual
functions for Helmholtz-type operators are not working due to wrong hardcoded
loop boundaries.

This is now fixed by using the number of shape functions.

These can be useful for mixed discretization of Helmholtz3D
type of integral equations.

Overall, the observed error is comparable to the standard discretizations.

Note, however, that it might be useful (TODO) adding unit tests for
individual matrix element assemblies.