Trusses with All Elements Hinged Fail with StabilityError

[smith120bh]

The following script for a simple king post truss fails with error anastruct.basic.FEMException: ('StabilityError', 'The eigenvalues of the stiffness matrix are non zero, which indicates a instable structure. Check your support conditions'):

Example Script

from anastruct.fem.system import SystemElements

sys = SystemElements(mesh=50)
sys.add_element(
    location=[[0.0, 0.0], [2.5, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[0.0, 0.0], [2.5, 2.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[2.5, 0.0], [5.0, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[2.5, 2.0], [2.5, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_element(
    location=[[2.5, 2.0], [5.0, 0.0]], EA=68300, EI=128, g=0, spring={1: 0, 2: 0}
)
sys.add_support_hinged(node_id=1)
sys.add_support_hinged(node_id=4)
sys.point_load(Fx=0, Fy=-20.0, node_id=3)
sys.solve()
sys.show_results()

Discussion

I see that the function ensure_single_hinge() tries to ensure that no node ends up with an unrestrained rotational degree of freedom, by ensuring that all nodes have at least one element with a non-hinged connection to the node. However, that function is broken due to a minor syntax error:
if spring is not None and 0 in spring:
should instead read
if spring is not None and 0 in spring.values():

I've tested a fix for this (along with a couple consequential downstream fixes that come up when this function is fixed). See below diffs to util.py, system.py, node.py, and postprocess.py

Diff of Proposed Fix

diff --git a/anastruct/fem/node.py b/anastruct/fem/node.py
index f905ca0..9dea756 100644
--- a/anastruct/fem/node.py
+++ b/anastruct/fem/node.py
@@ -1,5 +1,5 @@
 class Node:
-    def __init__(self, id, Fx=0, Fz=0, Ty=0, ux=0, uz=0, phi_y=0, vertex=None):
+    def __init__(self, id, Fx=0, Fz=0, Ty=0, ux=0, uz=0, phi_y=0, vertex=None, hinge=False):
         """
         :param id: ID of the node, integer
         :param Fx: Value of Fx
@@ -20,7 +20,7 @@ class Node:
         self.uz = uz
         self.phi_y = phi_y
         self.vertex = vertex
-        self.hinge = False
+        self.hinge = hinge
         self.elements = {}
 
     @property
diff --git a/anastruct/fem/postprocess.py b/anastruct/fem/postprocess.py
index d6b2d93..8750ee3 100644
--- a/anastruct/fem/postprocess.py
+++ b/anastruct/fem/postprocess.py
@@ -20,7 +20,7 @@ class SystemLevel:
 
         for el in self.system.element_map.values():
             # post processor element level
-            self.post_el.node_results(el)
+            self.post_el.node_results(self, el)
 
     def node_results_system(self):
         for k, v in self.system.node_element_map.items():
@@ -88,11 +88,15 @@ class ElementLevel:
         self.system = system
 
     @staticmethod
-    def node_results(element):
+    def node_results(self, element):
         """
         Determine node results on the element level.
 
         """
+        # Check for hinges
+        hinge1 = self.system.node_map[element.node_id1].hinge
+        hinge2 = self.system.node_map[element.node_id2].hinge
+        
         # Global coordinates system
         element.node_map[element.node_id1] = Node(
             id=element.node_id1,
@@ -101,10 +105,15 @@ class ElementLevel:
             Fz=element.element_force_vector[1]
             + element.element_primary_force_vector[1],
             Ty=element.element_force_vector[2]
-            + element.element_primary_force_vector[2],
+            + element.element_primary_force_vector[2]
+            if not hinge1
+            else 0,
             ux=element.element_displacement_vector[0],
             uz=element.element_displacement_vector[1],
-            phi_y=element.element_displacement_vector[2],
+            phi_y=element.element_displacement_vector[2]
+            if not hinge1
+            else 0,
+            hinge=hinge1,
         )
 
         element.node_map[element.node_id2] = Node(
@@ -114,10 +123,15 @@ class ElementLevel:
             Fz=element.element_force_vector[4]
             + element.element_primary_force_vector[4],
             Ty=element.element_force_vector[5]
-            + element.element_primary_force_vector[5],
+            + element.element_primary_force_vector[5]
+            if not hinge2
+            else 0,
             ux=element.element_displacement_vector[3],
             uz=element.element_displacement_vector[4],
-            phi_y=element.element_displacement_vector[5],
+            phi_y=element.element_displacement_vector[5]
+            if not hinge2
+            else 0,
+            hinge=hinge2,
         )
 
         # Local coordinate system. With inclined supports
diff --git a/anastruct/fem/system.py b/anastruct/fem/system.py
index 62d240c..41451c5 100644
--- a/anastruct/fem/system.py
+++ b/anastruct/fem/system.py
@@ -248,7 +248,7 @@ class SystemElements:
         system_components.util.append_node_id(
             self, point_1, point_2, node_id1, node_id2
         )
-        system_components.util.ensure_single_hinge(self, spring, node_id1, node_id2)
+        spring = system_components.util.ensure_single_hinge(self, spring, node_id1, node_id2)
 
         # add element
         element = Element(
@@ -634,7 +634,7 @@ class SystemElements:
 
         for id_ in node_id:
             id_ = _negative_index_to_id(id_, self.node_map.keys())
-            system_components.util.support_check(self, id_)
+            #system_components.util.support_check(self, id_)
 
             # add the support to the support list for the plotter
             self.supports_hinged.append(self.node_map[id_])
@@ -653,7 +653,7 @@ class SystemElements:
 
         for id_ in node_id:
             id_ = _negative_index_to_id(id_, self.node_map.keys())
-            system_components.util.support_check(self, id_)
+            #system_components.util.support_check(self, id_)
 
             if direction == "x":
                 direction = 2
diff --git a/anastruct/fem/system_components/util.py b/anastruct/fem/system_components/util.py
index e757a7e..8e2e212 100644
--- a/anastruct/fem/system_components/util.py
+++ b/anastruct/fem/system_components/util.py
@@ -4,31 +4,68 @@ from anastruct.basic import FEMException, angle_x_axis
 
 
 def ensure_single_hinge(system, spring, node_id1, node_id2):
-    if spring is not None and 0 in spring:
+    if spring is not None and 0 in spring.values():
         """
         Must be one rotational fixed element per node. Thus keep track of the hinges (k == 0).
         """
 
         for node in range(1, 3):
-            if spring[node] == 0:  # node is a hinged node
-                if node == 1:
-                    node_id = node_id1
-                else:
-                    node_id = node_id2
-
-                system.node_map[node_id].hinge = True
-
+            if node == 1:
+                node_id = node_id1
+            else:
+                node_id = node_id2
+            
+            if node in spring.keys() and spring[node] == 0:  # node is a hinged node
                 if len(system.node_map[node_id].elements) > 0:
-                    pass_hinge = not all(
-                        [
-                            el.hinge == node
-                            for el in system.node_map[node_id].elements.values()
-                        ]
-                    )
+                    hinges = []
+                    for el in system.node_map[node_id].elements.values():
+                        if node_id == el.node_id1:
+                            hinges.append(1 in el.springs.keys() and el.springs[1] == 0)
+                        elif node_id == el.node_id2:
+                            hinges.append(2 in el.springs.keys() and el.springs[2] == 0)
+                    pass_hinge = not all(hinges)
                 else:
                     pass_hinge = True
                 if not pass_hinge:
+                    system.node_map[node_id].hinge = True
                     del spring[node]  # too many hinges at that element.
+            
+            elif node not in spring.keys() or (node in spring.keys() and spring[node] != 0):
+                """
+                If a fixed element is added after a hinged element,
+                then add the removed spring release back in and set node as not hinged
+                """
+                system.node_map[node_id].hinge = False
+                if len(system.node_map[node_id].elements) > 0:
+                    for el in system.node_map[node_id].elements.values():
+                        if node_id == el.node_id1:
+                            system.element_map[el.id].springs.update({
+                                1: 0
+                            })
+                        elif node_id == el.node_id2:
+                            system.element_map[el.id].springs.update({
+                                2: 0
+                            })
+    else:
+        """
+        If a fixed element is added after a hinged element,
+        then add the removed spring release back in and set node as not hinged
+        """
+        if system.node_map[node_id1].hinge:
+            system.node_map[node_id1].hinge = False
+            for el in system.node_map[node_id1].elements.values():
+                system.element_map[el.id].springs.update({
+                    1: 0
+                })
+
+        if system.node_map[node_id2].hinge:
+            system.node_map[node_id2].hinge = False
+            for el in system.node_map[node_id2].elements.values():
+                system.element_map[el.id].springs.update({
+                    2: 0
+                })
+    
+    return spring
 
 
 def append_node_id(self, point_1, point_2, node_id1, node_id2):
@@ -75,7 +112,7 @@ def det_node_ids(system, point_1, point_2):
 def support_check(system, node_id):
     if system.node_map[node_id].hinge:
         raise FEMException(
-            "Flawed inputs", "You cannot add a support to a hinged node."
+            "Flawed inputs", "You cannot add a rotation-restraining support to a hinged node."
         )

[ritchie46]

Ah right, that seems like an edge case I haven't properly tested. I will take a look at your PR. Seems valid so far.

[ritchie46]

@smith120bh , is there any news on this?

[smith120bh]

@ritchie46 : So the diff above is complete and works from my testing. I am unable to create a pull request, as I get an error saying that I don't have write access to your repository any time I try!

[ritchie46]

@ritchie46 : So the diff above is complete and works from my testing. I am unable to create a pull request, as I get an error saying that I don't have write access to your repository any time I try!

Yeah you need to create a Pull Request from a fork from your repo. Could you try that? That saves me some copy pasting. ;)

[smith120bh]

Ah, thanks! I've done a lot of work on Github, but have never actually wanted to change more than a couple lines of code in someone else's repository before... But now done! You should be able to test and eventually merge that PR now :-)