*BUCKLE does not check its result: with a singular stiffness matrix it prints factors near the shift 1.0 that are not eigenvalues (exit 0, no message)

Analysis issues
1

Summary

Section 10.7.4 (Buckling calculations) of the 2.23 manual says the ARPACK buckling mode solves K U = lambda KG U “where [K] is symmetric and positive definite”. Golden rule 2 (section 4) asks the user to check that there are “no unrestrained rigid body modes”. So a singular K is a modelling error, and this report does not ask ccx to buckle such a model. It asks for a check.

What ccx does now with a singular K: it prints buckling factors clustered around the shift sigma = 1. They are not eigenvalues of the matrices ccx itself stores. The first factor of one and the same deck changes with the number of factors requested. Exit status 0, no warning, spooles.out empty. A tighter accuracy (second field of the *BUCKLE data line 0.0001 or 1e-8) gives the same values.

Section 10.7.4 also describes a plausibility loop: accept the first factor only if 5 sigma < factor < 50000 sigma, else set sigma = factor/500 and repeat, up to four times. It would flag these values. In arpackbu.c that loop is commented out, and nothing else checks the result.

A singular K is not always the user’s visible mistake. stub.inp is an S4 cantilever plate (clamped at x = 0, pushed at the tip by two forces of 1) with a short B31 stub in the plate’s plane at tip node 3. The manual (section 6.2.14) says beam and shell elements sharing a node “are always connected in a stiff way”, and ccx makes a knot in node 3 (stub.12d: “a KNOT was generated in node 3”). But the knot’s shell nodes lie on the normal through node 3, so the shell gives the knot no stiffness about the normal. The stub’s other end is free, so knot and stub can turn about z without straining anything. Stored with SOLVER=MATRIXSTORAGE, K has a smallest eigenvalue of -6.0e-9 and a largest of 7.5e8.

factors requested factors printed
1 4.987410
2 0.9065606, 1.104038
3 0.9945662, 1.005215, 1.030989
4 0.9970473, 0.9981509, 1.001781, 1.002897

None of the printed values below 4.987 is an eigenvalue of the stored K and KG. The only eigenvalue below 4.987 is 2e-6, from the free turn. An NLGEOM static run of the same model at 4.5 times the load is stable and deflects like the plate without the stub. Whether the stub’s node is loaded does not matter: with the stub on the unloaded node 5 (lines 7, 10., 11., 0. and 3, 5, 7), two factors come out as 1.868809 and 4.988736.

The same happens with a plain mechanism of the user’s own making (twist.inp below, a column whose base may twist): 22.30127 for one factor, 0.9208384 and 1.099216 for two, against 105.1482 for the clamped twin.

Version

ccx 2.23 (macOS arm64). Source: master 9ba5d76 (16 Oct 2026).

Already reported or fixed?

Not reported: I found no GitHub issue (open or closed) or pull request in Dhondtguido/CalculiX, and no calculix.discourse.group topic, that traces buckling factors near 1 to a singular stiffness matrix or asks for a check (searched 21 Oct 2026). Not fixed: arpackbu.c on master (9ba5d76) still sets sigma = 1 with the loop commented out.

Related:

New here: a small inline deck that shows where values near 1 come from (a singular K), that they are not eigenvalues of the stored matrices, and that the accuracy setting does not change them; and a request for a check.

Minimal input

stub.inp:

*NODE, NSET=NALL
1, 0., 0., 0.
2, 10., 0., 0.
3, 20., 0., 0.
4, 0., 10., 0.
5, 10., 10., 0.
6, 20., 10., 0.
7, 21., 0., 0.
*ELEMENT, TYPE=S4, ELSET=PLATE
1, 1, 2, 5, 4
2, 2, 3, 6, 5
*ELEMENT, TYPE=B31, ELSET=STUB
3, 3, 7
*MATERIAL, NAME=STEEL
*ELASTIC
210000., 0.3
*SHELL SECTION, ELSET=PLATE, MATERIAL=STEEL
0.2
*BEAM SECTION, ELSET=STUB, MATERIAL=STEEL, SECTION=RECT
0.2, 0.2
0., 0., 1.
*BOUNDARY
1, 1, 6
4, 1, 6
*STEP
*BUCKLE
2
*CLOAD
3, 1, -1.
6, 1, -1.
*END STEP
twist.inp (B31 column, base free to twist: a plain user error)

twist.inp:

*NODE, NSET=NALL
1, 0., 0., 0.
2, 25., 0., 0.
3, 50., 0., 0.
4, 75., 0., 0.
5, 100., 0., 0.
*ELEMENT, TYPE=B31, ELSET=COLUMN
1, 1, 2
2, 2, 3
3, 3, 4
4, 4, 5
*MATERIAL, NAME=STEEL
*ELASTIC
210000., 0.3
*BEAM SECTION, ELSET=COLUMN, MATERIAL=STEEL, SECTION=RECT
2., 3.
0., 0., 1.
*BOUNDARY
1, 1, 3
1, 5, 6
*STEP
*BUCKLE
2
*CLOAD
5, 1, -1.
*END STEP

The clamped twin replaces the two *BOUNDARY lines by 1, 1, 6.

NLGEOM check of stub.inp (replaces the *BUCKLE step) 
-*STEP
-*BUCKLE
-2
-*CLOAD
-3, 1, -1.
-6, 1, -1.
-*END STEP
+*STEP, NLGEOM, INC=100
+*STATIC
+0.05, 1., 0.001, 0.05
+*CLOAD
+3, 1, -4.5
+6, 1, -4.5
+3, 3, 0.001
+6, 3, 0.001
+*NODE PRINT, NSET=NALL
+U
+*END STEP

How to run

ccx stub and ccx twist. For the other rows, change the data line under *BUCKLE (2) to 1, 3 or 4; for the accuracy test, to 2, 0.0001 or 3, 1.e-8. The plate without the stub: delete node 7, element 3 with its *ELEMENT line, and the three *BEAM SECTION lines.

Dense check: change the step to *BUCKLE, SOLVER=MATRIXSTORAGE and solve the generalised problem K x = lambda KG x from stub.sti and stub.str (e.g. scipy.linalg.eigvals(K, KG) on the symmetric matrices built from the listed upper triangles).

Expected

Section 10.7.4 assumes K positive definite, so for a singular K the expected outcome is an *ERROR (or at least a *WARNING) that the stiffness matrix is singular or not positive definite, not a list of factors.

If factors are printed, they should be eigenvalues of K x = lambda KG x. For stub.inp the stored matrices have one eigenvalue of 2e-6 (the free turn: the null vector of K is the knot’s rotation about z, with the stub’s far nodes moving in y), then 4.98741 and 60.8287, and nothing in between. Other references for the plate:

  • *BUCKLE with one factor: 4.987410.
  • The plate without the stub: 4.988630 and 60.82335.
  • The knot’s rotation about the normal held (3, 6, 6 under *BOUNDARY): 4.98741 and 60.82871.
  • The NLGEOM run (a push of 0.001 in z at nodes 3 and 6): tip deflection uz(3) = 3.301987E-02 at 4.5 times the load, against 3.301897E-02 for the plate without the stub; at 1.5 times the load 4.644720E-03 against 4.644712E-03. Nothing buckles near a factor of 1.

For twist.inp: the same message (the base twist is free); the clamped twin gives 105.1482 and 236.3637.

Actual

Exit 0, no warning, spooles.out empty. For stub.inp, the table in the Summary; stdout shows sigma=1.000000,d[0]=0.906561 for two factors. stub.dat:

 MODE NO       BUCKLING
                FACTOR

      1   0.9065606E+00
      2   0.1104038E+01

With 2, 0.0001 and 2, 1.e-8 the two factors are the same. With three factors and accuracy 0.0001 or 1e-8: 0.9958893, 1.003983, 1.069215.

twist.inp:

factors requested factors printed
1 22.30127
2 0.9208384, 1.099216
3 0.7257664, 1.512910, 1.573933
4 0.8967395, 0.9064208, 1.115846, 1.140050

Where it goes wrong

arpackbu.c (master 9ba5d76; the 2.23 source has the same lines):

  • line 78: which[3]="LM"; line 588: iparam[6]=4; (ARPACK buckling mode, OP = (K - sigma KG)^-1 K, with K as the inner product); line 508: sigma=1.;.
  • lines 744-754: the plausibility check of section 10.7.4 is commented out:
    /* if((5.>d[0]/sigma)||(50000.<d[0]/sigma)){ … else{iconverged=0;}*/

ARPACK’s buckling mode uses K as its inner product, which needs K positive definite. With a singular K the values it returns are not eigenvalues of K x = lambda KG x (dense check above). They cluster around sigma = 1, and nothing checks them.

How stub.inp gets a singular K: the knot MPCs (fillknotmpc.f, lines 198-242) give each knot node the displacement of node 3 plus the cross product of the knot rotation with the node’s offset from node 3. For the shell’s own expanded nodes that offset lies along the normal, so the rotation about the normal drops out of their equations; only the stub’s nodes carry it, and nothing holds the stub’s far end. (The knot is labelled KNOT1, checked in a debugger: its node list also holds four nodes of the B31 expansion that no element uses and that are still at (0,0,0) at that point, so the node cloud is 3-D and line 158 picks the line formula. The plane formula, lines 243-314, has the same rotation terms.)

Workaround / fix

Fix: check the result before printing it. For example, compute the residual of K x - lambda KG x for each returned pair, and stop with an *ERROR that names a singular or not positive definite stiffness matrix when it is large. The documented loop (re-enabled by pull request 44) would reject a first factor below 5 sigma. I could not build ccx, so whether the loop alone then finds 4.987 or stops was not tested.

Workaround: remove the mechanism. For stub.inp, hold the knot’s rotation about the normal (3, 6, 6), hold the stub’s far end, or leave the stub out. For twist.inp, hold DOF 4. A quick test: run *BUCKLE once with one factor and once with two; if the first factor changes, the model probably has a zero-stiffness mode.