Hi Petr,
[with NGSolve I have so far only used the first approach below, so please take my suggestions rather as a possible guideline, which might or might not work because I could have overlooked some detail/aspect. I would be interested to hear about your success with any of the other method, if you try them....]
From my (engineering) experience, you only need a rather small penalty parameter to make your problem well-defined.
If the parameter is chosen small enough, the resulting perturbation of the solution *can* be negligible for practical purposes. I am not sure though whether this works well with iterative solvers. You might do a convergence study for decreasing the penalty parameter to get "a feeling" for that.
As a second options, you could try to let the FESpace parameter
"dirichlet_bbbnd" do the trick. It allows you to set Dirichlet BCs on nodes in 3d. How to name the corners of your geometry depends on the geometry backend (CSG, OCC...) that you use.
As a third option, you could manually create a "bitarray" based on fes.FreeDofs(), where you switch the bits corresponding to the displacements of three nodes to be constrained.
I hope this works also in the context of periodic dofs.
A fourth option is to constrain all 6 rigid body modes with one Lagrange multiplier per mode. The space for each multiplier would be a NumberFESpace. This is, however, without any linear algebra tricks, probably not the best approach performance-wise.
Best,
Matthias