I looked into this and the current solution is to use two separate nodes for the zero length elements rather than using the same node. I’ll probably work on getting your version of this example to work in a future minor release. Here’s how your example would change.
using GXBeam, StaticArrays, LinearAlgebra
stiff = [1000000 0 0 -.5 -1 -50000; 0 3000000 0 0 0 0; 0 0 3000000 0 0 0; 0 0 0 7 .1 -.02; 0 0 0 0 5 .1; 0 0 0 0 0 3000]
stiff = Symmetric(stiff)
mass = [.02 0 0 0 -5e-7 -1e-7; 0 .02 0 5e-7 0 1e-4; 0 0 .02 1e-7 -1e-4 0; 0 0 0 1e-5 1e-8 2e-10; 0 0 0 0 6e-7 9e-9; 0 0 0 0 0 1e-5]
mass = Symmetric(mass)
nodes = [SA_F64[0,0,0], SA_F64[0,.25,0], SA_F64[0,.25,0], SA_F64[0,.5,0], SA_F64[0,.5,0]]
start = 1:4
stop = 2:5
stiffness = [stiff,zeros(6,6),stiff,zeros(6,6)]
pointmass=Dict(2=>PointMass(mass),4=>PointMass(mass))
transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:4];
assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10, point_masses=pointmass);
I’d be curious to know how your lumped mass workarounds compare to what I’ve implemented. I don’t have a proper test for the lumped masses, but since they effectively just modify the beam element mass matrix there isn’t much to mess up with their implementation.
With regard to a citation, I don’t currently have a reference that goes with the code so I’ll get back to you on that.