Model produces results when unstable with a hinge

[weikequ]

Describe the bug
Define a simply supported beam with a hinge in the middle and a distributed load. Analyzing the structure using analyze() and it will run propery, but give erroneous reactions, moments, shears, and displacements. Additionally, running using analyze(sparse=False) will give a different set of results.

See video below.

To Reproduce
Here's the code to reproduce the behaviour:

from Pynite import FEModel3D

model = FEModel3D()

Iy = 0.000015234  # m^4 weak axis
Iz = 0.0000712  # m^4 strong axis
J = 0.000000241  # m^4
A = 0.006264504  # m^2
model.add_section("W10x33", A, Iy, Iz, J)

E = 200000000  # kPa
G = 77000000  # kPa
nu = 0.3  # Poisson's ratio
rho = 77  # Density (kN/m^3)
model.add_material("steel", E, G, nu, rho)

model.add_node("n1", 0, 0, 0)
model.add_node("n2", 1, 0, 0)
model.add_node("n3", 2, 0, 0)

model.add_member("m1", "n1", "n2", "steel", "W10x33")
model.add_member("m2", "n2", "n3", "steel", "W10x33")

model.def_support("n1", 1, 1, 1, 1, 1, 0)
model.def_support("n3", 0, 1, 1, 1, 1, 0)

model.add_member_dist_load("m1", "Fy", 10, 10)
model.add_member_dist_load("m2", "Fy", 10, 10)

model.def_releases("m2", Rzi=True)

model.analyze(sparse=False)

print(model.nodes["n1"].RxnFY)
print(model.nodes["n3"].RxnFY)
print(model.nodes["n2"].DY)

Expected behavior
An error should be thrown about the model being unstable.

Screenshots

Screen.Recording.2025-07-09.at.12.33.39.PM.mov

[angelmusonda]

Hi @weikequ,

You've raised an important issue, which you also brought up in #255. @JWock82 mentioned that he would take a closer look.

I'd like to offer a few comments.

Firstly, Pynite internally uses the "is" operator for Boolean comparisons instead of "==". This means that using 0 and 1 in methods that expect True and False may not behave as expected. The "is" operator checks for object identity, not value equality. As a result, in your model, all degrees of freedom at nodes n1 and n3 were fully restrained. However, this is beside the point, because even if you write:

model.def_support("n1", True, True, True, True, True, False)
model.def_support("n3", False, True, True, True, True, False)

you will still not get an instability error, for reasons explained below.

There are essentially two types of instability to consider: nodal and global.

Pynite checks for nodal instability by examining each node to verify that all six degrees of freedom are adequately restrained. In your case, this check passes — each node has sufficient constraints.

The issue lies with global instability. Pynite uses numpy.linalg.solve() for dense systems (when the sparse option is off) and scipy.sparse.linalg.spsolve() for sparse systems. It relies on these solvers to raise an error when a singular matrix is encountered, which would indicate a globally unstable structure.

However, this approach is unreliable:

scipy.sparse.linalg.spsolve() does not raise an error for singular matrices as far as I can see in the documentation.

numpy.linalg.solve() may fail to detect singularity due to floating-point precision issues.

A more robust approach is to explicitly check for singularity before solving the system. Computing the matrix rank would be one of the reliable methods to achieve this.

For the dense case, numpy.linalg.matrix_rank() can be used — it accepts a tolerance parameter and works well.

For the sparse case, one can perform an LU factorization and count the number of diagonal entries in the upper triangular matrix that exceed a given tolerance. This count effectively gives the rank of the matrix.

The computed rank can then be compared to the total number of degrees of freedom. If the rank is lower than the number of degrees of freedom, then the structure is unstable.