From 9eeae870bdb672841d32f139b808236d86c8c560 Mon Sep 17 00:00:00 2001 From: Ayberk Korkmaz Date: Wed, 16 Sep 2026 09:04:05 +0300 Subject: [PATCH] Fix node 2 rotational spring corrupting the constitutive matrix constitutive_matrix() reduces matrix[1][1] correctly for a node 1 spring, but for node 2 it was touching matrix[2][1] and matrix[1][2] (the moment/rotation coupling terms) instead of matrix[2][2]. A spring at node 2 never actually softened that node's own stiffness, and the off-diagonal terms ended up wrong and asymmetric, which a linear-elastic element's constitutive matrix should never be. The function's own comment above this block already states the correct form (c11 + 1/k1 on one diagonal, c22 + 1/k2 on the other, off-diagonal unchanged), the code just didn't match it for the node 2 case. Fix: matrix[2][2] gets the node 2 spring reduction, off-diagonal terms are left alone for both nodes. Added test_constitutive_matrix.py with a symmetry check, a check that each spring only touches its own diagonal term, and a check for the node-2-only case specifically. Verified all three fail on the old code with the exact wrong values, pass with the fix. Full suite passes, 214 tests (matplotlib needed to be installed for the plotter tests to run at all, otherwise they silently fall back to a null plotter and fail regardless of this change). --- anastruct/fem/elements.py | 4 +- tests/test_constitutive_matrix.py | 71 +++++++++++++++++++++++++++++++ 2 files changed, 72 insertions(+), 3 deletions(-) create mode 100644 tests/test_constitutive_matrix.py diff --git a/anastruct/fem/elements.py b/anastruct/fem/elements.py index b204c41..a97fece 100644 --- a/anastruct/fem/elements.py +++ b/anastruct/fem/elements.py @@ -365,10 +365,8 @@ def constitutive_matrix( # [ c21 , c21 + 1/k2 ]] if 1 in spring and spring[1] != 0: # not hinge matrix[1][1] = 1 / (1 / matrix[1][1] + 1 / spring[1]) - matrix[2][1] = 1 / (1 / matrix[2][1] + 1 / spring[1]) if 2 in spring and spring[2] != 0: # not hinge - matrix[2][1] = 1 / (1 / matrix[2][1] + 1 / spring[2]) - matrix[1][2] = 1 / (1 / matrix[1][2] + 1 / spring[2]) + matrix[2][2] = 1 / (1 / matrix[2][2] + 1 / spring[2]) return matrix diff --git a/tests/test_constitutive_matrix.py b/tests/test_constitutive_matrix.py new file mode 100644 index 0000000..4cbea3b --- /dev/null +++ b/tests/test_constitutive_matrix.py @@ -0,0 +1,71 @@ +import unittest + +import numpy as np + +from anastruct.fem.elements import constitutive_matrix + + +class TestConstitutiveMatrixSprings(unittest.TestCase): + def test_symmetric_with_springs_on_both_nodes(self): + # A linear-elastic element's constitutive matrix must be symmetric, + # spring release or not. + matrix = constitutive_matrix( + EA=1e9, + EI=15000.0, + l=4.0, + spring={1: 1000.0, 2: 2000.0}, + node_1_hinge=False, + node_2_hinge=False, + ) + np.testing.assert_allclose(matrix, matrix.T) + + def test_spring_reduces_only_its_own_node_diagonal_term(self): + # A rotational spring in series with the beam's own rotational + # stiffness should only soften that node's diagonal term (per the + # flexibility matrix in the function's own comment); the + # off-diagonal (moment/rotation coupling) term is unaffected. + EI, l = 15000.0, 4.0 + k1, k2 = 1000.0, 2000.0 + + matrix = constitutive_matrix( + EA=1e9, + EI=EI, + l=l, + spring={1: k1, 2: k2}, + node_1_hinge=False, + node_2_hinge=False, + ) + + c11_expected = 1 / (1 / (4 * EI / l) + 1 / k1) + c22_expected = 1 / (1 / (4 * EI / l) + 1 / k2) + + np.testing.assert_allclose(matrix[1][1], c11_expected) + np.testing.assert_allclose(matrix[2][2], c22_expected) + np.testing.assert_allclose(matrix[1][2], -2 * EI / l) + np.testing.assert_allclose(matrix[2][1], -2 * EI / l) + + def test_spring_at_node_2_only(self): + # Regression test: a spring at node 2 alone used to leave + # matrix[2][2] untouched and corrupt the off-diagonal terms instead. + EI, l = 15000.0, 4.0 + k2 = 2000.0 + + matrix = constitutive_matrix( + EA=1e9, + EI=EI, + l=l, + spring={2: k2}, + node_1_hinge=False, + node_2_hinge=False, + ) + + c22_expected = 1 / (1 / (4 * EI / l) + 1 / k2) + + np.testing.assert_allclose(matrix[1][1], 4 * EI / l) + np.testing.assert_allclose(matrix[2][2], c22_expected) + np.testing.assert_allclose(matrix[1][2], -2 * EI / l) + np.testing.assert_allclose(matrix[2][1], -2 * EI / l) + + +if __name__ == "__main__": + unittest.main()