matSpringDampMass

matSpringDampMass

Builds the system matrices for the granular eigenvalue problem M*x''(t) + Gamma*x'(t) + K*x(t) = 0: the stiffness matrix (Hessian) matSpring, the damping matrix matDamp, and the mass matrix matMass for a 2D or 3D packing with Hookean contacts. Degrees of freedom are ordered [x1, y1, (z1), x2, y2, (z2), ...], so the matrices are (dof x dof) with dof = N*dim.

Default Syntax

packing = data.packing;     % struct, loaded from a pack.m results .mat file
[K, Gamma, M] = matSpringDampMass(packing, dampingConstant, springConstant, mass)

With dampingConstant = 1, springConstant = 1, mass = 1 (all defaults). The packing struct is detected automatically: 2D vs 3D from the presence of vecPosZ, new format (vecPosX/scalBoxWidthX/...) vs old format via the oldPackingFormat option.

Options

opts.algorithm        = "cell";     % "cell" (default, linear scaling, fast for large N)
                                    % "particle" (all pairs, quadratic scaling, fine for small N)
opts.periodic         = false;      % logical, periodic boundary conditions
opts.sparseOutput     = false;      % logical, return sparse matrices
opts.uniformMass      = false;      % logical, all particles get the same mass
opts.oldPackingFormat = false;      % logical, packing uses x/y(/z), Dn, Lx/Ly(/Lz) fields
[K, Gamma, M] = matSpringDampMass(packing, dampingConstant, springConstant, mass, opts)

• algorithm = "cell" buckets particles into cells of width 4*max(radii) and only checks the neighboring cells — use for large packings. "particle" checks every pair — use for small ones.
• periodic must be true for the cell algorithm to run (the non-periodic cell branch is not implemented). For non-periodic 2D packings, wall contacts (particle within its radius of a wall) add K to the diagonal.
• sparseOutput returns sparse K/Gamma/M — needed before passing to polyeig for large systems.
• uniformMass gives every particle mass mass; otherwise mass = mass * (unit-ball volume) * r^dim.

Examples

data = load('data/eigenData/results_2D_iso_N100_P0.1_Seed5_gamma_3.00e-03.mat');

% Damped system matrices, dense
[K, Gamma, M] = matSpringDampMass(data.packing, 0.01, 100, 1);

% Sparse, periodic — ready for polyeig
opts.periodic      = true;
opts.sparseOutput  = true;
[K, Gamma, M] = matSpringDampMass(data.packing, 0.001, 100, 1, opts);

[eigenVectors, eigenValues] = polyeig(K, Gamma, M);

% Old-format packing struct (x/y/Dn/Lx/Ly)
packing = struct('x', x, 'y', y, 'Dn', Dn, 'Lx', Lx, 'Ly', Ly);
opts.oldPackingFormat = true;
opts.periodic         = true;
[K, Gamma, M] = matSpringDampMass(packing, 1, 100, 1, opts);

% Small packing, all-pairs check
opts2 = struct('algorithm', "particle");
[K, Gamma, M] = matSpringDampMass(packing, 1, 100, 1, opts2);