The direct stiffness method

1 Overview and assumptions

The program solves linear static problems for skeletal structures built from prismatic members joined at nodes:

Assumptions (all of them are limitations too, see §15):

1. Linear elastic material, Hooke's law, no yielding or cracking.
2. Small displacements: equilibrium is written on the undeformed geometry; no $P$–$\Delta$ or $P$–$\delta$ effects, no geometric stiffness. Superposition of load cases is exact.
3. Euler–Bernoulli beam theory: plane sections remain plane and normal to the deformed axis; shear deformation and rotary inertia are neglected. The exact solution of $EI\,v'''' = q$ for a prismatic member is a cubic in $v$ between loads, so the Hermite cubic element is exact for any load that is applied with the same shape functions (§7), and no mesh refinement within a member is ever needed.
4. Prismatic members: $E, A, I_z, I_y, J$ are constant along each member. Warping torsion is not modelled (St. Venant torsion only).
5. Joints are rigid unless a member end is released (§8). Rigid-body offsets and rigid zones are not modelled.

2 Coordinate systems and sign conventions

Global axes. $X$ to the right, $Y$ up, $Z$ towards the viewer (right-handed). 2D models lie in the $X$–$Y$ plane; beam models lie on the $X$ axis. Gravity acts in $-Y$. Positive rotations follow the right-hand rule (in 2D: counter-clockwise, about $+Z$).

Member local axes. Local $x$ runs from node $i$ to node $j$ with unit vector $\mathbf{e}_x = (\mathbf{x}_j - \mathbf{x}_i)/L$.

• 2D: $\mathbf{e}_y$ is $\mathbf{e}_x$ rotated $+90^\circ$, $\mathbf{e}_y = (-e_{x,2},\,e_{x,1},\,0)$, $\mathbf{e}_z = +Z$. A member pointing right has local $y$ = global $Y$.

• 3D, non-vertical member ($|\mathbf{e}_x\cdot\hat Y| \le 0.999$): $\mathbf{e}_y$ is the projection of $\hat Y$ on the plane normal to $\mathbf{e}_x$,

  $$\mathbf{e}_y = \frac{\hat Y - (\mathbf{e}_x\cdot\hat Y)\,\mathbf{e}_x}{\lVert \hat Y - (\mathbf{e}_x\cdot\hat Y)\,\mathbf{e}_x\rVert},\qquad \mathbf{e}_z = \mathbf{e}_x\times\mathbf{e}_y .$$

• 3D, vertical member ($|\mathbf{e}_x\cdot\hat Y| > 0.999$): the projection above is undefined, so $\mathbf{e}_y = -\hat X$ when $\mathbf{e}_x$ points up ($+Y$) and $\mathbf{e}_y = +\hat X$ when it points down, and again $\mathbf{e}_z = \mathbf{e}_x\times\mathbf{e}_y$. For both a rising and a falling column this gives $\mathbf{e}_z = +\hat Z$.

• Roll angle $\beta$ (3D only): after the rule above, $\mathbf{e}_y$ and $\mathbf{e}_z$ are rotated about $\mathbf{e}_x$,

  $$\mathbf{e}_y' = \cos\beta\,\mathbf{e}_y + \sin\beta\,\mathbf{e}_z,\qquad \mathbf{e}_z' = -\sin\beta\,\mathbf{e}_y + \cos\beta\,\mathbf{e}_z .$$

  A roll of $90^\circ$ therefore makes local $y$ coincide with the old local $z$: bending in the vertical plane of a horizontal member is then resisted by $I_y$ instead of $I_z$.

Section properties. $I_z$ is the second moment about local $z$ (bending in the local $x$–$y$ plane — the "strong axis" for a normally oriented section), $I_y$ about local $y$, $J$ the torsion constant.

Member end forces (the stiffness-method convention) are the components of the force and moment acting on the member at its ends, along the member's local axes, positive in the positive local directions (moments by the right-hand rule). The end force vector of a member is denoted $\mathbf{q}'$.

Internal forces are functions of $x$ obtained from the free body $0 \to x$:

• $N$ axial force, tension positive;
• $M_z$ bending moment about local $z$, positive-face convention (for a horizontal 2D member with local $y$ up it is sagging-positive);
• $V_y$ shear in local $y$ defined by $\boxed{V_y = \mathrm{d}M_z/\mathrm{d}x}$ (so a simply supported beam under gravity has positive shear at its left support);
• $M_y$, $V_z = \mathrm{d}M_y/\mathrm{d}x$ and torsion $T$ (positive-face, right-hand rule about $+x$) in 3D.

Displacement functions along a member (local components): axial $d_x$, transverse $d_y, d_z$, twist $r_x$, $r_z = \mathrm{d}d_y/\mathrm{d}x$ and $r_y = -\mathrm{d}d_z/\mathrm{d}x$. Curvatures are $\mathrm{d}^2 d_y/\mathrm{d}x^2 = M_z/(EI_z) + \kappa_z$ and $\mathrm{d}^2 d_z/\mathrm{d}x^2 = -M_y/(EI_y) + \kappa_y$, where $\kappa$ are initial (thermal) curvatures (§7.5).

Thermal loads. $\Delta T$ is the uniform temperature change, $\Delta T_y = T(+y\text{ face}) - T(-y\text{ face})$ over the depth $h_y$, and $\Delta T_z$ likewise over the width $h_z$. The free thermal curvature is $\kappa_z = -\alpha\,\Delta T_y / h_y$ (top hotter $\Rightarrow$ hogging) and $\kappa_y = -\alpha\,\Delta T_z/h_z$.

Diagram drawing side (tension or compression side) is a display setting, never a data convention.

3 Degrees of freedom and numbering

Every node has up to six dofs in the canonical order $[u_x, u_y, u_z, \theta_x, \theta_y, \theta_z]$ (labels ux uy uz rx ry rz). An analysis type uses a subset:

Internally every member is built once as the full 12-dof 3D frame element, with local dof order $[u_1\ v_1\ w_1\ \theta_{x1}\ \theta_{y1}\ \theta_{z1} \mid u_2\ v_2\ w_2\ \theta_{x2}\ \theta_{y2}\ \theta_{z2}]$, and each analysis type simply selects the rows and columns of its own dofs (elemSub in dofs.js). One code path serves all five types.

Numbering: free dofs first. If $n_f$ of the $n$ global dofs are unknown, the dofs are numbered $1\dots n_f$ for the free dofs and $n_f+1\dots n$ for the restrained ones (supports, plus dofs restrained automatically, §9.3), nodes in model order and labels in canonical order. The global vectors are then naturally partitioned,

(The code uses 0-based indices; the report shows 1-based "code numbers".)

4 Element stiffness matrices

The stiffness matrix $\mathbf{k}'$ of a member relates its local end displacements $\mathbf{u}'$ to the local end forces it needs to hold them: $\mathbf{q}' = \mathbf{k}'\mathbf{u}'$ (for an unloaded member). Column $c$ of $\mathbf{k}'$ is the force vector that produces $u'_c = 1$ and all other $u'_{c'} = 0$.

4.1 Axial and torsional bars

Linear displacement field, $\varepsilon = \mathrm{d}u/\mathrm{d}x = (u_2-u_1)/L$, $N = EA\,\varepsilon$:

4.2 Beam element, 4×4

The Hermite cubic shape functions for $(v_1,\theta_1,v_2,\theta_2)$ at $s = x/L$ are

With $v(x) = \mathbf{N}\mathbf{u}'$, the strain energy $\tfrac12\int EI\,(v'')^2\,\mathrm{d}x$ gives $k_{ab} = \int_0^L EI\,N_a''N_b''\,\mathrm{d}x$:

4.3 2D truss element (global 4×4)

An axial-only bar has the local matrix $(EA/L)\,[1,-1;-1,1]$ in the dofs $(u_1',u_2')$. With $c = \cos\phi = (x_j-x_i)/L$, $s = \sin\phi = (y_j-y_i)/L$ the axial displacements are $u_1' = c\,u_i + s\,v_i$ etc., and $\mathbf{k} = \mathbf{T}^T\mathbf{k}'\mathbf{T}$ becomes

4.4 2D frame element (6×6)

Superposing (4.1) and (4.2) in the dof order $[u_1, v_1, \theta_1, u_2, v_2, \theta_2]$:

Rotation to global axes uses

(The rotation $\theta_z$ is the same in both systems because $z$ is common.) The beam element is the same matrix with the axial rows removed and the axis along $X$.

4.5 3D truss element (global 6×6)

With the direction cosines $\boldsymbol\lambda = (l, m, n) = (\mathbf{x}_j-\mathbf{x}_i)/L$:

The member force is $N = (EA/L)\,\boldsymbol\lambda\cdot(\mathbf{u}_j - \mathbf{u}_i)$, tension positive.

4.6 3D frame element (12×12)

In the dof order $[u_1\ v_1\ w_1\ \theta_{x1}\ \theta_{y1}\ \theta_{z1}\ u_2\ v_2\ w_2\ \theta_{x2}\ \theta_{y2}\ \theta_{z2}]$ the four deformation modes are uncoupled in local axes. Writing $a=12EI/L^3$, $b = 6EI/L^2$, $c = 4EI/L$, $e = 2EI/L$, the non-zero entries (upper triangle; the matrix is symmetric) are

• axial: $k_{1,1} = k_{7,7} = EA/L,\ k_{1,7} = -EA/L$;
• torsion: $k_{4,4} = k_{10,10} = GJ/L,\ k_{4,10} = -GJ/L$;
• bending about $z$ ($v,\theta_z$; $I_z$): $k_{2,2} = k_{8,8} = \tfrac{12EI_z}{L^3}$, $k_{2,6} = k_{2,12} = \tfrac{6EI_z}{L^2}$, $k_{2,8} = -\tfrac{12EI_z}{L^3}$, $k_{6,6} = k_{12,12} = \tfrac{4EI_z}{L}$, $k_{6,8} = k_{8,12} = -\tfrac{6EI_z}{L^2}$, $k_{6,12} = \tfrac{2EI_z}{L}$;
• bending about $y$ ($w,\theta_y$; $I_y$): the same pattern with $I_y$, but because $\theta_y = -\mathrm{d}w/\mathrm{d}x$ the $w$–$\theta$ coupling terms change sign: $k_{3,3} = k_{9,9} = \tfrac{12EI_y}{L^3}$, $k_{3,5} = k_{3,11} = -\tfrac{6EI_y}{L^2}$, $k_{3,9} = -\tfrac{12EI_y}{L^3}$, $k_{5,5} = k_{11,11} = \tfrac{4EI_y}{L}$, $k_{5,9} = k_{9,11} = +\tfrac{6EI_y}{L^2}$, $k_{5,11} = \tfrac{2EI_y}{L}$.

(indices are 1-based here.) The transformation is $\mathbf{T} = \mathrm{diag}(\boldsymbol\lambda,\boldsymbol\lambda,\boldsymbol\lambda,\boldsymbol\lambda)$ where the rows of the $3\times3$ matrix $\boldsymbol\lambda$ are the local unit vectors,

The axes $\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$ follow the local-axis rule of §2, including the vertical-member and roll rules.

5 Assembly

Let member $m$ join nodes $i$ and $j$ and let its code number vector $\mathbf{c}_m = [\,g(i,1),\dots,g(i,n_d),\ g(j,1),\dots,g(j,n_d)\,]$ list the global dof index of each element dof (this is the "code" column of the report). The global stiffness matrix and the assembled load vectors are

(spring stiffnesses are added to the diagonal $K_{gg}$ of the dof they act on, §9.1). Because the free dofs are numbered first, $\mathbf{K}_{ff}$, $\mathbf{K}_{fs}$, $\mathbf{K}_{sf}$, $\mathbf{K}_{ss}$ are simply the leading/trailing blocks of $\mathbf{K}$.

The solver stores $\mathbf{K}_{ff}$ in skyline (profile) form after a reverse Cuthill–McKee reordering, and factorises it once as $\mathbf{L}\mathbf{D}\mathbf{L}^T$. Every load case and every combination is then one additional right-hand side, so the stiffness matrix is never re-factored.

6 Partitioned equations, displacements and reactions

Equilibrium of the whole structure, with $\mathbf{P}$ the joint loads and $\mathbf{FEF}$ the assembled fixed-end forces of the member loads (§7), reads $\mathbf{K}\mathbf{U} + \mathbf{FEF} = \mathbf{P} + \mathbf{R}$ with $\mathbf{R}$ the support reactions (non-zero only on restrained dofs). Partitioning:

Hence

$\mathbf{U}_s$ holds the prescribed support displacements (zero unless a settlement is specified, §10). Reactions are the forces the supports exert on the structure. For a restrained node with an elastic spring on a free dof the "reaction" is $-k_s u$.

A structure is stable when $\mathbf{K}_{ff}$ is positive definite. A zero or negative pivot during factorisation signals a mechanism (§13).

7 Fixed-end forces of member loads

7.1 Work-equivalent loads

Consider a member whose two ends are fully fixed ($\mathbf{u}' = \mathbf{0}$) under a distributed transverse load $q(x)$ (positive in $+y$). The fixed-end force vector $\mathbf{q}_F$ is the set of end forces needed to hold the ends — they are the member-end reactions. By the principle of virtual work with the same Hermite displacement field as the stiffness matrix,

for distributed loads, point forces $P_k$ at $x = a_k$ (positive $+y$) and concentrated couples $M_k$ (positive counter-clockwise). The sign is negative because $\mathbf q_F$ are the forces on the member by the fixed ends: a downward load ($q<0$) gives upward end forces. In the structure the member loads are reversed into joint loads: a member with end forces $\mathbf q_F$ in the fixed state is equivalent to the loading $-\mathbf q_F$ on the joints, which is why (6.1) contains $-\mathbf{FEF}_f$.

Because (4.1) are the exact solution of the homogeneous beam equation, (7.1) is exact for Euler–Bernoulli members, whatever the load. The code integrates distributed loads with a 3-point Gauss rule (exact for the quartic integrand "cubic × linear"), evaluates $N_a(a)$ for point loads and $N_a'(a)$ for couples, and handles axial and torsional loads with linear shape functions.

7.2 Closed forms (downward loads, magnitudes $P$, $w$; $a$ from node $i$, $b = L-a$)

Writing $\mathbf q_F = [V_i,\ M_i,\ V_j,\ M_j]$ (beam dofs; $V$ up, $M$ counter-clockwise):

Derivation of the point-load row. $\mathbf N(a)$ from (4.1) with $s=a/L$: $N_1 = 1-3s^2+2s^3 = b^2(3a+b)/L^3$, $N_2 = L s(1-s)^2 = ab^2/L^2$, $N_3 = s^2(3-2s) = a^2(a+3b)/L^3$, $N_4 = -L s^2(1-s) = -a^2b/L^2$. Then $\mathbf q_F = -\mathbf N^T(a)\,(-P) = P\,\mathbf N^T(a)$ — the classic fixed-end values. The couple row follows from $N_1' = -6ab/L^3$, $N_2' = b(b-2a)/L^2$, $N_3' = 6ab/L^3$, $N_4' = a(a-2b)/L^2$ and $\mathbf q_F = -M\,\mathbf N'^T(a)$. For a UDL, $\int_0^L N_1\,\mathrm{d}x = L/2$, $\int_0^L N_2\,\mathrm{d}x = L^2/12$, giving the table entries. Equilibrium of the member is a useful check on every row, e.g. for the point load: $V_i + V_j = P$ and $M_i + M_j + V_jL - Pa = 0$.

Loads that are not applied along local axes are first resolved into local components $q_x, q_y, q_z$ (global loads with $\hat Y$ etc.; "projected" loads are per unit projected length, so the resolved intensity is $w\cdot L_x/L$ for a load in $Y$ and a member with horizontal extent $L_x$). Partial and trapezoidal loads are integrated over their actual extent $[x_1, x_2]$.

7.3 Moment loads

A concentrated couple $M$ about local $z$ at $x=a$ gives the last row of the table. In 3D, a couple about local $y$ couples to $(w,\theta_y)$ with the sign of $\theta_y = -\mathrm dw/\mathrm dx$ taken into account; a couple about $x$ is carried by the torsion dofs with linear shape functions.

7.4 Axial loads, self-weight

An axial distributed load $q_x$ gives $\mathbf q_F = -\int \mathbf N_a^T q_x\,\mathrm dx$ with the linear functions $N = [1-s,\ s]$; a uniform $q_x$ gives $-q_xL/2$ at each end. Self-weight is the load $w = \rho\,g\,A$ per unit length acting in global $-Y$ over the full member length; its total is $\rho g A L$ whatever the inclination.

7.5 Temperature changes and lack of fit

Initial strains are treated like loads: with both ends fixed, the stress-producing part of the strain is the negative of the free strain. A member with a free thermal strain $\varepsilon_0 = \alpha\Delta T$ and a free thermal curvature $\kappa_z = -\alpha\Delta T_y/h_y$, held fixed, develops

(compression for a temperature rise; a sagging moment when the $+y$ face is the hotter one). In terms of the stored fixed-end forces, $N(0) = -q_{u1}$ and $M_z(0) = -q_{\theta z1}$ give

A lack of fit $\Delta$ (member fabricated too long, $\Delta>0$) is a free strain $\varepsilon_0 = \Delta/L$, so a member forced between fixed points carries

In an unrestrained member a thermal strain produces displacement but no force; the displacement functions of §12 include $\varepsilon_0$ and $\kappa$ so that deflections are exact.

8 Member end releases (static condensation)

Let the dofs of an element be split into retained ($r$) and released ($c$, "condensed") ones. A released dof has zero end force, so the corresponding row of $\mathbf q = \mathbf k\mathbf u + \mathbf q_F$ reads $0 = \mathbf k_{cr}\mathbf u_r + \mathbf k_{cc}\mathbf u_c + \mathbf f_c$ and

Substituting back into the retained rows gives the condensed stiffness and fixed-end forces,

The condensed matrix keeps its size; rows/columns of released dofs become zero. The code performs one dof at a time (Gaussian elimination), which handles any combination of releases; $\mathbf u_c$ is recovered with (8.1) for the diagrams. Example — hinge at node $j$ of a beam ($c = \theta_2$):

A truss (axial-only) member is the same operation with all bending and torsion dofs released: its matrix is $(EA/L)[1,-1;-1,1]$, and transverse member loads are condensed to simply-supported end reactions applied at its end nodes (a stiffness of arbitrary positive value is used for the condensation, which does not influence the result). A torsion release at both ends has no stiffness at all; such a dof is simply left without stiffness.

9 Springs and inclined supports

9.1 Elastic springs

A spring of stiffness $k_s$ acting on a free dof $g$ adds $k_s$ to $K_{gg}$; it is an extra stiffness element grounded at the node. Its reaction is $-k_s U_g$. Springs may carry any of the six dofs (N/m or N·m/rad).

9.2 Inclined supports and the nodal rotation $\boldsymbol\Gamma$

A roller on an inclined plane restrains only the displacement normal to the plane. Such a node uses a support frame rotated by $\vartheta$ (angle) about $Z$; its dofs are $\mathbf u_s = \boldsymbol\Gamma\,\mathbf u$ with

The support-frame $x$ axis is $(\cos\vartheta,\sin\vartheta)$, the $y$ axis $(-\sin\vartheta,\cos\vartheta)$. Restraints are applied in the support frame: "$u_y$ restrained" gives a reaction along the support-frame $y$ axis, i.e. normal to the plane. Every member end at such a node uses the modified transformation $\tilde{\mathbf T} = \mathbf T\,\boldsymbol\Gamma^T$ so that $\mathbf u' = \mathbf T\,\boldsymbol\Gamma^T\mathbf u_s$ and $\mathbf k = \tilde{\mathbf T}^T\mathbf k'\tilde{\mathbf T}$. The solved $\mathbf U_f$, $\mathbf R$ are in the support frame at those nodes and are rotated back ($\mathbf u = \boldsymbol\Gamma^T\mathbf u_s$) for reporting; reactionsLocal keeps the support-frame components. Inclined supports exist only in the 2D truss and frame analyses.

9.3 Auto-restrained dofs

A rotational dof that has no stiffness — e.g. a node where only axial-only members or fully hinged members meet — would make $\mathbf K_{ff}$ singular although it has no physical meaning. Such dofs (and all dofs of nodes not connected to any member) are restrained automatically and reported as information (auto-restrained). Translational mechanisms are never hidden this way (§13).

10 Prescribed displacements (settlement)

A support displacement $\bar{\mathbf U}_s\neq\mathbf 0$ (settlement, imposed rotation) enters through $\mathbf U_s$ in (6.1) and (6.2). A settlement $\Delta$ of one support of a determinate structure causes rigid-body motion and no member forces, while in an indeterminate structure it induces forces: for a fixed–fixed beam whose right end sinks by $\Delta$, $q_{v1} = 12EI\Delta/L^3$, $q_{\theta 1} = q_{\theta 2} = 6EI\Delta/L^2$, i.e. moments $6EI\Delta/L^2$ at each end of opposite sense (hogging at the left end, sagging at the right end).

11 Member end forces

After solving, each member's local end displacements are $\mathbf u' = \mathbf T\,\tilde{\mathbf u}$ (the member's rows of $\mathbf U$, in the support frame where applicable) and

with $\mathbf k'$ and $\mathbf q_F$ already condensed for any releases. Released end displacements are recovered with (8.1). The member end forces in global axes are $\mathbf q = \mathbf T^T\mathbf q'$.

12 Exact internal-force and deflection functions

The end forces and the member loads are all that is needed to write equilibrium of the free body from $0$ to $x$ (with all loads between $0$ and $x$ included):

with the sums over loads at positions $\le x$. By construction $V_y = \mathrm dM_z/\mathrm dx$ and $V_z = \mathrm dM_y/\mathrm dx$. Within each interval between breakpoints (member ends, point loads and couples, start/end of distributed loads) every function is an exact polynomial in $s = x - x_0$: $N$ and $V$ are one degree above the load, $M$ two. Point forces and couples produce jumps in $V$ and $M$ that are applied at the breakpoints; a distributed trapezoid gives a quadratic $V$ and a cubic $M$. The segments are stored as coefficient arrays, so extremes are found by root-finding of the derivative, not by sampling.

Deflections follow by integrating the curvature and strain, starting from the member's own end displacements $\mathbf u'$ at $x=0$:

For a simply supported beam under a UDL this yields $d_y(L/2) = -5wL^4/(384EI)$; for a cantilever with tip load $P$ (downward) $d_y(L) = -PL^3/(3EI)$, in agreement with the closed forms. Points of interest — the extreme values of every function, the zeros of $V$ (location of the extreme moment) and the largest deflection — are found as the roots of the derivative polynomial in each segment, plus the breakpoints themselves.

For a truss member carrying transverse loads, the diagrams are those of a simply supported member under the same load (end reactions applied to the joints).

13 Self-checks

Three checks always run; failures become warnings or errors.

1. Global equilibrium. The resultant of the applied loads (joint loads plus the resultants of all member loads, including self-weight) and the reactions, about the global origin, $\sum\mathbf F_\text{applied} + \sum\mathbf R \approx \mathbf 0$ and $\sum\mathbf M_\text{applied} + \sum\mathbf M_R \approx \mathbf 0$, with a tolerance relative to the magnitude of the loads ($10^{-6}$).
2. End compatibility. The internal forces of §12 evaluated at $x=L$ must equal the $j$-end forces: $N(L) = q_{u2}$, $V_y(L) = -q_{v2}$, $M_z(L) = +q_{\theta z2}$, and the integrated deflection at $x=L$ must equal the $j$-end displacement. This ties the displacement functions to the end forces and verifies the sign conventions end to end.
3. Instability. During factorisation of $\mathbf K_{ff}$ a zero or negative pivot (relative tolerance $10^{-11}$) means that a dof has no stiffness. The solver reports "the structure is unstable (a mechanism)", naming the node and the dof (and flags inclined support frames), instead of returning garbage.

Further model checks (unknown nodes or sections, zero-length members, dofs a load cannot act on) are reported as errors or warnings by the model validation.

14 Load combinations and envelopes

The structure is linear, so the response to a combination $\sum_k \gamma_k\,\mathrm{LC}_k$ is the factored sum of the case responses. It is computed by applying the same factors to the joint loads, fixed-end forces and prescribed displacements, and solving one more right-hand side with the existing factorisation; displacements, reactions, end forces and the piecewise-polynomial diagrams therefore equal $\sum_k\gamma_k$(case result), to round-off.

The envelope is computed over the combinations (or over the load cases when no combination is defined): at every station $x$ of every member the maximum and minimum of $N, V, T, M, d$ over all combinations, and for every support the maximum and minimum reaction component. Hence at every station the envelope maximum is at least, and the minimum at most, the value of every combination.

15 Limitations and references

Not modelled: geometric non-linearity and second-order effects, buckling, shear deformation (Timoshenko beams), warping torsion and thin-walled sections, tapered members, rigid end zones and joint offsets, non-linear or tension-only/compression-only elements, material non-linearity, dynamic and seismic response, and soil–structure interaction other than linear springs. Results are only as good as the idealisation: the user remains responsible for code checks and design.

References

1. A. Kassimali, Matrix Analysis of Structures, 2nd ed., Cengage Learning. (Notation for the stiffness method, member code numbers, partitioning, fixed-end forces, releases, inclined supports.)
2. W. McGuire, R. H. Gallagher, R. D. Ziemian, Matrix Structural Analysis, 2nd ed., Wiley. (3D frames, local axes and roll angle, static condensation, work-equivalent loads.)
3. R. C. Hibbeler, Structural Analysis, 10th ed., Pearson. (Closed-form fixed-end moments and deflection formulas used as benchmarks.)
4. J. S. Przemieniecki, Theory of Matrix Structural Analysis, McGraw-Hill. (Element matrices.)

The numerical results of the JavaScript solver are verified in tests/solver.test.js (closed-form benchmarks, derived by hand in the comments) and tests/crosscheck.test.js (agreement with the independent Python implementation python/stiffstruct.py on 49 models).