From 2f3ec33bf91c482aa87aad64e3080d1f081a52df Mon Sep 17 00:00:00 2001 From: Vincent Gao Date: Mon, 1 Jun 2026 17:12:35 +0200 Subject: [PATCH] Assemble system matrix over actual elements, not contiguous ids insert_node() removes the split element and appends two new ones, leaving element_map with non-contiguous keys (e.g. {1, 3, 4}). assemble_system_matrix iterated 'for i in range(len(element_map)): element_map[i + 1]', which assumed ids 1..N and raised KeyError during validation/solve. Iterate over element_map.values() instead; assembly is an order-independent accumulation into the system matrix indexed by node id, so this is equivalent for the contiguous case and correct when ids have gaps. Fixes #308 --- anastruct/fem/system_components/assembly.py | 7 +++- tests/test_e2e.py | 45 +++++++++++++++++++++ 2 files changed, 50 insertions(+), 2 deletions(-) diff --git a/anastruct/fem/system_components/assembly.py b/anastruct/fem/system_components/assembly.py index 7cf830ec..416bdc3b 100644 --- a/anastruct/fem/system_components/assembly.py +++ b/anastruct/fem/system_components/assembly.py @@ -223,8 +223,11 @@ def assemble_system_matrix( # # thus with appending numbers in the system matrix: column = row - for i in range(len(system.element_map)): - element = system.element_map[i + 1] + # Iterate over the actual elements rather than assuming the element ids are + # a contiguous 1..N sequence: remove_element (e.g. via insert_node) can leave + # gaps in element_map, and assembly is an order-independent accumulation into + # the system matrix indexed by node id (see #308). + for element in system.element_map.values(): element_matrix = element.stiffness_matrix # n1 and n2 are starting indexes of the rows and the columns for node 1 and node 2 diff --git a/tests/test_e2e.py b/tests/test_e2e.py index d753523c..4c5b3b24 100644 --- a/tests/test_e2e.py +++ b/tests/test_e2e.py @@ -266,6 +266,51 @@ def it_retains_supports_and_loads(SS_20): assert SS_20.element_map[4].q_load == approx((1.6, 2.205258)) assert SS_20.element_map[5].q_load == approx((2.205258, 3)) + def context_insert_node_then_solve(): + # insert_node deletes the split element and appends two new ones, leaving + # element_map with non-contiguous ids. The structure must still assemble + # and solve, giving the same result as the equivalent model built with + # that node defined up front (#308). + inserted = SystemElements(EA=15000, EI=5000, mesh=101) + inserted.add_element(location=[[0, 0], [10, 4]]) + inserted.add_element(location=[[10, 4], [10, 0]]) + inserted.add_support_hinged(node_id=1) + inserted.add_support_hinged(node_id=3) + inserted.q_load(q=5, element_id=1, direction="element") + inserted.point_load(node_id=2, Fx=4, Fy=0) + inserted.insert_node(element_id=2, factor=0.25) + inserted.point_load(node_id=4, Fx=8, Fy=0) + inserted.solve() + + unsplit = SystemElements(EA=15000, EI=5000, mesh=101) + unsplit.add_element(location=[[0, 0], [10, 4]]) + unsplit.add_element(location=[[10, 4], [10, 3]]) + unsplit.add_element(location=[[10, 3], [10, 0]]) + unsplit.add_support_hinged(node_id=1) + unsplit.add_support_hinged(node_id=4) + unsplit.q_load(q=5, element_id=1, direction="element") + unsplit.point_load(node_id=2, Fx=4, Fy=0) + unsplit.point_load(node_id=3, Fx=8, Fy=0) + unsplit.solve() + + # the element ids are non-contiguous after the split + assert list(inserted.element_map.keys()) == [1, 3, 4] + + def displacements_by_coordinate(system): + return { + (round(node.vertex.x, 6), round(node.vertex.y, 6)): ( + system.get_node_displacements(node_id)["ux"], + system.get_node_displacements(node_id)["uy"], + ) + for node_id, node in system.node_map.items() + } + + inserted_disp = displacements_by_coordinate(inserted) + unsplit_disp = displacements_by_coordinate(unsplit) + assert set(inserted_disp) == set(unsplit_disp) + for coordinate, displacement in inserted_disp.items(): + assert displacement == approx(unsplit_disp[coordinate]) + def context_find_node_id(): # find_node_id() function using Example 8