Running it with debug mode, there is a crash, so any potential solution output is completely unreliable.

edit: This might be my own fault, let me get back to this later

Tracking the error backwards, I see we end up with NAN values in the matrix itself, a block of these on the last DOFs/equations in the system matrix.

As far as I can tell, that should reflect the DOFs in Node 17 (since 18-20 are fully prescribed).
I'm not sure what's the cause yet.

A few details to sort out,

The domain type is deprecated, and is (almost) not used.  Instead, the necessary DOFs are generated by what the connected elements need, node-per-node.

You can also specify directly what dofs nodes should contain, by specifying the "dofidmask" (you can do this for any node).

Nothing is connected to node 20 right now, so no DOFs are created there. You specify that 17 is connected to it, but there are no DOFs to connect to.

Solution would be to add "dofidmask" to node the master nodes as well (or every node, if you want, it makes no difference in this case)

But (!) what are you trying to achieve here?
Having a slave node to a prescribed node.. with no other elements connected to it, is pointless. It's the identical effect as just prescribing the slave node directly.

That used to be what the domain type did, but in mixed problems (like some shell elements, some solid elements) or incompressible problems (linear pressure, quadratic velocities), XFEM, contact problems, etc. etc.  you would have different dofids in each node.

The "eigen_beam3d" file uses dofcondense in a very pointless way, which is why you needed to do unnecessary things like this. So it's not the ideal problem to translate into slave nodes.
Instead, eigen_beam3d should have just used BCs to achieve this with the existing nodes.

eigen_beam3d.out
eigen vibration analysis of simple suported beam
EigenValueDynamic nroot 4 rtolv 1.e-6 nmodules 1
errorcheck
domain 3dShell
OutputManager tstep_all dofman_all element_all
ndofman 18 nelem 16 ncrosssect 1 nmat 1 nbc 1 nic 0 nltf 1 nset 2
node 1 coords 3 0.   0.    0.00
node 2 coords 3 0.   0.0   0.25
node 3 coords 3 0.   0.0   0.50
node 4 coords 3 0.0  0.0   0.75
node 5 coords 3 0.   0.0   1.00
node 6 coords 3 0.   0.0   1.25
node 7 coords 3 0.   0.0   1.50
node 8 coords 3 0.0  0.0   1.75
node 9 coords 3 0.   0.0   2.00
node 10 coords 3 0.   0.0   2.25
node 11 coords 3 0.   0.0   2.50
node 12 coords 3 0.0  0.0   2.75
node 13 coords 3 0.   0.0   3.00
node 14 coords 3 0.   0.0   3.25
node 15 coords 3 0.   0.0   3.50
node 16 coords 3 0.0  0.0   3.75
node 17 coords 3 0.   0.0   4.00
node 18 coords 3 1.0  0.0   0.00
#
Beam3d 1 nodes 2 1 2 refNode 18
Beam3d 2 nodes 2 2 3 refNode 18
Beam3d 3 nodes 2 3 4 refNode 18
Beam3d 4 nodes 2 4 5 refNode 18
Beam3d 5 nodes 2 5 6 refNode 18
Beam3d 6 nodes 2 6 7 refNode 18
Beam3d 7 nodes 2 7 8 refNode 18
Beam3d 8 nodes 2 8 9 refNode 18
Beam3d 9 nodes 2 9 10 refNode 18
Beam3d 10 nodes 2 10 11 refNode 18
Beam3d 11 nodes 2 11 12 refNode 18
Beam3d 12 nodes 2 12 13 refNode 18
Beam3d 13 nodes 2 13 14 refNode 18
Beam3d 14 nodes 2 14 15 refNode 18
Beam3d 15 nodes 2 15 16 refNode 18
Beam3d 16 nodes 2 16 17 refNode 18
#
#
SimpleCS 1 area 0.06 Iy 0.00045 Iz 0.0002 Ik 0.000498461  beamShearCoeff 1.e60 material 1 set 1
IsoLE 1 d 25.0 E 25.e6 n 0.15 tAlpha 1.2e-5
BoundaryCondition 1 loadTimeFunction 1 dofs 4 1 2 3 4 values 4 0. 0. 0. 0. set 2
ConstantFunction 1 f(t) 1.
Set 1 elementranges {(1 16)}
Set 2 nodes 2 1 17
#
#
#%BEGIN_CHECK% tolerance 1.e-2
## check eigen values
#EIGVAL tStep 1 EigNum 1 value 1.28049596e+03
#EIGVAL tStep 1 EigNum 2 value 2.85378786e+03
#EIGVAL tStep 1 EigNum 3 value 2.10785640e+04
#%END_CHECK%

This is how the file should have been. It's simpler, and clearer (more direct).
You don't need dofstocondense or slave nodes in order to have momentless boundary conditions.

Dea bp,
thanks for the code. Unfortunately we cannot find the file bctracker.h in the repo.

We will test it as soon as we can compile the code undeer windows (see here)

Thanks bp. Now I can compile and beam3d seems to work properly, but not with dostributed loads. Consider for example this input:

beam3dsubsoil02.out
test
LinearStatic nsteps 1
domain 3dShell
OutputManager tstep_all dofman_all element_all
ndofman 4 nelem 2 ncrosssect 2 nmat 2 nbc 2 nic 0 nltf 1 nset 3
node 1 coords 3 0.  0.  0.
node 2 coords 3 1800  0.  0.
node 3 coords 3 3600 0. 0.
node 4 coords 3 3600 1.  0.
#
Beam3d 1 nodes 2 1 2 refNode 4 subsoilmat 2 mat 1
Beam3d 2 nodes 2 2 3 refNode 4 subsoilmat 2 mat 1
#
SimpleCS 1 area 80000 Iy 1066666666.6666665 Iz 266666666.66666663 Ik 160000000.0 beamShearCoeff 1.e30 material 1 set 1
SimpleCS 2 
IsoLE 1 d 1. E 2.e5 n 0.15 tAlpha 1.2e-5
winkler 2 d 0.0 c1 6 0 0 20. 0 0 0 global 1
BoundaryCondition 1 loadTimeFunction 1 dofs 4 1 2  4   6 values 4 0.  0. 0. 0. set 2
#ConstantEdgeLoad 2 loadTimeFunction 1 dofs 6 1 2 3 4 5 6 Components 6 0.0 0.0 20000000.0 0.0 0.0 0.0 loadType 3 set 3
#NodalLoad 2 loadTimeFunction 1  Components 6 0.0 0.0 1.e5 0.0 0.0 0.0 
ConstantEdgeLoad 2 loadTimeFunction 1 Components 1 -1 DOFs 1 3 set 3
PeakFunction 1 t 1.0 f(t) 1.
Set 1 elementranges {(1 2)}
Set 2 nodes 3 1 2 3
Set 3 elementranges {(1 2)}

ps.
I would be surely useful to have the same modification on beam2d, but also in triangular shells tr_shell01 (tria1platesubsoil) and mitc4shell (or quad1mindlinshell?) (quad1platesubsoil) if you can and if it is possible.

Thanks!

Just a final question on that: if I want to use shell elements as elastic foundation, can I overlap the elements tria1platesubsoil or quad1platesubsoil to the original shells?

Thanks for the clarification.

EDIT: from official GIT repo:
git cherry-pick b9bf8b8f13713700f7e2dcfa0134defc6bec7e31

fatal: bad object b9bf8b8f13713700f7e2dcfa0134defc6bec7e31   <<<<<<<<<<<<<<<<<<<<<

is it an error from my side?

EDIT2: repo issue solved.

Going back to the beam3d, we noticed that you use Beam3d::computeNmatrixAt to get the stiffness contribution from elastic soil. Can you clarify the formulation used? We need to correct the bem module to extract the proper force diagrams alogn the beam.

In winkler material, I noticed there's no dependence on width or height of the beam. Which are the units for the parameter c1? It seems they are force/Length^2 ...

Anyway, the input parameters are not clear. Can you help on that? Winkler subsoil constant usually has force/Length^3 ...

The Winkler  model assumes that the foundation resists to lateral displacement. There is no shear stiffness of the foundation assumed, so surrounding foundation is unaffected by the deformation under the beam.
The resisting pressure p = Kw, where K is modulus of foundation [N/m^3] and w is corresponding lateral deflection. The differential dx of the beam with area dA=dx*b (b being the resisting width of the beam) acts like a linear spring with stiffness c = p*dx*b/w = K*w*dx*b/w = K*b*dx.
The force in the spring is c*w and according to the principle of virtual work the total work on the beam done by spring force on corresponding displacement w is E = 1/2 \int_0^l c*w^2 dx = 1/2 \int_0^l K*b*w^2 dx, thus after substituting the approximation for lateral displacement the stiffness matrix of the foundation is
K = \int_0^l N_w^T K*b N_w dx, where N_w is lateral displacement interpolation matrix.

The oofem on input expects the foundation stiffness multiplied by the corresponding resisting width K*b [N/m^2], representing the stiffness per length. The formulation is however extended that additional spring stiffness can be applied to arbitrary DOF on element, for example in longitudinal direction or distributed rotational spring stiffness.
Hope that it will help.

Thank you for reporting the problem.
I managed to fix the bug in reaction forces evaluation for beam3d with active subsoil.
Regarding the output of contact pressure: currently there is no support for this. However, it is relatively straightforward and I will try to add this once finding some spare time.

Borek, thanks for your previous replies. You wrote:

Regarding this implementation, we added the code reported here directly in out "bem" module, trying to integrate the contact pressure calculation.

We think that such modifications should be integrated directly in the beam2d/3d elements, and we are not sure about our implementation (which is only a trial, we currently obtain a constant soil pressure for a simple supported beam with subsoil, which is theoretically wrong).

Can you help on that?
many thanks in advance