From cb6d6e31d348ed5ff2c18a1d1a8ae132bdea3282 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Mon, 10 Jul 2023 20:01:52 -0400
Subject: [PATCH 01/20] Option for not rendering nodes added

---
 PyNite/Visualization.py | 117 ++++++++++++++++++++++------------------
 README.md               |   3 ++
 2 files changed, 68 insertions(+), 52 deletions(-)

diff --git a/PyNite/Visualization.py b/PyNite/Visualization.py
index 4f9f1828..dd47285c 100644
--- a/PyNite/Visualization.py
+++ b/PyNite/Visualization.py
@@ -23,6 +23,7 @@ def __init__(self, model):
         self.annotation_size = 5
         self.deformed_shape = False
         self.deformed_scale = 30
+        self.render_nodes = True
         self.render_loads = True
         self.color_map = None
         self.combo_name = 'Combo 1'
@@ -64,6 +65,9 @@ def set_deformed_shape(self, deformed_shape=False):
     def set_deformed_scale(self, scale=30):
         self.deformed_scale = scale
     
+    def set_render_nodes(self, render_nodes=True):
+        self.render_nodes = render_nodes
+
     def set_render_loads(self, render_loads=True):
         self.render_loads = render_loads
     
@@ -208,15 +212,18 @@ def update(self, reset_camera=True):
             self.render_loads = False
             warnings.warn('Unable to render load combination. No load combinations defined.', UserWarning)
         
-        # Create a visual node for each node in the model
-        vis_nodes = []
-        for node in self.model.Nodes.values():
-            vis_nodes.append(VisNode(node, self.annotation_size))
-        
-        # Create a visual auxiliary node for each auxiliary node in the model
-        vis_aux_nodes = []
-        for aux_node in self.model.AuxNodes.values():
-            vis_aux_nodes.append(VisNode(aux_node, self.annotation_size, color='red'))
+        # Check if nodes are to be rendered
+        if self.render_nodes == True:
+
+            # Create a visual node for each node in the model
+            vis_nodes = []
+            for node in self.model.Nodes.values():
+                vis_nodes.append(VisNode(node, self.annotation_size))
+            
+            # Create a visual auxiliary node for each auxiliary node in the model
+            vis_aux_nodes = []
+            for aux_node in self.model.AuxNodes.values():
+                vis_aux_nodes.append(VisNode(aux_node, self.annotation_size, color='red'))
         
         # Create a visual spring for each spring in the model
         vis_springs = []
@@ -262,52 +269,55 @@ def update(self, reset_camera=True):
                 # Set the text to follow the camera as the user interacts. This will
                 # require a reset of the camera (see below)
                 vis_member.lblActor.SetCamera(renderer.GetActiveCamera())
+        
+        # Check if nodes are to be rendered
+        if self.render_nodes == True:
+            
+            # Combine the polydata from each node
 
-        # Combine the polydata from each node
-
-        # Create an append filter for combining node polydata
-        node_polydata = vtk.vtkAppendPolyData()
+            # Create an append filter for combining node polydata
+            node_polydata = vtk.vtkAppendPolyData()
 
-        for vis_node in vis_nodes:
-            
-            # Add the node's polydata
-            node_polydata.AddInputData(vis_node.polydata.GetOutput())
+            for vis_node in vis_nodes:
+                
+                # Add the node's polydata
+                node_polydata.AddInputData(vis_node.polydata.GetOutput())
 
-            if self.labels == True:
+                if self.labels == True:
+                    
+                    # Add the actor for the node label
+                    renderer.AddActor(vis_node.lblActor)
                 
-                # Add the actor for the node label
-                renderer.AddActor(vis_node.lblActor)
+                    # Set the text to follow the camera as the user interacts. This will
+                    # require a reset of the camera (see below)
+                    vis_node.lblActor.SetCamera(renderer.GetActiveCamera())
             
-                # Set the text to follow the camera as the user interacts. This will
-                # require a reset of the camera (see below)
-                vis_node.lblActor.SetCamera(renderer.GetActiveCamera())
-        
-        # Update the node polydata in the append filter
-        node_polydata.Update()
-        
-        # Create a mapper and actor for the nodes
-        node_mapper = vtk.vtkPolyDataMapper()
-        node_mapper.SetInputConnection(node_polydata.GetOutputPort())
-        node_actor = vtk.vtkActor()
-        node_actor.SetMapper(node_mapper)
-        
-        # Add the node actor to the renderer
-        renderer.AddActor(node_actor)
+            # Update the node polydata in the append filter
+            node_polydata.Update()
+            
+            # Create a mapper and actor for the nodes
+            node_mapper = vtk.vtkPolyDataMapper()
+            node_mapper.SetInputConnection(node_polydata.GetOutputPort())
+            node_actor = vtk.vtkActor()
+            node_actor.SetMapper(node_mapper)
+            
+            # Add the node actor to the renderer
+            renderer.AddActor(node_actor)
 
-        # Add actors for each auxiliary node
-        for vis_aux_node in vis_aux_nodes:
-        
-            # Add the actor for the auxiliary node
-            renderer.AddActor(vis_aux_node.actor)
+            # Add actors for each auxiliary node
+            for vis_aux_node in vis_aux_nodes:
+            
+                # Add the actor for the auxiliary node
+                renderer.AddActor(vis_aux_node.actor)
 
-            if self.labels == True:
-                
-                # Add the actor for the auxiliary node label
-                renderer.AddActor(vis_aux_node.lblActor)
-        
-                # Set the text to follow the camera as the user interacts. This will
-                # require a reset of the camera (see below)
-                vis_aux_node.lblActor.SetCamera(renderer.GetActiveCamera())
+                if self.labels == True:
+                    
+                    # Add the actor for the auxiliary node label
+                    renderer.AddActor(vis_aux_node.lblActor)
+            
+                    # Set the text to follow the camera as the user interacts. This will
+                    # require a reset of the camera (see below)
+                    vis_aux_node.lblActor.SetCamera(renderer.GetActiveCamera())
 
         # Render the deformed shape if requested
         if self.deformed_shape == True:
@@ -1475,11 +1485,14 @@ def _DeformedShape(model, renderer, scale_factor, annotation_size, combo_name):
     # Create an append filter to add all the shape polydata to
     append_filter = vtk.vtkAppendPolyData()
     
-    # Add the deformed nodes to the append filter
-    for node in model.Nodes.values():
+    # Check if nodes are to be rendered
+    if renderer.render_nodes == True:
         
-        vis_node = VisDeformedNode(node, scale_factor, annotation_size, combo_name)
-        append_filter.AddInputData(vis_node.source.GetOutput())
+        # Add the deformed nodes to the append filter
+        for node in model.Nodes.values():
+            
+            vis_node = VisDeformedNode(node, scale_factor, annotation_size, combo_name)
+            append_filter.AddInputData(vis_node.source.GetOutput())
         
     # Add the springs to the append filter
     for spring in model.Springs.values():
diff --git a/README.md b/README.md
index 7f747b71..5604924d 100644
--- a/README.md
+++ b/README.md
@@ -61,6 +61,9 @@ Here's a list of projects that use PyNite:
 * Phaenotyp (https://github.com/bewegende-Architektur/Phaenotyp) (https://youtu.be/shloSw9HjVI)
 
 # What's New?
+v0.0.79
+* Add the option to turn off nodes during visualization.
+
 v0.0.78
 * Corrections to tension/compression only support springs. v0.0.76 and v0.0.77 were not working as expected. 3rd time's a charm.
 

From 54cbb1ed3d3c0d882f64a4b937d605cd9906af60 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 15 Jul 2023 16:49:11 -0400
Subject: [PATCH 02/20] Fix for cylinder meshes about X and Z axes

---
 PyNite/Mesh.py | 51 +++++++++++++++++++++++++++-----------------------
 1 file changed, 28 insertions(+), 23 deletions(-)

diff --git a/PyNite/Mesh.py b/PyNite/Mesh.py
index e864fb79..eec1f9f0 100644
--- a/PyNite/Mesh.py
+++ b/PyNite/Mesh.py
@@ -1203,7 +1203,7 @@ class CylinderMesh(Mesh):
     A mesh of quadrilaterals forming a cylinder.
 
     The mesh is formed with the local y-axis of the elements pointed toward
-    the base of the 
+    the base of the cylinder 
 
     Parameters
     ----------
@@ -1240,40 +1240,53 @@ class CylinderMesh(Mesh):
         The type of element to use for the mesh: 'Quad' or 'Rect'
     """
 
-    def __init__(self, mesh_size, radius, height, thickness, material, model, kx_mod=1, ky_mod=1,
-                 origin=[0, 0, 0], axis='Y', start_node='N1', start_element='Q1',
-                 num_elements=None, element_type='Quad'):
+    def __init__(self, mesh_size, radius, height, thickness, material, model, kx_mod=1, ky_mod=1,origin=[0, 0, 0], axis='Y', start_node='N1', start_element='Q1', num_elements=None, element_type='Quad'):
 
+        # Inherit properties and methods from the parent `Mesh` class
         super().__init__(thickness, material, model, kx_mod, ky_mod, start_node, start_element)
 
+        # Define a few new additional class properties related to cylinders
         self.radius = radius
         self.h = height
         self.mesh_size = mesh_size
+        self.origin = origin
+        self.axis = axis
 
+        # Check if the user has requested a specific number of elements for each course of plates. This can be useful for ensuring the mesh matches up with other meshes.
         if num_elements == None:
+            # Calculate the number of elements if the user hasn't specified
             self.num_elements = int(round(2*pi*radius/mesh_size, 0))
         else:
+            # Use the user specified number of elements
             self.num_elements = num_elements
-        
-        self.origin = origin
-        self.axis = axis
 
+        # Check which type of element the user has requested (rectangular plate or quad)
         self.element_type = element_type
 
+        # Generate the mesh
         self.generate()
     
     def generate(self):
         
+        # Get the mesh thickness and the material name
         thickness = self.thickness
         material = self.material
 
         mesh_size = self.mesh_size  # Desired mesh size
-        num_elements = self.num_elements  # Number of quadrilaterals in the ring
-        n = self.num_elements
+        num_elements = self.num_elements  # Number of quadrilaterals in each course of the ring
+        n = self.num_elements  # Total number of elements in the mesh (initialized for a single ring at the moment)
 
         radius = self.radius
         h = self.h
-        y = self.origin[1]
+
+        # Set the cylinder base's local y-coordinate
+        if self.axis == 'Y':
+            y = self.origin[1]  
+        elif self.axis == 'X':
+            y = self.origin[0]  
+        elif self.axis == 'Z':
+             y = self.origin[2]
+
         n = int(self.start_node[1:])
         q = int(self.start_element[1:])
 
@@ -1292,31 +1305,23 @@ def generate(self):
         
             # Create a mesh of nodes for the ring
             if self.axis == 'Y':
-                ring = CylinderRingMesh(radius, h_y, num_elements, thickness, material, self.model, 1, 1, [0, y, 0],
-                                        self.axis, 'N' + str(n), 'Q' + str(q), element_type)
+                ring = CylinderRingMesh(radius, h_y, num_elements, thickness, material, self.model, 1, 1, [0, y, 0], self.axis, 'N' + str(n), 'Q' + str(q), element_type)
             elif self.axis == 'X':
-                ring = CylinderRingMesh(radius, h_y, num_elements, thickness, material, self.model, 1, 1, [y, 0, 0],
-                                        self.axis, 'N' + str(n), 'Q' + str(q), element_type)
+                ring = CylinderRingMesh(radius, h_y, num_elements, thickness, material, self.model, 1, 1, [y, 0, 0], self.axis, 'N' + str(n), 'Q' + str(q), element_type)
             elif self.axis == 'Z':
-                ring = CylinderRingMesh(radius, h_y, num_elements, thickness, material, self.model, 1, 1, [0, 0, y],
-                                        self.axis, 'N' + str(n), 'Q' + str(q), element_type)
+                ring = CylinderRingMesh(radius, h_y, num_elements, thickness, material, self.model, 1, 1, [0, 0, y], self.axis, 'N' + str(n), 'Q' + str(q), element_type)
 
             n += num_elements
             q += num_elements
         
-            # Add the newly generated nodes and elements to the overall mesh. Note that if duplicate
-            # keys exist, the `.update()` method will overwrite them with the newly generated key value
-            # pairs. This works in our favor by automatically eliminating duplicate nodes at the shared
-            # boundaries between rings.
+            # Add the newly generated nodes and elements to the overall mesh. Note that if duplicate keys exist, the `.update()` method will overwrite them with the newly generated key value pairs. This works in our favor by automatically eliminating duplicate nodes at the shared boundaries between rings.
             self.nodes.update(ring.nodes)
             self.elements.update(ring.elements)
 
             # Prepare to move to the next ring
             y += h_y
         
-        # After calling the `.update()` method some elements are still attached to the duplicate
-        # nodes that are no longer in the dictionary. Attach these plates to the nodes that are
-        # still in the dictionary instead. 
+        # After calling the `.update()` method some elements are still attached to the duplicate nodes that are no longer in the dictionary. Attach these plates to the nodes that are still in the dictionary instead. 
         for element in self.elements.values():
             element.i_node = self.nodes[element.i_node.name]
             element.j_node = self.nodes[element.j_node.name]

From c63780a66e04af39343af13fe09207ead2573551 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 15 Jul 2023 17:16:21 -0400
Subject: [PATCH 03/20] Update README.md

---
 README.md | 26 +++-----------------------
 1 file changed, 3 insertions(+), 23 deletions(-)

diff --git a/README.md b/README.md
index 5604924d..15e4522e 100644
--- a/README.md
+++ b/README.md
@@ -62,7 +62,8 @@ Here's a list of projects that use PyNite:
 
 # What's New?
 v0.0.79
-* Add the option to turn off nodes during visualization.
+* Added the option to turn off nodes during visualization.
+* Bug fix for meshing cylinders about the global X or Z axis.
 
 v0.0.78
 * Corrections to tension/compression only support springs. v0.0.76 and v0.0.77 were not working as expected. 3rd time's a charm.
@@ -99,25 +100,4 @@ v0.0.71
 * Simplified internal code for finding unique names for objects.
 
 v0.0.70
-* Array output of member force diagrams and displacement diagrams has been added.
-
-v0.0.69
-* Bug fix for rotational springs. Exceptions were being thrown due to an inconsistent variable name.
-* Cleared out old branches from the repository that were no longer being used.
-* Updated CI to check against python 3.10 and 3.11. Removed CI for python 3.6 as it's no longer supported by the latest version of github actions.
-* Subtle changes to the logo to make it look a little more "pythonic".
-* Bug fix for rendering screenshots. The ability to interact with the render window was being disabled after the first screenshot had been taken, forcing subsequent screenshots to use the same view as the first one.
-
-v0.0.68
-* Bug fix for member distributed loads on physical members. Added a unit test to check for this error going forward. This bug only affected physical members (new as of v0.0.67) that had distributed loads and internal nodes.
-
-v0.0.67
-* Added physical members. Members now automatically detect internal nodes and subdivide themselves and their loads.
-* Refactoring: deprecated old method names for member results. You may now have some errors show up if you still try to get member results using the old method names.
-* Bug fix for P-Delta analysis. Global displacements were correct, but member internal forces were neglecting the geometric stiffness matrix. The impact of this bug was minimal, since the strain induced by correct global displacements was still being considered prior to this update. You should see a slight change to member P-Delta results.
-* Code simplification for P-Delta analysis.
-
-v0.0.66
-* Code simplification and bug fix for merging duplicate nodes.
-* When nodes are merged, support conditions for the deleted node are now assigned to the remaining node.
-* Added a linear solver for faster analysis of simple models. If you don't need P-Delta analysis or tension/compression-only analysis this solver saves time by only assembling the global stiffness matrix once.
+* Array output of member force diagrams and displacement diagrams has been added.
\ No newline at end of file

From 43a8db7a2f7017614fc7063c1c63d0d3eb0d14ad Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 15 Jul 2023 17:43:40 -0400
Subject: [PATCH 04/20] Added 'default' theme to rendering

---
 PyNite/Visualization.py | 68 ++++++++++++++++++++++++++---------------
 1 file changed, 43 insertions(+), 25 deletions(-)

diff --git a/PyNite/Visualization.py b/PyNite/Visualization.py
index dd47285c..6f042c31 100644
--- a/PyNite/Visualization.py
+++ b/PyNite/Visualization.py
@@ -31,6 +31,7 @@ def __init__(self, model):
         self.labels = True
         self.scalar_bar = False
         self.scalar_bar_text_size = 24
+        self.theme = 'default'
 
         # Initialize VTK objects
         self.renderer = vtk.vtkRenderer()
@@ -215,15 +216,18 @@ def update(self, reset_camera=True):
         # Check if nodes are to be rendered
         if self.render_nodes == True:
 
+            if self.theme == 'print': color = 'black'
+            else: color = None
+
             # Create a visual node for each node in the model
             vis_nodes = []
             for node in self.model.Nodes.values():
-                vis_nodes.append(VisNode(node, self.annotation_size))
+                vis_nodes.append(VisNode(node, self.annotation_size, color))
             
             # Create a visual auxiliary node for each auxiliary node in the model
             vis_aux_nodes = []
             for aux_node in self.model.AuxNodes.values():
-                vis_aux_nodes.append(VisNode(aux_node, self.annotation_size, color='red'))
+                vis_aux_nodes.append(VisNode(aux_node, self.annotation_size, color))
         
         # Create a visual spring for each spring in the model
         vis_springs = []
@@ -325,7 +329,7 @@ def update(self, reset_camera=True):
 
         # Render the loads if requested
         if (self.combo_name != None or self.case != None) and self.render_loads != False:
-            _RenderLoads(self.model, renderer, self.annotation_size, self.combo_name, self.case)
+            _RenderLoads(self.model, renderer, self.annotation_size, self.combo_name, self.case, self.theme)
         
         # Render the plates and quads, if present
         if self.model.Quads or self.model.Plates:
@@ -347,9 +351,7 @@ def RenderModel(model, annotation_size=5, deformed_shape=False, deformed_scale=3
     render_model(model, annotation_size, deformed_shape, deformed_scale, render_loads, color_map,
                  True, combo_name, case, labels, screenshot)
 
-def render_model(model, annotation_size=5, deformed_shape=False, deformed_scale=30,
-                 render_loads=True, color_map=None, scalar_bar=True, combo_name='Combo 1', case=None, labels=True,
-                 screenshot=None):
+def render_model(model, annotation_size=5, deformed_shape=False, deformed_scale=30, render_loads=True, color_map=None, scalar_bar=True, combo_name='Combo 1', case=None, labels=True, screenshot=None, theme='default'):
     '''
     Renders a finite element model using VTK.
     
@@ -848,18 +850,21 @@ def __init__(self, node, annotation_size=5, color=None):
       
         # Add color to the actors
         if color == 'red':
-            self.actor.GetProperty().SetColor(255, 0, 0) # Red
-            self.lblActor.GetProperty().SetColor(255, 0, 0) # Red
+            self.actor.GetProperty().SetColor(255, 0, 0)  # Red
+            self.lblActor.GetProperty().SetColor(255, 0, 0)  # Red
         elif color == 'yellow':
-            self.actor.GetProperty().SetColor(255, 255, 0) # Yellow
-            self.lblActor.GetProperty().SetColor(255, 255, 0) # Yellow
+            self.actor.GetProperty().SetColor(255, 255, 0)  # Yellow
+            self.lblActor.GetProperty().SetColor(255, 255, 0)  # Yellow
+        elif color == 'black':
+            self.actor.GetProperty().SetColor(0, 0, 0)  # Black
+            self.lblActor.GetProperty().SetColor(0, 0, 0)  # Black
         
         # Set the mapper for the node's actor
         self.actor.SetMapper(mapper)
 
 class VisSpring():
     
-    def __init__(self, spring, nodes, annotation_size=5):
+    def __init__(self, spring, nodes, annotation_size=5, color=None):
 
         # Generate a line source for the spring
         line = vtk.vtkLineSource()
@@ -888,7 +893,8 @@ def __init__(self, spring, nodes, annotation_size=5):
 
         # Set up an actor for the spring
         self.actor = vtk.vtkActor()
-        self.actor.GetProperty().SetColor(255, 0, 255) # Magenta
+        if color is None: self.actor.GetProperty().SetColor(255, 0, 255)  # Magenta
+        elif color == 'black': self.actor.GetProperty().SetColor(0, 0, 0)  # Black
         self.actor.SetMapper(mapper)
 
         # Create the text for the spring label
@@ -909,7 +915,7 @@ def __init__(self, spring, nodes, annotation_size=5):
 class VisMember():
 
     # Constructor
-    def __init__(self, member, nodes, annotation_size=5):
+    def __init__(self, member, nodes, annotation_size=5, color=None):
     
         # Generate a line for the member
         line = vtk.vtkLineSource()
@@ -953,6 +959,11 @@ def __init__(self, member, nodes, annotation_size=5):
         self.lblActor.SetScale(annotation_size, annotation_size, annotation_size)
         self.lblActor.SetPosition((Xi+Xj)/2, (Yi+Yj)/2, (Zi+Zj)/2)
 
+        # Adjust the color of the member
+        if color == 'black':
+            self.actor.GetProperty().SetColor(0, 0, 0)  # Black
+            self.lblActor.GetProperty().SetColor(0, 0, 0) # Black
+
 # Converts a node object into a node in its deformed position for the viewer
 class VisDeformedNode():
     
@@ -1083,7 +1094,7 @@ class VisPtLoad():
     Creates a point load for the viewer
     '''
     
-    def __init__(self, position, direction, length, label_text=None, annotation_size=5):
+    def __init__(self, position, direction, length, label_text=None, annotation_size=5, color=None):
         '''
         Constructor.
       
@@ -1142,7 +1153,8 @@ def __init__(self, position, direction, length, label_text=None, annotation_size
         mapper = vtk.vtkPolyDataMapper()
         mapper.SetInputConnection(self.polydata.GetOutputPort())
         self.actor = vtk.vtkActor()
-        self.actor.GetProperty().SetColor(0, 255, 0) # Green
+        if color is None: self.actor.GetProperty().SetColor(0, 255, 0) # Green
+        elif color == 'black': self.actor.GetProperty().SetColor(0, 0, 0)  # Black
         self.actor.SetMapper(mapper)
       
         # Create the label if needed
@@ -1163,14 +1175,15 @@ def __init__(self, position, direction, length, label_text=None, annotation_size
             self.lblActor.SetPosition(position[0] - (length - 0.6*annotation_size)*unitVector[0], \
                                       position[1] - (length - 0.6*annotation_size)*unitVector[1], \
                                       position[2] - (length - 0.6*annotation_size)*unitVector[2])
-            self.lblActor.GetProperty().SetColor(0, 255, 0) # Green
+            if color is None: self.lblActor.GetProperty().SetColor(0, 255, 0)  # Green
+            elif color == 'black': self.lblActor.GetProperty().SetColor(0, 0, 0)  # Black
       
 class VisDistLoad():
     '''
     Creates a distributed load for the viewer
     '''
     
-    def __init__(self, position1, position2, direction, length1, length2, label_text1, label_text2, annotation_size=5):
+    def __init__(self, position1, position2, direction, length1, length2, label_text1, label_text2, annotation_size=5, color=None):
         '''
         Constructor.
         '''
@@ -1231,7 +1244,8 @@ def __init__(self, position1, position2, direction, length1, length2, label_text
         mapper = vtk.vtkPolyDataMapper()
         mapper.SetInputConnection(self.polydata.GetOutputPort())
         self.actor = vtk.vtkActor()
-        self.actor.GetProperty().SetColor(0, 255, 0) # Green
+        if color is None: self.actor.GetProperty().SetColor(0, 255, 0)  # Green
+        elif color == 'black': self.actor.GetProperty().SetColor(0, 0, 0)  # Black
         self.actor.SetMapper(mapper)
       
         # Get the actors for the labels
@@ -1242,7 +1256,7 @@ class VisMoment():
     Creates a concentrated moment for the viewer
     '''
     
-    def __init__(self, center, direction, radius, label_text=None, annotation_size=5):
+    def __init__(self, center, direction, radius, label_text=None, annotation_size=5, color=None):
         '''
         Constructor.
       
@@ -1305,7 +1319,8 @@ def __init__(self, center, direction, radius, label_text=None, annotation_size=5
         self.lblActor.SetPosition(Xc + v3[0]*(radius + 0.25*annotation_size), \
                                   Yc + v3[1]*(radius + 0.25*annotation_size), \
                                   Zc + v3[2]*(radius + 0.25*annotation_size))
-        self.lblActor.GetProperty().SetColor(0, 255, 0)  # Green
+        if color is None: self.lblActor.GetProperty().SetColor(0, 255, 0)  # Green
+        elif color == 'black': self.lblActor.GetProperty().SetColor(0, 0, 0)  # Black
 
 class VisAreaLoad():
     '''
@@ -1458,7 +1473,7 @@ def _PrepContour(model, stress_type='Mx', combo_name='Combo 1'):
             if node.contour != []:
                 node.contour = sum(node.contour)/len(node.contour)
 
-def _DeformedShape(model, renderer, scale_factor, annotation_size, combo_name):
+def _DeformedShape(model, renderer, scale_factor, annotation_size, combo_name, color=None):
     '''
     Renders the deformed shape of a model.
     
@@ -1516,11 +1531,12 @@ def _DeformedShape(model, renderer, scale_factor, annotation_size, combo_name):
     mapper = vtk.vtkPolyDataMapper()
     mapper.SetInputConnection(append_filter.GetOutputPort())
     actor = vtk.vtkActor()
-    actor.GetProperty().SetColor(255, 255, 0)  # Yellow
+    if color is None: actor.GetProperty().SetColor(255, 255, 0)  # Yellow
+    elif color == 'black': actor.GetProperty().SetColor(0, 0, 0)  # Black
     actor.SetMapper(mapper)
     renderer.AddActor(actor)
 
-def _RenderLoads(model, renderer, annotation_size, combo_name, case):
+def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='default'):
 
     # Create an append filter to store all the polydata in. This will allow us to use fewer actors to
     # display all the loads, which will greatly improve rendering speed as the user interacts. VTK
@@ -1741,7 +1757,8 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case):
     load_mapper = vtk.vtkPolyDataMapper()
     load_mapper.SetInputConnection(polydata.GetOutputPort())
     load_actor = vtk.vtkActor()
-    load_actor.GetProperty().SetColor(0, 255, 0)  # Green
+    if theme != 'print': load_actor.GetProperty().SetColor(0, 255, 0)  # Green
+    else: load_actor.GetProperty().SetColor(0, 0, 0)  # Black
     load_actor.SetMapper(load_mapper)
     renderer.AddActor(load_actor)
     
@@ -1749,7 +1766,8 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case):
     polygon_mapper = vtk.vtkPolyDataMapper()
     polygon_mapper.SetInputData(polygon_polydata)
     polygon_actor = vtk.vtkActor()
-    polygon_actor.GetProperty().SetColor(0, 255, 0)  # Green
+    if theme != 'print': polygon_actor.GetProperty().SetColor(0, 255, 0)  # Green
+    else: polygon_actor.GetProperty().SetColor(128, 128, 128)  # Grey
     # polygon_actor.GetProperty().SetOpacity(0.5)      # 50% opacity
     polygon_actor.SetMapper(polygon_mapper)
     renderer.AddActor(polygon_actor)

From bbff7a1f672397822fc418aa4bd06eb19774a068 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 15 Jul 2023 20:26:49 -0400
Subject: [PATCH 05/20] Move TC checks out of `FEModel3D`

---
 PyNite/Analysis.py  | 127 +++++++++++++++++++++
 PyNite/FEModel3D.py | 264 ++------------------------------------------
 2 files changed, 138 insertions(+), 253 deletions(-)
 create mode 100644 PyNite/Analysis.py

diff --git a/PyNite/Analysis.py b/PyNite/Analysis.py
new file mode 100644
index 00000000..d70845f1
--- /dev/null
+++ b/PyNite/Analysis.py
@@ -0,0 +1,127 @@
+
+def _check_TC_convergence(model, combo_name='Combo 1', log=True):
+    
+    # Assume the model has converged until we find out otherwise
+    convergence = True
+    
+    # Provide an update to the console if requested by the user
+    if log: print('- Checking for tension/compression-only support spring convergence')
+    for node in model.Nodes.values():
+
+        # Check convergence of tension/compression-only spring supports and activate/deactivate them as necessary
+        if node.spring_DX[1] is not None:
+            if ((node.spring_DX[1] == '-' and node.DX[combo_name] > 0)
+            or (node.spring_DX[1] == '+' and node.DX[combo_name] < 0)):
+                # Check if the spring is switching from active to inactive
+                if node.spring_DX[2] == True: convergence = False
+                # Make sure the spring is innactive
+                node.spring_DX[2] = False
+            elif ((node.spring_DX[1] == '-' and node.DX[combo_name] < 0)
+            or (node.spring_DX[1] == '+' and node.DX[combo_name] > 0)):
+                # Check if the spring is switching from inactive to active
+                if node.spring_DX[2] == False: convergence = False
+                # Make sure the spring is active
+                node.spring_DX[2] = True
+        if node.spring_DY[1] is not None:
+            if ((node.spring_DY[1] == '-' and node.DY[combo_name] > 0)
+            or (node.spring_DY[1] == '+' and node.DY[combo_name] < 0)):
+                # Check if the spring is switching from active to inactive
+                if node.spring_DY[2] == True: convergence = False
+                # Make sure the spring is innactive
+                node.spring_DY[2] = False
+            elif ((node.spring_DY[1] == '-' and node.DY[combo_name] < 0)
+            or (node.spring_DY[1] == '+' and node.DY[combo_name] > 0)):
+                # Check if the spring is switching from inactive to active
+                if node.spring_DY[2] == False: convergence = False
+                # Make sure the spring is active
+                node.spring_DY[2] = True
+        if node.spring_DZ[1] is not None:
+            if ((node.spring_DZ[1] == '-' and node.DZ[combo_name] > 0)
+            or (node.spring_DZ[1] == '+' and node.DZ[combo_name] < 0)):
+                # Check if the spring is switching from active to inactive
+                if node.spring_DZ[2] == True: convergence = False
+                # Make sure the spring is innactive
+                node.spring_DZ[2] = False
+            elif ((node.spring_DZ[1] == '-' and node.DZ[combo_name] < 0)
+            or (node.spring_DZ[1] == '+' and node.DZ[combo_name] > 0)):
+                # Check if the spring is switching from inactive to active
+                if node.spring_DZ[2] == False: convergence = False
+                # Make sure the spring is active
+                node.spring_DZ[2] = True
+        if node.spring_RX[1] is not None:
+            if ((node.spring_RX[1] == '-' and node.RX[combo_name] > 0)
+            or (node.spring_RX[1] == '+' and node.RX[combo_name] < 0)):
+                # Check if the spring is switching from active to inactive
+                if node.spring_RX[2] == True: convergence = False
+                # Make sure the spring is innactive
+                node.spring_RX[2] = False
+            elif ((node.spring_RX[1] == '-' and node.RX[combo_name] < 0)
+            or (node.spring_RX[1] == '+' and node.RX[combo_name] > 0)):
+                # Check if the spring is switching from inactive to active
+                if node.spring_RX[2] == False: convergence = False
+                # Make sure the spring is active
+                node.spring_RX[2] = True
+        if node.spring_RY[1] is not None:
+            if ((node.spring_RY[1] == '-' and node.RY[combo_name] > 0)
+            or (node.spring_RY[1] == '+' and node.RY[combo_name] < 0)):
+                # Check if the spring is switching from active to inactive
+                if node.spring_RY[2] == True: convergence = False
+                # Make sure the spring is innactive
+                node.spring_RY[2] = False
+            elif ((node.spring_RY[1] == '-' and node.RY[combo_name] < 0)
+            or (node.spring_RY[1] == '+' and node.RY[combo_name] > 0)):
+                # Check if the spring is switching from inactive to active
+                if node.spring_RY[2] == False: convergence = False
+                # Make sure the spring is active
+                node.spring_RY[2] = True
+        if node.spring_RZ[1] is not None:
+            if ((node.spring_RZ[1] == '-' and node.RZ[combo_name] > 0)
+            or (node.spring_RZ[1] == '+' and node.RZ[combo_name] < 0)):
+                # Check if the spring is switching from active to inactive
+                if node.spring_RZ[2] == True: convergence = False
+                # Make sure the spring is innactive
+                node.spring_RZ[2] = False
+            elif ((node.spring_RZ[1] == '-' and node.RZ[combo_name] < 0)
+            or (node.spring_RZ[1] == '+' and node.RZ[combo_name] > 0)):
+                # Check if the spring is switching from inactive to active
+                if node.spring_RZ[2] == False: convergence = False
+                # Make sure the spring is active
+                node.spring_RZ[2] = True
+    
+    # TODO: Adjust the code below to allow elements to reactivate on subsequent iterations if deformations at element nodes indicate the member goes back into an active state. This will lead to a less conservative and more realistic analysis. Nodal springs (above) already do this.
+
+    # Check tension/compression-only springs
+    if log: print('- Checking for tension/compression-only spring convergence')
+    for spring in model.Springs.values():
+
+        if spring.active[combo_name] == True:
+
+            # Check if tension-only conditions exist
+            if spring.tension_only == True and spring.axial(combo_name) > 0:
+                spring.active[combo_name] = False
+                convergence = False
+            
+            # Check if compression-only conditions exist
+            elif spring.comp_only == True and spring.axial(combo_name) < 0:
+                spring.active[combo_name] = False
+                convergence = False
+
+    # Check tension/compression only members
+    if log: print('- Checking for tension/compression-only member convergence')
+    for phys_member in model.Members.values():
+
+        # Only run the tension/compression only check if the member is still active
+        if phys_member.active[combo_name] == True:
+
+            # Check if tension-only conditions exist
+            if phys_member.tension_only == True and phys_member.max_axial(combo_name) > 0:
+                phys_member.active[combo_name] = False
+                convergence = False
+
+            # Check if compression-only conditions exist
+            elif phys_member.comp_only == True and phys_member.min_axial(combo_name) < 0:
+                phys_member.active[combo_name] = False
+                convergence = False
+
+    # Return whether the TC analysis has converged
+    return convergence
\ No newline at end of file
diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index 9d9ea9b2..3c854eda 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -16,6 +16,7 @@
 from PyNite.Plate3D import Plate3D
 from PyNite.LoadCombo import LoadCombo
 from PyNite.Mesh import Mesh, RectangleMesh, AnnulusMesh, FrustrumMesh, CylinderMesh
+from PyNite import Analysis
 
 # %%
 class FEModel3D():
@@ -1971,7 +1972,12 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
 
             # Iterate until convergence or divergence occurs
             while convergence == False and divergence == False:
-
+                
+                # Check for tension/compression-only divergence
+                if iter_count > max_iter:
+                    divergence = True
+                    raise Exception('Model diverged during tension/compression-only analysis')
+                
                 # Get the partitioned global stiffness matrix K11, K12, K21, K22
                 if sparse == True:
                     K11, K12, K21, K22 = self._partition(self.K(combo.name, log, check_stability, sparse).tolil(), D1_indices, D2_indices)
@@ -2052,131 +2058,9 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
                     node.RY[combo.name] = D[node.ID*6 + 4, 0]
                     node.RZ[combo.name] = D[node.ID*6 + 5, 0]
                 
-                # Check for divergence
-                if iter_count > max_iter:
-                    divergence = True
-                    raise Exception('Model diverged during tension/compression-only analysis')
-                
-                # Assume the model has converged (to be checked below)
-                convergence = True
+                # Check for tension/compression-only convergence
+                convergence = Analysis._check_TC_convergence(self, combo.name, log=log)
 
-                # Check tension/compression-only spring supports
-                if log: print('- Checking for tension/compression-only support spring convergence')
-                for node in self.Nodes.values():
-                    
-                    # Check convergence of tension/compression-only spring supports and activate/deactivate them as necessary
-                    if node.spring_DX[1] is not None:
-                        if ((node.spring_DX[1] == '-' and node.DX[combo.name] > 0)
-                        or (node.spring_DX[1] == '+' and node.DX[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_DX[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_DX[2] = False
-                        elif ((node.spring_DX[1] == '-' and node.DX[combo.name] < 0)
-                        or (node.spring_DX[1] == '+' and node.DX[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_DX[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_DX[2] = True
-                    if node.spring_DY[1] is not None:
-                        if ((node.spring_DY[1] == '-' and node.DY[combo.name] > 0)
-                        or (node.spring_DY[1] == '+' and node.DY[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_DY[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_DY[2] = False
-                        elif ((node.spring_DY[1] == '-' and node.DY[combo.name] < 0)
-                        or (node.spring_DY[1] == '+' and node.DY[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_DY[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_DY[2] = True
-                    if node.spring_DZ[1] is not None:
-                        if ((node.spring_DZ[1] == '-' and node.DZ[combo.name] > 0)
-                        or (node.spring_DZ[1] == '+' and node.DZ[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_DZ[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_DZ[2] = False
-                        elif ((node.spring_DZ[1] == '-' and node.DZ[combo.name] < 0)
-                        or (node.spring_DZ[1] == '+' and node.DZ[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_DZ[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_DZ[2] = True
-                    if node.spring_RX[1] is not None:
-                        if ((node.spring_RX[1] == '-' and node.RX[combo.name] > 0)
-                        or (node.spring_RX[1] == '+' and node.RX[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_RX[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_RX[2] = False
-                        elif ((node.spring_RX[1] == '-' and node.RX[combo.name] < 0)
-                        or (node.spring_RX[1] == '+' and node.RX[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_RX[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_RX[2] = True
-                    if node.spring_RY[1] is not None:
-                        if ((node.spring_RY[1] == '-' and node.RY[combo.name] > 0)
-                        or (node.spring_RY[1] == '+' and node.RY[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_RY[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_RY[2] = False
-                        elif ((node.spring_RY[1] == '-' and node.RY[combo.name] < 0)
-                        or (node.spring_RY[1] == '+' and node.RY[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_RY[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_RY[2] = True
-                    if node.spring_RZ[1] is not None:
-                        if ((node.spring_RZ[1] == '-' and node.RZ[combo.name] > 0)
-                        or (node.spring_RZ[1] == '+' and node.RZ[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_RZ[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_RZ[2] = False
-                        elif ((node.spring_RZ[1] == '-' and node.RZ[combo.name] < 0)
-                        or (node.spring_RZ[1] == '+' and node.RZ[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_RZ[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_RZ[2] = True
-
-                # Check tension/compression-only springs
-                if log: print('- Checking for tension/compression-only spring convergence')
-                for spring in self.Springs.values():
-
-                    if spring.active[combo.name] == True:
-
-                        # Check if tension-only conditions exist
-                        if spring.tension_only == True and spring.axial(combo.name) > 0:
-                            spring.active[combo.name] = False
-                            convergence = False
-                        
-                        # Check if compression-only conditions exist
-                        elif spring.comp_only == True and spring.axial(combo.name) < 0:
-                            spring.active[combo.name] = False
-                            convergence = False
-
-                # Check tension/compression only members
-                if log: print('- Checking for tension/compression-only member convergence')
-                for phys_member in self.Members.values():
-
-                    # Only run the tension/compression only check if the member is still active
-                    if phys_member.active[combo.name] == True:
-
-                        # Check if tension-only conditions exist
-                        if phys_member.tension_only == True and phys_member.max_axial(combo.name) > 0:
-                            phys_member.active[combo.name] = False
-                            convergence = False
-
-                        # Check if compression-only conditions exist
-                        elif phys_member.comp_only == True and phys_member.min_axial(combo.name) < 0:
-                            phys_member.active[combo.name] = False
-                            convergence = False
-                
                 if convergence == False:
                     if log: print('- Tension/compression-only analysis did not converge. Adjusting stiffness matrix and reanalyzing.')
                 else:
@@ -2553,140 +2437,14 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
                     node.RZ[combo.name] = D[node.ID*6 + 5, 0]
                 
                 # Assume the model has converged (to be checked below)
-                convergence_TC = True
-                
-                # Check tension/compression-only spring supports
-                if log: print('- Checking for tension/compression-only support spring convergence')
-                for node in self.Nodes.values():
-                    
-                    # Check convergence of tension/compression-only spring supports and activate/deactivate them as necessary
-                    if node.spring_DX[1] is not None:
-                        if ((node.spring_DX[1] == '-' and node.DX[combo.name] > 0)
-                        or (node.spring_DX[1] == '+' and node.DX[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_DX[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_DX[2] = False
-                        elif ((node.spring_DX[1] == '-' and node.DX[combo.name] < 0)
-                        or (node.spring_DX[1] == '+' and node.DX[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_DX[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_DX[2] = True
-                    if node.spring_DY[1] is not None:
-                        if ((node.spring_DY[1] == '-' and node.DY[combo.name] > 0)
-                        or (node.spring_DY[1] == '+' and node.DY[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_DY[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_DY[2] = False
-                        elif ((node.spring_DY[1] == '-' and node.DY[combo.name] < 0)
-                        or (node.spring_DY[1] == '+' and node.DY[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_DY[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_DY[2] = True
-                    if node.spring_DZ[1] is not None:
-                        if ((node.spring_DZ[1] == '-' and node.DZ[combo.name] > 0)
-                        or (node.spring_DZ[1] == '+' and node.DZ[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_DZ[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_DZ[2] = False
-                        elif ((node.spring_DZ[1] == '-' and node.DZ[combo.name] < 0)
-                        or (node.spring_DZ[1] == '+' and node.DZ[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_DZ[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_DZ[2] = True
-                    if node.spring_RX[1] is not None:
-                        if ((node.spring_RX[1] == '-' and node.RX[combo.name] > 0)
-                        or (node.spring_RX[1] == '+' and node.RX[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_RX[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_RX[2] = False
-                        elif ((node.spring_RX[1] == '-' and node.RX[combo.name] < 0)
-                        or (node.spring_RX[1] == '+' and node.RX[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_RX[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_RX[2] = True
-                    if node.spring_RY[1] is not None:
-                        if ((node.spring_RY[1] == '-' and node.RY[combo.name] > 0)
-                        or (node.spring_RY[1] == '+' and node.RY[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_RY[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_RY[2] = False
-                        elif ((node.spring_RY[1] == '-' and node.RY[combo.name] < 0)
-                        or (node.spring_RY[1] == '+' and node.RY[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_RY[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_RY[2] = True
-                    if node.spring_RZ[1] is not None:
-                        if ((node.spring_RZ[1] == '-' and node.RZ[combo.name] > 0)
-                        or (node.spring_RZ[1] == '+' and node.RZ[combo.name] < 0)):
-                            # Check if the spring is switching from active to inactive
-                            if node.spring_RZ[2] == True: convergence = False
-                            # Make sure the spring is innactive
-                            node.spring_RZ[2] = False
-                        elif ((node.spring_RZ[1] == '-' and node.RZ[combo.name] < 0)
-                        or (node.spring_RZ[1] == '+' and node.RZ[combo.name] > 0)):
-                            # Check if the spring is switching from inactive to active
-                            if node.spring_RZ[2] == False: convergence = False
-                            # Make sure the spring is active
-                            node.spring_RZ[2] = True
-
-                # Check for tension/compression-only springs that need to be deactivated
-                if log: print('- Checking for tension/compression-only spring convergence')
-                for spring in self.Springs.values():
-
-                    # Only run the tension/compression only check if the spring is still active
-                    if spring.active[combo.name] == True:
-
-                        # Check if tension-only conditions exist
-                        if spring.tension_only == True and spring.axial(combo.name) > 0:
-                            
-                            spring.active[combo.name] = False
-                            convergence_TC = False
-
-                        # Check if compression-only conditions exist
-                        elif spring.comp_only == True and spring.axial(combo.name) < 0:
-                            
-                            spring.active[combo.name] = False
-                            convergence_TC = False
-                
-                # Check for tension/compression-only members that need to be deactivated
-                if log: print('- Checking for tension/compression-only member convergence')
-                for phys_member in self.Members.values():
-
-                    # Only run the tension/compression only check if the member is still active
-                    if phys_member.active[combo.name] == True:
-
-                        # Check if tension-only conditions exist
-                        if phys_member.tension_only == True and phys_member.max_axial(combo.name) > 0:
-                            
-                            phys_member.active[combo.name] = False
-                            for member in phys_member.sub_members.values():
-                                member.active[combo.name] = False
-                            convergence_TC = False
-
-                        # Check if compression-only conditions exist
-                        elif phys_member.comp_only == True and phys_member.min_axial(combo.name) < 0:
-                            
-                            phys_member.active[combo.name] = False
-                            for member in phys_member.sub_members.values():
-                                member.active[combo.name] = False
-                            convergence_TC = False
+                convergence_TC = Analysis._check_TC_convergence(self, combo.name, log)
                 
                 # Report on convergence of tension/compression only analysis
                 if convergence_TC == False:
                     
                     if log:
                         print('- Tension/compression-only analysis did not converge on this iteration')
-                        print('- Stiffness matrix will be adjusted for newly deactivated elements')
+                        print('- Stiffness matrix will be adjusted')
                         print('- P-Delta analysis will be restarted')
                     
                     # Increment the tension/compression-only iteration count

From 8378bf1aac18526c2a10732f85c3898dada29278 Mon Sep 17 00:00:00 2001
From: Craig <33323539+JWock82@users.noreply.github.com>
Date: Sat, 22 Jul 2023 13:44:52 -0400
Subject: [PATCH 06/20] Revert "Allow Selective Analysis Of Load Combos"

---
 PyNite/FEModel3D.py | 107 +++++++++++---------------------------------
 1 file changed, 27 insertions(+), 80 deletions(-)

diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index 0fc11639..3c854eda 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -1903,7 +1903,7 @@ def _partition(self, unp_matrix, D1_indices, D2_indices):
             m22 = unp_matrix[D2_indices, :][:, D2_indices]
             return m11, m12, m21, m22
 
-    def analyze(self, log=False, check_stability=True, check_statics=False, max_iter=30, sparse=True,combos:list=None):
+    def analyze(self, log=False, check_stability=True, check_statics=False, max_iter=30, sparse=True):
         """Performs first-order static analysis. Iterations are performed if tension-only members or compression-only members are present.
 
         :param log: Prints the analysis log to the console if set to True. Default is False.
@@ -1918,8 +1918,6 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
         :type sparse: bool, optional
         :raises Exception: _description_
         :raises Exception: _description_
-        :param combos: provide a subset of combos for analysis that will be run. If None or no input is provided all combos will be run.
-        :type combos: list, optional 
         """
 
         if log:
@@ -1944,20 +1942,12 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
         # Activate all springs and members for all load combinations
         for spring in self.Springs.values():
             for combo_name in self.LoadCombos.keys():
-                #Check if specific combos and if this combo is selected
-                if combos is not None and combo_name in combos:
-                    spring.active[combo_name] = True
-                elif combos is None:
-                    spring.active[combo_name] = True                
+                spring.active[combo_name] = True
         
         # Activate all physical members for all load combinations
         for phys_member in self.Members.values():
             for combo_name in self.LoadCombos.keys():
-                #Check if specific combos and if this combo is selected
-                if combos is not None and combo_name in combos:
-                    phys_member.active[combo_name] = True
-                elif combos is None:
-                    phys_member.active[combo_name] = True
+                phys_member.active[combo_name] = True
         
         # Assign an internal ID to all nodes and elements in the model
         self._renumber()
@@ -1969,12 +1959,7 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
         D2 = atleast_2d(D2).T
 
         # Step through each load combination
-        for combo_name,combo in self.LoadCombos.items():
-
-            #Check if specific combos and if this combo is selected
-            if combos is not None and combo_name not in combos:
-                continue     
-
+        for combo in self.LoadCombos.values():
 
             if log:
                 print('')
@@ -2085,7 +2070,7 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
                 iter_count += 1
 
         # Calculate reactions
-        self._calc_reactions(combos=combos)
+        self._calc_reactions()
 
         if log:
             print('')     
@@ -2094,12 +2079,12 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
 
         # Check statics if requested
         if check_statics == True:
-            self._check_statics(combos=combos)
+            self._check_statics()
         
         # Flag the model as solved
         self.solution = 'Linear TC'
 
-    def analyze_linear(self, log=False, check_stability=True, check_statics=False, sparse=True, combos=None):
+    def analyze_linear(self, log=False, check_stability=True, check_statics=False, sparse=True):
         """Performs first-order static analysis. This analysis procedure is much faster since it only assembles the global stiffness matrix once, rather than once for each load combination. It is not appropriate when non-linear behavior such as tension/compression only analysis or P-Delta analysis are required.
 
         :param log: Prints the analysis log to the console if set to True. Default is False.
@@ -2111,8 +2096,6 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
         :param sparse: Indicates whether the sparse matrix solver should be used. A matrix can be considered sparse or dense depening on how many zero terms there are. Structural stiffness matrices often contain many zero terms. The sparse solver can offer faster solutions for such matrices. Using the sparse solver on dense matrices may lead to slower solution times. Be sure ``scipy`` is installed to use the sparse solver. Default is True.
         :type sparse: bool, optional
         :raises Exception: Occurs when a singular stiffness matrix is found. This indicates an unstable structure has been modeled.
-        :param combos: provide a subset of combos for analysis that will be run. If None or no input is provided all combos will be run.
-        :type combos: list, optional 
         """
 
         if log:
@@ -2134,23 +2117,15 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
             if mesh.is_generated == False:
                 mesh.generate()
         
-        # Activate all springs and members for all load combinations
+        # Activate all springs for all load combinations
         for spring in self.Springs.values():
             for combo_name in self.LoadCombos.keys():
-                #Check if specific combos and if this combo is selected
-                if combos is not None and combo_name in combos:
-                    spring.active[combo_name] = True
-                elif combos is None:
-                    spring.active[combo_name] = True                
+                spring.active[combo_name] = True
         
         # Activate all physical members for all load combinations
         for phys_member in self.Members.values():
             for combo_name in self.LoadCombos.keys():
-                #Check if specific combos and if this combo is selected
-                if combos is not None and combo_name in combos:
-                    phys_member.active[combo_name] = True
-                elif combos is None:
-                    phys_member.active[combo_name] = True
+                phys_member.active[combo_name] = True
         
         # Assign an internal ID to all nodes and elements in the model
         self._renumber()
@@ -2169,11 +2144,7 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
             K11, K12, K21, K22 = self._partition(self.K(combo_name, log, check_stability, sparse), D1_indices, D2_indices)
 
         # Step through each load combination
-        for combo_name,combo in self.LoadCombos.items():
-
-            #Check if specific combos and if this combo is selected
-            if combos is not None and combo_name not in combos:
-                continue            
+        for combo in self.LoadCombos.values():
 
             if log:
                 print('')
@@ -2254,7 +2225,7 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
                 node.RZ[combo.name] = D[node.ID*6 + 5, 0]
 
         # Calculate reactions
-        self._calc_reactions(combos=combos)
+        self._calc_reactions()
 
         if log:
             print('')     
@@ -2263,12 +2234,12 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
 
         # Check statics if requested
         if check_statics == True:
-            self._check_statics(combos=combos)
+            self._check_statics()
         
         # Flag the model as solved
         self.solution = 'Linear'
 
-    def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01, sparse=True,combos=None):
+    def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01, sparse=True):
         """Performs second order (P-Delta) analysis. This type of analysis is appropriate for most models using beams, columns and braces. Second order analysis is usually required by material specific codes. The analysis is iterative and takes longer to solve. Models with slender members and/or members with combined bending and axial loads will generally have more significant P-Delta effects. P-Delta effects in plates/quads are not considered.
 
         :param log: Prints updates to the console if set to True. Default is False.
@@ -2283,8 +2254,6 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         :type sparse: bool, optional
         :raises ValueError: Occurs when there is a singularity in the stiffness matrix, which indicates an unstable structure.
         :raises Exception: Occurs when a model fails to converge.
-        :param combos: provide a subset of combos for analysis that will be run. If None or no input is provided all combos will be run.
-        :type combos: list, optional 
         """
         
         if log:
@@ -2306,23 +2275,16 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
             if mesh.is_generated == False:
                 mesh.generate()
         
-        # Activate all springs and members for all load combinations
+        # Activate all springs for all load combinations. They can be turned inactive
+        # during the course of the tension/compression-only analysis
         for spring in self.Springs.values():
             for combo_name in self.LoadCombos.keys():
-                #Check if specific combos and if this combo is selected
-                if combos is not None and combo_name in combos:
-                    spring.active[combo_name] = True
-                elif combos is None:
-                    spring.active[combo_name] = True                
-        
+                spring.active[combo_name] = True
+                
         # Activate all physical members for all load combinations
         for phys_member in self.Members.values():
             for combo_name in self.LoadCombos.keys():
-                #Check if specific combos and if this combo is selected
-                if combos is not None and combo_name in combos:
-                    phys_member.active[combo_name] = True
-                elif combos is None:
-                    phys_member.active[combo_name] = True
+                phys_member.active[combo_name] = True
                
         # Assign an internal ID to all nodes and elements in the model
         self._renumber()
@@ -2334,11 +2296,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         D2 = array(D2, ndmin=2).T
 
         # Step through each load combination
-        for combo_name,combo in self.LoadCombos.items():
-
-            #Check if specific combos and if this combo is selected
-            if combos is not None and combo_name not in combos:
-                continue   
+        for combo in self.LoadCombos.values():
             
             if log:
                 print('')
@@ -2535,7 +2493,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
                 iter_count_PD += 1
         
         # Calculate reactions
-        self._calc_reactions(combos=combos)
+        self._calc_reactions()
 
         if log:
             print('')
@@ -2545,7 +2503,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         # Flag the model as solved
         self.solution = 'P-Delta'
 
-    def _calc_reactions(self, log=False,combos=None):
+    def _calc_reactions(self, log=False):
         """
         Calculates reactions internally once the model is solved.
 
@@ -2553,8 +2511,6 @@ def _calc_reactions(self, log=False,combos=None):
         ----------
         log : bool, optional
             Prints updates to the console if set to True. Default is False.
-        combos: list, optional
-            Selects which combos to evaluate
         """
 
         # Print a status update to the console
@@ -2564,12 +2520,7 @@ def _calc_reactions(self, log=False,combos=None):
         for node in self.Nodes.values():
             
             # Step through each load combination
-            for combo_name,combo in self.LoadCombos.items():
-
-                #Check if specific combos and if this combo is selected
-                if combos is not None and combo_name not in combos:
-                    continue
-
+            for combo in self.LoadCombos.values():
                 
                 # Initialize reactions for this node and load combination
                 node.RxnFX[combo.name] = 0.0
@@ -2865,14 +2816,14 @@ def _check_stability(self, K):
 
         return
     
-    def _check_statics(self,combos=None):
+    def _check_statics(self):
         '''
         Checks static equilibrium and prints results to the console.
 
         Parameters
         ----------
-        combos: list, optional
-            Selects which combos to evaluate
+        precision : number
+            The number of decimal places to carry the results to.
         '''
 
         print('+----------------+')
@@ -2887,11 +2838,7 @@ def _check_statics(self,combos=None):
         statics_table.field_names = ['Load Combination', 'Sum FX', 'Sum RX', 'Sum FY', 'Sum RY', 'Sum FZ', 'Sum RZ', 'Sum MX', 'Sum RMX', 'Sum MY', 'Sum RMY', 'Sum MZ', 'Sum RMZ']
 
         # Step through each load combination
-        for combo_name,combo in self.LoadCombos.items():
-
-            #Check if specific combos and if this combo is selected
-            if combos is not None and combo_name not in combos:
-                continue
+        for combo in self.LoadCombos.values():
 
             # Initialize force and moment summations to zero
             SumFX, SumFY, SumFZ = 0.0, 0.0, 0.0

From ef6ab49e8cab160e68e5fd0439c1b7d62a02be16 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 22 Jul 2023 13:45:50 -0400
Subject: [PATCH 07/20] Moved analysis methods out of `FEModel3D`

---
 PyNite/Analysis.py  | 402 +++++++++++++++++++++++++++++++++++++++++-
 PyNite/FEModel3D.py | 417 +-------------------------------------------
 2 files changed, 410 insertions(+), 409 deletions(-)

diff --git a/PyNite/Analysis.py b/PyNite/Analysis.py
index d70845f1..548ec358 100644
--- a/PyNite/Analysis.py
+++ b/PyNite/Analysis.py
@@ -1,3 +1,57 @@
+from math import isclose
+
+def _check_stability(model, K):
+    """
+    Identifies nodal instabilities in a model's stiffness matrix.
+    """
+
+    # Initialize the `unstable` flag to `False`
+    unstable = False
+
+    # Step through each diagonal term in the stiffness matrix
+    for i in range(K.shape[0]):
+        
+        # Determine which node this term belongs to
+        node = [node for node in model.Nodes.values() if node.ID == int(i/6)][0]
+
+        # Determine which degree of freedom this term belongs to
+        dof = i%6
+
+        # Check to see if this degree of freedom is supported
+        if dof == 0:
+            supported = node.support_DX
+        elif dof == 1:
+            supported = node.support_DY
+        elif dof == 2:
+            supported = node.support_DZ
+        elif dof == 3:
+            supported = node.support_RX
+        elif dof == 4:
+            supported = node.support_RY
+        elif dof == 5:
+            supported = node.support_RZ
+
+        # Check if the degree of freedom on this diagonal is unstable
+        if isclose(K[i, i], 0) and not supported:
+
+            # Flag the model as unstable
+            unstable = True
+
+            # Identify which direction this instability effects
+            if i%6 == 0: direction = 'for translation in the global X direction.'
+            if i%6 == 1: direction = 'for translation in the global Y direction.'
+            if i%6 == 2: direction = 'for translation in the global Z direction.'
+            if i%6 == 3: direction = 'for rotation about the global X axis.'
+            if i%6 == 4: direction = 'for rotation about the global Y axis.'
+            if i%6 == 5: direction = 'for rotation about the global Z axis.'
+
+            # Print a message to the console
+            print('* Nodal instability detected: node ' + node.name + ' is unstable ' + direction)
+
+    if unstable:
+        raise Exception('Unstable node(s). See console output for details.')
+
+    return
 
 def _check_TC_convergence(model, combo_name='Combo 1', log=True):
     
@@ -124,4 +178,350 @@ def _check_TC_convergence(model, combo_name='Combo 1', log=True):
                 convergence = False
 
     # Return whether the TC analysis has converged
-    return convergence
\ No newline at end of file
+    return convergence
+
+def _calc_reactions(model, log=False):
+    """
+    Calculates reactions internally once the model is solved.
+
+    Parameters
+    ----------
+    log : bool, optional
+        Prints updates to the console if set to True. Default is False.
+    """
+
+    # Print a status update to the console
+    if log: print('- Calculating reactions')
+
+    # Calculate the reactions node by node
+    for node in model.Nodes.values():
+        
+        # Step through each load combination
+        for combo in model.LoadCombos.values():
+            
+            # Initialize reactions for this node and load combination
+            node.RxnFX[combo.name] = 0.0
+            node.RxnFY[combo.name] = 0.0
+            node.RxnFZ[combo.name] = 0.0
+            node.RxnMX[combo.name] = 0.0
+            node.RxnMY[combo.name] = 0.0
+            node.RxnMZ[combo.name] = 0.0
+
+            # Determine if the node has any supports
+            if (node.support_DX or node.support_DY or node.support_DZ 
+            or  node.support_RX or node.support_RY or node.support_RZ):
+
+                # Sum the spring end forces at the node
+                for spring in model.Springs.values():
+
+                    if spring.i_node == node and spring.active[combo.name] == True:
+                        
+                        # Get the spring's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        spring_F = spring.F(combo.name)
+
+                        node.RxnFX[combo.name] += spring_F[0, 0]
+                        node.RxnFY[combo.name] += spring_F[1, 0]
+                        node.RxnFZ[combo.name] += spring_F[2, 0]
+                        node.RxnMX[combo.name] += spring_F[3, 0]
+                        node.RxnMY[combo.name] += spring_F[4, 0]
+                        node.RxnMZ[combo.name] += spring_F[5, 0]
+
+                    elif spring.j_node == node and spring.active[combo.name] == True:
+                    
+                        # Get the spring's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        spring_F = spring.F(combo.name)
+                    
+                        node.RxnFX[combo.name] += spring_F[6, 0]
+                        node.RxnFY[combo.name] += spring_F[7, 0]
+                        node.RxnFZ[combo.name] += spring_F[8, 0]
+                        node.RxnMX[combo.name] += spring_F[9, 0]
+                        node.RxnMY[combo.name] += spring_F[10, 0]
+                        node.RxnMZ[combo.name] += spring_F[11, 0]
+
+                # Step through each physical member in the model
+                for phys_member in model.Members.values():
+
+                    # Sum the sub-member end forces at the node
+                    for member in phys_member.sub_members.values():
+                        
+                        if member.i_node == node and phys_member.active[combo.name] == True:
+                        
+                            # Get the member's global force matrix
+                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                            member_F = member.F(combo.name)
+
+                            node.RxnFX[combo.name] += member_F[0, 0]
+                            node.RxnFY[combo.name] += member_F[1, 0]
+                            node.RxnFZ[combo.name] += member_F[2, 0]
+                            node.RxnMX[combo.name] += member_F[3, 0]
+                            node.RxnMY[combo.name] += member_F[4, 0]
+                            node.RxnMZ[combo.name] += member_F[5, 0]
+
+                        elif member.j_node == node and phys_member.active[combo.name] == True:
+                        
+                            # Get the member's global force matrix
+                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                            member_F = member.F(combo.name)
+                        
+                            node.RxnFX[combo.name] += member_F[6, 0]
+                            node.RxnFY[combo.name] += member_F[7, 0]
+                            node.RxnFZ[combo.name] += member_F[8, 0]
+                            node.RxnMX[combo.name] += member_F[9, 0]
+                            node.RxnMY[combo.name] += member_F[10, 0]
+                            node.RxnMZ[combo.name] += member_F[11, 0]
+
+                # Sum the plate forces at the node
+                for plate in model.Plates.values():
+
+                    if plate.i_node == node:
+
+                        # Get the plate's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        plate_F = plate.F(combo.name)
+                
+                        node.RxnFX[combo.name] += plate_F[0, 0]
+                        node.RxnFY[combo.name] += plate_F[1, 0]
+                        node.RxnFZ[combo.name] += plate_F[2, 0]
+                        node.RxnMX[combo.name] += plate_F[3, 0]
+                        node.RxnMY[combo.name] += plate_F[4, 0]
+                        node.RxnMZ[combo.name] += plate_F[5, 0]
+
+                    elif plate.j_node == node:
+
+                        # Get the plate's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        plate_F = plate.F(combo.name)
+                
+                        node.RxnFX[combo.name] += plate_F[6, 0]
+                        node.RxnFY[combo.name] += plate_F[7, 0]
+                        node.RxnFZ[combo.name] += plate_F[8, 0]
+                        node.RxnMX[combo.name] += plate_F[9, 0]
+                        node.RxnMY[combo.name] += plate_F[10, 0]
+                        node.RxnMZ[combo.name] += plate_F[11, 0]
+
+                    elif plate.m_node == node:
+
+                        # Get the plate's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        plate_F = plate.F(combo.name)
+                
+                        node.RxnFX[combo.name] += plate_F[12, 0]
+                        node.RxnFY[combo.name] += plate_F[13, 0]
+                        node.RxnFZ[combo.name] += plate_F[14, 0]
+                        node.RxnMX[combo.name] += plate_F[15, 0]
+                        node.RxnMY[combo.name] += plate_F[16, 0]
+                        node.RxnMZ[combo.name] += plate_F[17, 0]
+
+                    elif plate.n_node == node:
+
+                        # Get the plate's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        plate_F = plate.F(combo.name)
+                
+                        node.RxnFX[combo.name] += plate_F[18, 0]
+                        node.RxnFY[combo.name] += plate_F[19, 0]
+                        node.RxnFZ[combo.name] += plate_F[20, 0]
+                        node.RxnMX[combo.name] += plate_F[21, 0]
+                        node.RxnMY[combo.name] += plate_F[22, 0]
+                        node.RxnMZ[combo.name] += plate_F[23, 0]
+
+                # Sum the quad forces at the node
+                for quad in model.Quads.values():
+
+                    if quad.m_node == node:
+
+                        # Get the quad's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        quad_F = quad.F(combo.name)
+
+                        node.RxnFX[combo.name] += quad_F[0, 0]
+                        node.RxnFY[combo.name] += quad_F[1, 0]
+                        node.RxnFZ[combo.name] += quad_F[2, 0]
+                        node.RxnMX[combo.name] += quad_F[3, 0]
+                        node.RxnMY[combo.name] += quad_F[4, 0]
+                        node.RxnMZ[combo.name] += quad_F[5, 0]
+
+                    elif quad.n_node == node:
+
+                        # Get the quad's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        quad_F = quad.F(combo.name)
+                
+                        node.RxnFX[combo.name] += quad_F[6, 0]
+                        node.RxnFY[combo.name] += quad_F[7, 0]
+                        node.RxnFZ[combo.name] += quad_F[8, 0]
+                        node.RxnMX[combo.name] += quad_F[9, 0]
+                        node.RxnMY[combo.name] += quad_F[10, 0]
+                        node.RxnMZ[combo.name] += quad_F[11, 0]
+
+                    elif quad.i_node == node:
+
+                        # Get the quad's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        quad_F = quad.F(combo.name)
+                
+                        node.RxnFX[combo.name] += quad_F[12, 0]
+                        node.RxnFY[combo.name] += quad_F[13, 0]
+                        node.RxnFZ[combo.name] += quad_F[14, 0]
+                        node.RxnMX[combo.name] += quad_F[15, 0]
+                        node.RxnMY[combo.name] += quad_F[16, 0]
+                        node.RxnMZ[combo.name] += quad_F[17, 0]
+
+                    elif quad.j_node == node:
+
+                        # Get the quad's global force matrix
+                        # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
+                        quad_F = quad.F(combo.name)
+                
+                        node.RxnFX[combo.name] += quad_F[18, 0]
+                        node.RxnFY[combo.name] += quad_F[19, 0]
+                        node.RxnFZ[combo.name] += quad_F[20, 0]
+                        node.RxnMX[combo.name] += quad_F[21, 0]
+                        node.RxnMY[combo.name] += quad_F[22, 0]
+                        node.RxnMZ[combo.name] += quad_F[23, 0]
+                
+                # Sum the joint loads applied to the node
+                for load in node.NodeLoads:
+
+                    for case, factor in combo.factors.items():
+                        
+                        if load[2] == case:
+
+                            if load[0] == 'FX':
+                                node.RxnFX[combo.name] -= load[1]*factor
+                            elif load[0] == 'FY':
+                                node.RxnFY[combo.name] -= load[1]*factor
+                            elif load[0] == 'FZ':
+                                node.RxnFZ[combo.name] -= load[1]*factor
+                            elif load[0] == 'MX':
+                                node.RxnMX[combo.name] -= load[1]*factor
+                            elif load[0] == 'MY':
+                                node.RxnMY[combo.name] -= load[1]*factor
+                            elif load[0] == 'MZ':
+                                node.RxnMZ[combo.name] -= load[1]*factor
+            
+            # Calculate reactions due to active spring supports at the node
+            elif node.spring_DX[0] != None and node.spring_DX[2] == True:
+                sign = node.spring_DX[1]
+                k = node.spring_DX[0]
+                if sign != None: k = float(sign + str(k))
+                DX = node.DX[combo.name]
+                node.RxnFX[combo.name] += k*DX
+            elif node.spring_DY[0] != None and node.spring_DY[2] == True:
+                sign = node.spring_DY[1]
+                k = node.spring_DY[0]
+                if sign != None: k = float(sign + str(k))
+                DY = node.DY[combo.name]
+                node.RxnFY[combo.name] += k*DY
+            elif node.spring_DZ[0] != None and node.spring_DZ[2] == True:
+                sign = node.spring_DZ[1]
+                k = node.spring_DZ[0]
+                if sign != None: k = float(sign + str(k))
+                DZ = node.DZ[combo.name]
+                node.RxnFZ[combo.name] += k*DZ
+            elif node.spring_RX[0] != None and node.spring_RX[2] == True:
+                sign = node.spring_RX[1]
+                k = node.spring_RX[0]
+                if sign != None: k = float(sign + str(k))
+                RX = node.RX[combo.name]
+                node.RxnMX[combo.name] += k*RX
+            elif node.spring_RY[0] != None and node.spring_RY[2] == True:
+                sign = node.spring_RY[1]
+                k = node.spring_RY[0]
+                if sign != None: k = float(sign + str(k))
+                RY = node.RY[combo.name]
+                node.RxnMY[combo.name] += k*RY
+            elif node.spring_RZ[0] != None and node.spring_RZ[2] == True:
+                sign = node.spring_RZ[1]
+                k = node.spring_RZ[0]
+                if sign != None: k = float(sign + str(k))
+                RZ = node.RZ[combo.name]
+                node.RxnMZ[combo.name] += k*RZ
+
+def _check_statics(model):
+    '''
+    Checks static equilibrium and prints results to the console.
+
+    Parameters
+    ----------
+    precision : number
+        The number of decimal places to carry the results to.
+    '''
+
+    print('+----------------+')
+    print('| Statics Check: |')
+    print('+----------------+')
+    print('')
+
+    from prettytable import PrettyTable
+
+    # Start a blank table and create a header row
+    statics_table = PrettyTable()
+    statics_table.field_names = ['Load Combination', 'Sum FX', 'Sum RX', 'Sum FY', 'Sum RY', 'Sum FZ', 'Sum RZ', 'Sum MX', 'Sum RMX', 'Sum MY', 'Sum RMY', 'Sum MZ', 'Sum RMZ']
+
+    # Step through each load combination
+    for combo in model.LoadCombos.values():
+
+        # Initialize force and moment summations to zero
+        SumFX, SumFY, SumFZ = 0.0, 0.0, 0.0
+        SumMX, SumMY, SumMZ = 0.0, 0.0, 0.0
+        SumRFX, SumRFY, SumRFZ = 0.0, 0.0, 0.0
+        SumRMX, SumRMY, SumRMZ = 0.0, 0.0, 0.0
+
+        # Get the global force vector and the global fixed end reaction vector
+        P = model.P(combo.name)
+        FER = model.FER(combo.name)
+
+        # Step through each node and sum its forces
+        for node in model.Nodes.values():
+
+            # Get the node's coordinates
+            X = node.X
+            Y = node.Y
+            Z = node.Z
+
+            # Get the nodal forces
+            FX = P[node.ID*6+0][0] - FER[node.ID*6+0][0]
+            FY = P[node.ID*6+1][0] - FER[node.ID*6+1][0]
+            FZ = P[node.ID*6+2][0] - FER[node.ID*6+2][0]
+            MX = P[node.ID*6+3][0] - FER[node.ID*6+3][0]
+            MY = P[node.ID*6+4][0] - FER[node.ID*6+4][0]
+            MZ = P[node.ID*6+5][0] - FER[node.ID*6+5][0]
+
+            # Get the nodal reactions
+            RFX = node.RxnFX[combo.name]
+            RFY = node.RxnFY[combo.name]
+            RFZ = node.RxnFZ[combo.name]
+            RMX = node.RxnMX[combo.name]
+            RMY = node.RxnMY[combo.name]
+            RMZ = node.RxnMZ[combo.name]
+
+            # Sum the global forces
+            SumFX += FX
+            SumFY += FY
+            SumFZ += FZ
+            SumMX += MX - FY*Z + FZ*Y
+            SumMY += MY + FX*Z - FZ*X
+            SumMZ += MZ - FX*Y + FY*X
+
+            # Sum the global reactions
+            SumRFX += RFX
+            SumRFY += RFY
+            SumRFZ += RFZ
+            SumRMX += RMX - RFY*Z + RFZ*Y
+            SumRMY += RMY + RFX*Z - RFZ*X
+            SumRMZ += RMZ - RFX*Y + RFY*X   
+
+        # Add the results to the table
+        statics_table.add_row([combo.name, '{:.3g}'.format(SumFX), '{:.3g}'.format(SumRFX),
+                                            '{:.3g}'.format(SumFY), '{:.3g}'.format(SumRFY),
+                                            '{:.3g}'.format(SumFZ), '{:.3g}'.format(SumRFZ),
+                                            '{:.3g}'.format(SumMX), '{:.3g}'.format(SumRMX),
+                                            '{:.3g}'.format(SumMY), '{:.3g}'.format(SumRMY),
+                                            '{:.3g}'.format(SumMZ), '{:.3g}'.format(SumRMZ)])
+
+    # Print the static check table
+    print(statics_table)
+    print('')
\ No newline at end of file
diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index 3c854eda..b9267515 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -1631,8 +1631,8 @@ def K(self, combo_name='Combo 1', log=False, check_stability=True, sparse=True):
         # Check that there are no nodal instabilities
         if check_stability:
             if log: print('- Checking nodal stability')
-            if sparse: self._check_stability(K.tocsr())
-            else: self._check_stability(K)
+            if sparse: Analysis._check_stability(self, K.tocsr())
+            else: Analysis._check_stability(self, K)
 
         # Return the global stiffness matrix
         return K    
@@ -2070,7 +2070,7 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
                 iter_count += 1
 
         # Calculate reactions
-        self._calc_reactions()
+        Analysis._calc_reactions(self)
 
         if log:
             print('')     
@@ -2079,7 +2079,7 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
 
         # Check statics if requested
         if check_statics == True:
-            self._check_statics()
+            Analysis._check_statics(self)
         
         # Flag the model as solved
         self.solution = 'Linear TC'
@@ -2225,7 +2225,7 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
                 node.RZ[combo.name] = D[node.ID*6 + 5, 0]
 
         # Calculate reactions
-        self._calc_reactions()
+        Analysis._calc_reactions(self)
 
         if log:
             print('')     
@@ -2234,7 +2234,7 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
 
         # Check statics if requested
         if check_statics == True:
-            self._check_statics()
+            Analysis._check_statics(self)
         
         # Flag the model as solved
         self.solution = 'Linear'
@@ -2331,7 +2331,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
                                                        
                     # Check that the structure is stable
                     if log: print('- Checking stability')
-                    self._check_stability(K11)
+                    Analysis._check_stability(self, K11)
 
                     # Assemble the force matrices
                     FER1, FER2 = self._partition(self.FER(combo.name), D1_indices, D2_indices)  # Fixed end reactions
@@ -2493,7 +2493,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
                 iter_count_PD += 1
         
         # Calculate reactions
-        self._calc_reactions()
+        Analysis._calc_reactions(self)
 
         if log:
             print('')
@@ -2502,406 +2502,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         
         # Flag the model as solved
         self.solution = 'P-Delta'
-
-    def _calc_reactions(self, log=False):
-        """
-        Calculates reactions internally once the model is solved.
-
-        Parameters
-        ----------
-        log : bool, optional
-            Prints updates to the console if set to True. Default is False.
-        """
-
-        # Print a status update to the console
-        if log: print('- Calculating reactions')
-
-        # Calculate the reactions node by node
-        for node in self.Nodes.values():
-            
-            # Step through each load combination
-            for combo in self.LoadCombos.values():
-                
-                # Initialize reactions for this node and load combination
-                node.RxnFX[combo.name] = 0.0
-                node.RxnFY[combo.name] = 0.0
-                node.RxnFZ[combo.name] = 0.0
-                node.RxnMX[combo.name] = 0.0
-                node.RxnMY[combo.name] = 0.0
-                node.RxnMZ[combo.name] = 0.0
-
-                # Determine if the node has any supports
-                if (node.support_DX or node.support_DY or node.support_DZ 
-                or  node.support_RX or node.support_RY or node.support_RZ):
-
-                    # Sum the spring end forces at the node
-                    for spring in self.Springs.values():
-
-                        if spring.i_node == node and spring.active[combo.name] == True:
-                            
-                            # Get the spring's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            spring_F = spring.F(combo.name)
-
-                            node.RxnFX[combo.name] += spring_F[0, 0]
-                            node.RxnFY[combo.name] += spring_F[1, 0]
-                            node.RxnFZ[combo.name] += spring_F[2, 0]
-                            node.RxnMX[combo.name] += spring_F[3, 0]
-                            node.RxnMY[combo.name] += spring_F[4, 0]
-                            node.RxnMZ[combo.name] += spring_F[5, 0]
-
-                        elif spring.j_node == node and spring.active[combo.name] == True:
-                        
-                            # Get the spring's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            spring_F = spring.F(combo.name)
-                        
-                            node.RxnFX[combo.name] += spring_F[6, 0]
-                            node.RxnFY[combo.name] += spring_F[7, 0]
-                            node.RxnFZ[combo.name] += spring_F[8, 0]
-                            node.RxnMX[combo.name] += spring_F[9, 0]
-                            node.RxnMY[combo.name] += spring_F[10, 0]
-                            node.RxnMZ[combo.name] += spring_F[11, 0]
-
-                    # Step through each physical member in the model
-                    for phys_member in self.Members.values():
-
-                        # Sum the sub-member end forces at the node
-                        for member in phys_member.sub_members.values():
-                            
-                            if member.i_node == node and phys_member.active[combo.name] == True:
-                            
-                                # Get the member's global force matrix
-                                # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                                member_F = member.F(combo.name)
-
-                                node.RxnFX[combo.name] += member_F[0, 0]
-                                node.RxnFY[combo.name] += member_F[1, 0]
-                                node.RxnFZ[combo.name] += member_F[2, 0]
-                                node.RxnMX[combo.name] += member_F[3, 0]
-                                node.RxnMY[combo.name] += member_F[4, 0]
-                                node.RxnMZ[combo.name] += member_F[5, 0]
-
-                            elif member.j_node == node and phys_member.active[combo.name] == True:
-                            
-                                # Get the member's global force matrix
-                                # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                                member_F = member.F(combo.name)
-                            
-                                node.RxnFX[combo.name] += member_F[6, 0]
-                                node.RxnFY[combo.name] += member_F[7, 0]
-                                node.RxnFZ[combo.name] += member_F[8, 0]
-                                node.RxnMX[combo.name] += member_F[9, 0]
-                                node.RxnMY[combo.name] += member_F[10, 0]
-                                node.RxnMZ[combo.name] += member_F[11, 0]
-
-                    # Sum the plate forces at the node
-                    for plate in self.Plates.values():
-
-                        if plate.i_node == node:
-
-                            # Get the plate's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            plate_F = plate.F(combo.name)
-                    
-                            node.RxnFX[combo.name] += plate_F[0, 0]
-                            node.RxnFY[combo.name] += plate_F[1, 0]
-                            node.RxnFZ[combo.name] += plate_F[2, 0]
-                            node.RxnMX[combo.name] += plate_F[3, 0]
-                            node.RxnMY[combo.name] += plate_F[4, 0]
-                            node.RxnMZ[combo.name] += plate_F[5, 0]
-
-                        elif plate.j_node == node:
-
-                            # Get the plate's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            plate_F = plate.F(combo.name)
-                    
-                            node.RxnFX[combo.name] += plate_F[6, 0]
-                            node.RxnFY[combo.name] += plate_F[7, 0]
-                            node.RxnFZ[combo.name] += plate_F[8, 0]
-                            node.RxnMX[combo.name] += plate_F[9, 0]
-                            node.RxnMY[combo.name] += plate_F[10, 0]
-                            node.RxnMZ[combo.name] += plate_F[11, 0]
-
-                        elif plate.m_node == node:
-
-                            # Get the plate's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            plate_F = plate.F(combo.name)
-                    
-                            node.RxnFX[combo.name] += plate_F[12, 0]
-                            node.RxnFY[combo.name] += plate_F[13, 0]
-                            node.RxnFZ[combo.name] += plate_F[14, 0]
-                            node.RxnMX[combo.name] += plate_F[15, 0]
-                            node.RxnMY[combo.name] += plate_F[16, 0]
-                            node.RxnMZ[combo.name] += plate_F[17, 0]
-
-                        elif plate.n_node == node:
-
-                            # Get the plate's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            plate_F = plate.F(combo.name)
-                    
-                            node.RxnFX[combo.name] += plate_F[18, 0]
-                            node.RxnFY[combo.name] += plate_F[19, 0]
-                            node.RxnFZ[combo.name] += plate_F[20, 0]
-                            node.RxnMX[combo.name] += plate_F[21, 0]
-                            node.RxnMY[combo.name] += plate_F[22, 0]
-                            node.RxnMZ[combo.name] += plate_F[23, 0]
-
-                    # Sum the quad forces at the node
-                    for quad in self.Quads.values():
-
-                        if quad.m_node == node:
-
-                            # Get the quad's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            quad_F = quad.F(combo.name)
-
-                            node.RxnFX[combo.name] += quad_F[0, 0]
-                            node.RxnFY[combo.name] += quad_F[1, 0]
-                            node.RxnFZ[combo.name] += quad_F[2, 0]
-                            node.RxnMX[combo.name] += quad_F[3, 0]
-                            node.RxnMY[combo.name] += quad_F[4, 0]
-                            node.RxnMZ[combo.name] += quad_F[5, 0]
-
-                        elif quad.n_node == node:
-
-                            # Get the quad's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            quad_F = quad.F(combo.name)
-                    
-                            node.RxnFX[combo.name] += quad_F[6, 0]
-                            node.RxnFY[combo.name] += quad_F[7, 0]
-                            node.RxnFZ[combo.name] += quad_F[8, 0]
-                            node.RxnMX[combo.name] += quad_F[9, 0]
-                            node.RxnMY[combo.name] += quad_F[10, 0]
-                            node.RxnMZ[combo.name] += quad_F[11, 0]
-
-                        elif quad.i_node == node:
-
-                            # Get the quad's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            quad_F = quad.F(combo.name)
-                    
-                            node.RxnFX[combo.name] += quad_F[12, 0]
-                            node.RxnFY[combo.name] += quad_F[13, 0]
-                            node.RxnFZ[combo.name] += quad_F[14, 0]
-                            node.RxnMX[combo.name] += quad_F[15, 0]
-                            node.RxnMY[combo.name] += quad_F[16, 0]
-                            node.RxnMZ[combo.name] += quad_F[17, 0]
-
-                        elif quad.j_node == node:
-
-                            # Get the quad's global force matrix
-                            # Storing it as a local variable eliminates the need to rebuild it every time a term is needed                    
-                            quad_F = quad.F(combo.name)
-                    
-                            node.RxnFX[combo.name] += quad_F[18, 0]
-                            node.RxnFY[combo.name] += quad_F[19, 0]
-                            node.RxnFZ[combo.name] += quad_F[20, 0]
-                            node.RxnMX[combo.name] += quad_F[21, 0]
-                            node.RxnMY[combo.name] += quad_F[22, 0]
-                            node.RxnMZ[combo.name] += quad_F[23, 0]
-                    
-                    # Sum the joint loads applied to the node
-                    for load in node.NodeLoads:
-
-                        for case, factor in combo.factors.items():
-                            
-                            if load[2] == case:
-
-                                if load[0] == 'FX':
-                                    node.RxnFX[combo.name] -= load[1]*factor
-                                elif load[0] == 'FY':
-                                    node.RxnFY[combo.name] -= load[1]*factor
-                                elif load[0] == 'FZ':
-                                    node.RxnFZ[combo.name] -= load[1]*factor
-                                elif load[0] == 'MX':
-                                    node.RxnMX[combo.name] -= load[1]*factor
-                                elif load[0] == 'MY':
-                                    node.RxnMY[combo.name] -= load[1]*factor
-                                elif load[0] == 'MZ':
-                                    node.RxnMZ[combo.name] -= load[1]*factor
-                
-                # Calculate reactions due to active spring supports at the node
-                elif node.spring_DX[0] != None and node.spring_DX[2] == True:
-                    sign = node.spring_DX[1]
-                    k = node.spring_DX[0]
-                    if sign != None: k = float(sign + str(k))
-                    DX = node.DX[combo.name]
-                    node.RxnFX[combo.name] += k*DX
-                elif node.spring_DY[0] != None and node.spring_DY[2] == True:
-                    sign = node.spring_DY[1]
-                    k = node.spring_DY[0]
-                    if sign != None: k = float(sign + str(k))
-                    DY = node.DY[combo.name]
-                    node.RxnFY[combo.name] += k*DY
-                elif node.spring_DZ[0] != None and node.spring_DZ[2] == True:
-                    sign = node.spring_DZ[1]
-                    k = node.spring_DZ[0]
-                    if sign != None: k = float(sign + str(k))
-                    DZ = node.DZ[combo.name]
-                    node.RxnFZ[combo.name] += k*DZ
-                elif node.spring_RX[0] != None and node.spring_RX[2] == True:
-                    sign = node.spring_RX[1]
-                    k = node.spring_RX[0]
-                    if sign != None: k = float(sign + str(k))
-                    RX = node.RX[combo.name]
-                    node.RxnMX[combo.name] += k*RX
-                elif node.spring_RY[0] != None and node.spring_RY[2] == True:
-                    sign = node.spring_RY[1]
-                    k = node.spring_RY[0]
-                    if sign != None: k = float(sign + str(k))
-                    RY = node.RY[combo.name]
-                    node.RxnMY[combo.name] += k*RY
-                elif node.spring_RZ[0] != None and node.spring_RZ[2] == True:
-                    sign = node.spring_RZ[1]
-                    k = node.spring_RZ[0]
-                    if sign != None: k = float(sign + str(k))
-                    RZ = node.RZ[combo.name]
-                    node.RxnMZ[combo.name] += k*RZ
-
-    def _check_stability(self, K):
-        """
-        Identifies nodal instabilities in the model's stiffness matrix.
-        """
-
-        # Initialize the `unstable` flag to `False`
-        unstable = False
-
-        # Step through each diagonal term in the stiffness matrix
-        for i in range(K.shape[0]):
-            
-            # Determine which node this term belongs to
-            node = [node for node in self.Nodes.values() if node.ID == int(i/6)][0]
-
-            # Determine which degree of freedom this term belongs to
-            dof = i%6
-
-            # Check to see if this degree of freedom is supported
-            if dof == 0:
-                supported = node.support_DX
-            elif dof == 1:
-                supported = node.support_DY
-            elif dof == 2:
-                supported = node.support_DZ
-            elif dof == 3:
-                supported = node.support_RX
-            elif dof == 4:
-                supported = node.support_RY
-            elif dof == 5:
-                supported = node.support_RZ
-
-            # Check if the degree of freedom on this diagonal is unstable
-            if isclose(K[i, i], 0) and not supported:
-
-                # Flag the model as unstable
-                unstable = True
-
-                # Identify which direction this instability effects
-                if i%6 == 0: direction = 'for translation in the global X direction.'
-                if i%6 == 1: direction = 'for translation in the global Y direction.'
-                if i%6 == 2: direction = 'for translation in the global Z direction.'
-                if i%6 == 3: direction = 'for rotation about the global X axis.'
-                if i%6 == 4: direction = 'for rotation about the global Y axis.'
-                if i%6 == 5: direction = 'for rotation about the global Z axis.'
-
-                # Print a message to the console
-                print('* Nodal instability detected: node ' + node.name + ' is unstable ' + direction)
-
-        if unstable:
-            raise Exception('Unstable node(s). See console output for details.')
-
-        return
-    
-    def _check_statics(self):
-        '''
-        Checks static equilibrium and prints results to the console.
-
-        Parameters
-        ----------
-        precision : number
-            The number of decimal places to carry the results to.
-        '''
-
-        print('+----------------+')
-        print('| Statics Check: |')
-        print('+----------------+')
-        print('')
-
-        from prettytable import PrettyTable
-
-        # Start a blank table and create a header row
-        statics_table = PrettyTable()
-        statics_table.field_names = ['Load Combination', 'Sum FX', 'Sum RX', 'Sum FY', 'Sum RY', 'Sum FZ', 'Sum RZ', 'Sum MX', 'Sum RMX', 'Sum MY', 'Sum RMY', 'Sum MZ', 'Sum RMZ']
-
-        # Step through each load combination
-        for combo in self.LoadCombos.values():
-
-            # Initialize force and moment summations to zero
-            SumFX, SumFY, SumFZ = 0.0, 0.0, 0.0
-            SumMX, SumMY, SumMZ = 0.0, 0.0, 0.0
-            SumRFX, SumRFY, SumRFZ = 0.0, 0.0, 0.0
-            SumRMX, SumRMY, SumRMZ = 0.0, 0.0, 0.0
-
-            # Get the global force vector and the global fixed end reaction vector
-            P = self.P(combo.name)
-            FER = self.FER(combo.name)
-
-            # Step through each node and sum its forces
-            for node in self.Nodes.values():
-
-                # Get the node's coordinates
-                X = node.X
-                Y = node.Y
-                Z = node.Z
-
-                # Get the nodal forces
-                FX = P[node.ID*6+0][0] - FER[node.ID*6+0][0]
-                FY = P[node.ID*6+1][0] - FER[node.ID*6+1][0]
-                FZ = P[node.ID*6+2][0] - FER[node.ID*6+2][0]
-                MX = P[node.ID*6+3][0] - FER[node.ID*6+3][0]
-                MY = P[node.ID*6+4][0] - FER[node.ID*6+4][0]
-                MZ = P[node.ID*6+5][0] - FER[node.ID*6+5][0]
-
-                # Get the nodal reactions
-                RFX = node.RxnFX[combo.name]
-                RFY = node.RxnFY[combo.name]
-                RFZ = node.RxnFZ[combo.name]
-                RMX = node.RxnMX[combo.name]
-                RMY = node.RxnMY[combo.name]
-                RMZ = node.RxnMZ[combo.name]
-
-                # Sum the global forces
-                SumFX += FX
-                SumFY += FY
-                SumFZ += FZ
-                SumMX += MX - FY*Z + FZ*Y
-                SumMY += MY + FX*Z - FZ*X
-                SumMZ += MZ - FX*Y + FY*X
-
-                # Sum the global reactions
-                SumRFX += RFX
-                SumRFY += RFY
-                SumRFZ += RFZ
-                SumRMX += RMX - RFY*Z + RFZ*Y
-                SumRMY += RMY + RFX*Z - RFZ*X
-                SumRMZ += RMZ - RFX*Y + RFY*X   
-
-            # Add the results to the table
-            statics_table.add_row([combo.name, '{:.3g}'.format(SumFX), '{:.3g}'.format(SumRFX),
-                                               '{:.3g}'.format(SumFY), '{:.3g}'.format(SumRFY),
-                                               '{:.3g}'.format(SumFZ), '{:.3g}'.format(SumRFZ),
-                                               '{:.3g}'.format(SumMX), '{:.3g}'.format(SumRMX),
-                                               '{:.3g}'.format(SumMY), '{:.3g}'.format(SumRMY),
-                                               '{:.3g}'.format(SumMZ), '{:.3g}'.format(SumRMZ)])
-
-        # Print the static check table
-        print(statics_table)
-        print('')
-    
+       
     def _renumber(self):
         """
         Assigns node and element ID numbers to be used internally by the program. Numbers are

From 63afe56f5ea7dbc7a8e6278dafdc1128a14d94d7 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 22 Jul 2023 14:31:00 -0400
Subject: [PATCH 08/20] Added load combination tags

---
 PyNite/Analysis.py  | 20 ++++++++++++----
 PyNite/FEModel3D.py | 58 +++++++++++++++++++++++++++++++++------------
 PyNite/LoadCombo.py |  9 ++++---
 3 files changed, 63 insertions(+), 24 deletions(-)

diff --git a/PyNite/Analysis.py b/PyNite/Analysis.py
index 548ec358..5dddd3d2 100644
--- a/PyNite/Analysis.py
+++ b/PyNite/Analysis.py
@@ -180,7 +180,7 @@ def _check_TC_convergence(model, combo_name='Combo 1', log=True):
     # Return whether the TC analysis has converged
     return convergence
 
-def _calc_reactions(model, log=False):
+def _calc_reactions(model, log=False, combo_tags=None):
     """
     Calculates reactions internally once the model is solved.
 
@@ -193,11 +193,17 @@ def _calc_reactions(model, log=False):
     # Print a status update to the console
     if log: print('- Calculating reactions')
 
+    # Identify which load combinations to evaluate
+    if combo_tags is None:
+        combo_list = model.LoadCombos.values()
+    else:
+        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tag in combo_tags]
+
     # Calculate the reactions node by node
     for node in model.Nodes.values():
         
         # Step through each load combination
-        for combo in model.LoadCombos.values():
+        for combo in combo_list:
             
             # Initialize reactions for this node and load combination
             node.RxnFX[combo.name] = 0.0
@@ -440,7 +446,7 @@ def _calc_reactions(model, log=False):
                 RZ = node.RZ[combo.name]
                 node.RxnMZ[combo.name] += k*RZ
 
-def _check_statics(model):
+def _check_statics(model, combo_tags=None):
     '''
     Checks static equilibrium and prints results to the console.
 
@@ -461,8 +467,14 @@ def _check_statics(model):
     statics_table = PrettyTable()
     statics_table.field_names = ['Load Combination', 'Sum FX', 'Sum RX', 'Sum FY', 'Sum RY', 'Sum FZ', 'Sum RZ', 'Sum MX', 'Sum RMX', 'Sum MY', 'Sum RMY', 'Sum MZ', 'Sum RMZ']
 
+    # Identify which load combinations to evaluate
+    if combo_tags is None:
+        combo_list = model.LoadCombos.values()
+    else:
+        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tag in combo_tags]
+
     # Step through each load combination
-    for combo in model.LoadCombos.values():
+    for combo in combo_list:
 
         # Initialize force and moment summations to zero
         SumFX, SumFY, SumFZ = 0.0, 0.0, 0.0
diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index b9267515..988bdf82 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -1021,19 +1021,19 @@ def def_releases(self, Member, Dxi=False, Dyi=False, Dzi=False, Rxi=False, Ryi=F
         # Flag the model as unsolved
         self.solution = None
 
-    def add_load_combo(self, name, factors, combo_type='strength'):
+    def add_load_combo(self, name, factors, combo_tag='strength'):
         """Adds a load combination to the model.
 
         :param name: A unique name for the load combination (e.g. '1.2D+1.6L+0.5S' or 'Gravity Combo').
         :type name: str
         :param factors: A dictionary containing load cases and their corresponding factors (e.g. {'D':1.2, 'L':1.6, 'S':0.5}).
         :type factors: dict
-        :param combo_type: A description of the type of load combination (e.g. 'strength', 'service'). This has no effect on the analysis. It can be used to mark special combinations for easier filtering through them later on. Defaults to 'service'.
-        :type combo_type: str, optional
+        :param combo_tag: A description of the type of load combination (e.g. 'strength', 'service'). This has no effect on the analysis. It can be used to mark special combinations for easier filtering through them later on. Defaults to 'service'.
+        :type combo_tag: str, optional
         """            
 
         # Create a new load combination object
-        new_combo = LoadCombo(name, combo_type, factors)
+        new_combo = LoadCombo(name, combo_tag, factors)
 
         # Add the load combination to the dictionary of load combinations
         self.LoadCombos[name] = new_combo
@@ -1903,7 +1903,7 @@ def _partition(self, unp_matrix, D1_indices, D2_indices):
             m22 = unp_matrix[D2_indices, :][:, D2_indices]
             return m11, m12, m21, m22
 
-    def analyze(self, log=False, check_stability=True, check_statics=False, max_iter=30, sparse=True):
+    def analyze(self, log=False, check_stability=True, check_statics=False, max_iter=30, sparse=True, combo_tags=None):
         """Performs first-order static analysis. Iterations are performed if tension-only members or compression-only members are present.
 
         :param log: Prints the analysis log to the console if set to True. Default is False.
@@ -1958,8 +1958,17 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
         # Convert D2 from a list to a vector
         D2 = atleast_2d(D2).T
 
+        # Identify which load combinations to evaluate
+        if combo_tags is None:
+            combo_list = self.LoadCombos.values()
+        else:
+            combo_list = []
+            for combo in self.LoadCombos.values():
+                if any([tag in combo.combo_tags for tag in combo_tags]):
+                    combo_list.append(combo)
+
         # Step through each load combination
-        for combo in self.LoadCombos.values():
+        for combo in combo_list:
 
             if log:
                 print('')
@@ -2070,7 +2079,7 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
                 iter_count += 1
 
         # Calculate reactions
-        Analysis._calc_reactions(self)
+        Analysis._calc_reactions(self, log, combo_tags)
 
         if log:
             print('')     
@@ -2079,12 +2088,12 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
 
         # Check statics if requested
         if check_statics == True:
-            Analysis._check_statics(self)
+            Analysis._check_statics(self, combo_tags)
         
         # Flag the model as solved
         self.solution = 'Linear TC'
 
-    def analyze_linear(self, log=False, check_stability=True, check_statics=False, sparse=True):
+    def analyze_linear(self, log=False, check_stability=True, check_statics=False, sparse=True, combo_tags=None):
         """Performs first-order static analysis. This analysis procedure is much faster since it only assembles the global stiffness matrix once, rather than once for each load combination. It is not appropriate when non-linear behavior such as tension/compression only analysis or P-Delta analysis are required.
 
         :param log: Prints the analysis log to the console if set to True. Default is False.
@@ -2137,14 +2146,24 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
         D2 = atleast_2d(D2).T
 
         # Get the partitioned global stiffness matrix K11, K12, K21, K22
+        # Note that for linear analysis the stiffness matrix can be obtained for any load combination, as it's the same for all of them
         combo_name = list(self.LoadCombos.keys())[0]
         if sparse == True:
             K11, K12, K21, K22 = self._partition(self.K(combo_name, log, check_stability, sparse).tolil(), D1_indices, D2_indices)
         else:
             K11, K12, K21, K22 = self._partition(self.K(combo_name, log, check_stability, sparse), D1_indices, D2_indices)
 
+        # Identify which load combinations to evaluate
+        if combo_tags is None:
+            combo_list = self.LoadCombos.values()
+        else:
+            combo_list = []
+            for combo in self.LoadCombos.values():
+                if any([tag in combo.combo_tags for tag in combo_tags]):
+                    combo_list.append(combo)
+
         # Step through each load combination
-        for combo in self.LoadCombos.values():
+        for combo in combo_list:
 
             if log:
                 print('')
@@ -2225,7 +2244,7 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
                 node.RZ[combo.name] = D[node.ID*6 + 5, 0]
 
         # Calculate reactions
-        Analysis._calc_reactions(self)
+        Analysis._calc_reactions(self, log, combo_tags)
 
         if log:
             print('')     
@@ -2234,12 +2253,12 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
 
         # Check statics if requested
         if check_statics == True:
-            Analysis._check_statics(self)
+            Analysis._check_statics(self, combo_tags)
         
         # Flag the model as solved
         self.solution = 'Linear'
 
-    def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01, sparse=True):
+    def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01, sparse=True, combo_tags=None):
         """Performs second order (P-Delta) analysis. This type of analysis is appropriate for most models using beams, columns and braces. Second order analysis is usually required by material specific codes. The analysis is iterative and takes longer to solve. Models with slender members and/or members with combined bending and axial loads will generally have more significant P-Delta effects. P-Delta effects in plates/quads are not considered.
 
         :param log: Prints updates to the console if set to True. Default is False.
@@ -2295,8 +2314,17 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         # Convert D2 from a list to a matrix
         D2 = array(D2, ndmin=2).T
 
+        # Identify which load combinations to evaluate
+        if combo_tags is None:
+            combo_list = self.LoadCombos.values()
+        else:
+            combo_list = []
+            for combo in self.LoadCombos.values():
+                if any([tag in combo.combo_tags for tag in combo_tags]):
+                    combo_list.append(combo)
+
         # Step through each load combination
-        for combo in self.LoadCombos.values():
+        for combo in combo_list:
             
             if log:
                 print('')
@@ -2493,7 +2521,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
                 iter_count_PD += 1
         
         # Calculate reactions
-        Analysis._calc_reactions(self)
+        Analysis._calc_reactions(self, log, combo_tags)
 
         if log:
             print('')
diff --git a/PyNite/LoadCombo.py b/PyNite/LoadCombo.py
index 56b82f4f..07896290 100644
--- a/PyNite/LoadCombo.py
+++ b/PyNite/LoadCombo.py
@@ -2,20 +2,19 @@ class LoadCombo():
     """A class that stores all the information necessary to define a load combination.
     """
 
-    def __init__(self, name, combo_type=None, factors={}):
+    def __init__(self, name, combo_tag=None, factors={}):
         """Initializes a new load combination.
 
         :param name: A unique name for the load combination.
         :type name: str
-        :param combo_type: The type of load combination. This can be any string you would like to use to categorize your load combinations. It is useful for separating load combinations into strength, service, or overstrength combinations as often required by building codes. This parameter has no effect on the analysis. It's simply a tool to help you filter results later on. Defaults to `None`.
-        :type combo_type: str, optional
+        :param combo_tag: A list of tags for the load combination. This is a list of any strings you would like to use to categorize your load combinations. It is useful for separating load combinations into strength, service, or overstrength combinations as often required by building codes. This parameter has no effect on the analysis, but it can be used to restrict analysis to only the load combinations with the tags you specify.
+        :type combo_tag: list, optional
         :param factors: A dictionary of load case names (`keys`) followed by their load factors (`items`). For example, the load combination 1.2D+1.6L would be represented as follows: `{'D': 1.2, 'L': 1.6}`. Defaults to {}.
         :type factors: dict, optional
         """
         
         self.name = name             # A unique user-defined name for the load combination
-        self.combo_type = combo_type # Used to track what type of load combination (e.g. strength or serviceability)
-                                     # The 'combo_type' is just a placeholder for future features for now
+        self.combo_tag = combo_tag   # Used to categorize the load combination (e.g. strength or serviceability)
         self.factors = factors       # A dictionary containing each load case name and associated load factor
     
     def AddLoadCase(self, case_name, factor):

From f108b21fe8b9f0d86bff566005e5d1381ff89d49 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Tue, 25 Jul 2023 11:14:35 -0400
Subject: [PATCH 09/20] Changed * to @ now that arrays are used

---
 PyNite/Plate3D.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/PyNite/Plate3D.py b/PyNite/Plate3D.py
index 0992bf76..26bce7fb 100644
--- a/PyNite/Plate3D.py
+++ b/PyNite/Plate3D.py
@@ -618,19 +618,19 @@ def shear(self, x, y, combo_name='Combo 1'):
         # Calculate the derivatives of the plate moments needed to compute shears
         dMx_dx = (Db*array([[0, 0, 0, 0, 0, 0, -6, 0, 0, 0, -6*y, 0],
                              [0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*y],
-                             [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]])*a)[0]
+                             [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]])@a)[0]
 
         dMxy_dy = (Db*array([[0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*x, 0],
                               [0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, -6*x],
-                              [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]])*a)[2]
+                              [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]])@a)[2]
         
         dMy_dy = (Db*array([[0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*x, 0],
                              [0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, -6*x],
-                             [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]])*a)[1]
+                             [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]])@a)[1]
 
         dMxy_dx = (Db*array([[0, 0, 0, 0, 0, 0, -6, 0, 0, 0, -6*y, 0],
                               [0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*y],
-                              [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]])*a)[2]
+                              [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]])@a)[2]
         
         # Calculate internal shears
         Qx = (dMx_dx + dMxy_dy)[0, 0]

From c0671958cfea6a02a60a10a724e43aa9a6ac913d Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Tue, 25 Jul 2023 12:02:42 -0400
Subject: [PATCH 10/20] Changed * to @ for numpy array multiplication

---
 PyNite/Plate3D.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/PyNite/Plate3D.py b/PyNite/Plate3D.py
index 26bce7fb..dbb3443b 100644
--- a/PyNite/Plate3D.py
+++ b/PyNite/Plate3D.py
@@ -616,25 +616,25 @@ def shear(self, x, y, combo_name='Combo 1'):
         a = self._a(combo_name)
 
         # Calculate the derivatives of the plate moments needed to compute shears
-        dMx_dx = (Db*array([[0, 0, 0, 0, 0, 0, -6, 0, 0, 0, -6*y, 0],
+        dMx_dx = (Db @ array([[0, 0, 0, 0, 0, 0, -6, 0, 0, 0, -6*y, 0],
                              [0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*y],
-                             [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]])@a)[0]
+                             [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]]) @ a)[0]
 
-        dMxy_dy = (Db*array([[0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*x, 0],
+        dMxy_dy = (Db @ array([[0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*x, 0],
                               [0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, -6*x],
-                              [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]])@a)[2]
+                              [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]]) @ a)[2]
         
-        dMy_dy = (Db*array([[0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*x, 0],
+        dMy_dy = (Db @ array([[0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*x, 0],
                              [0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, -6*x],
-                             [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]])@a)[1]
+                             [0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*y]]) @ a)[1]
 
-        dMxy_dx = (Db*array([[0, 0, 0, 0, 0, 0, -6, 0, 0, 0, -6*y, 0],
+        dMxy_dx = (Db @ array([[0, 0, 0, 0, 0, 0, -6, 0, 0, 0, -6*y, 0],
                               [0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -6*y],
-                              [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]])@a)[2]
+                              [0, 0, 0, 0, 0, 0, 0, -4, 0, 0, -12*x, 0]]) @ a)[2]
         
         # Calculate internal shears
-        Qx = (dMx_dx + dMxy_dy)[0, 0]
-        Qy = (dMy_dy + dMxy_dx)[0, 0]
+        Qx = (dMx_dx + dMxy_dy)[0]
+        Qy = (dMy_dy + dMxy_dx)[0]
 
         # Return internal shears
         return array([[Qx], 

From 03147465a009a8e70427e9336caaaf07fa7d7968 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Tue, 25 Jul 2023 21:15:48 -0400
Subject: [PATCH 11/20] Improvements for 'print' theme

---
 PyNite/Visualization.py | 78 +++++++++++++++++++++++++++++------------
 1 file changed, 55 insertions(+), 23 deletions(-)

diff --git a/PyNite/Visualization.py b/PyNite/Visualization.py
index 6f042c31..40f52fef 100644
--- a/PyNite/Visualization.py
+++ b/PyNite/Visualization.py
@@ -267,6 +267,7 @@ def update(self, reset_camera=True):
             renderer.AddActor(vis_member.actor)
 
             if self.labels == True:
+                
                 # Add the actor for the member label
                 renderer.AddActor(vis_member.lblActor)
             
@@ -304,6 +305,10 @@ def update(self, reset_camera=True):
             node_mapper.SetInputConnection(node_polydata.GetOutputPort())
             node_actor = vtk.vtkActor()
             node_actor.SetMapper(node_mapper)
+
+            # While the `VisNode` object we created earlier does store the color information for the node, it will be lost when we combine all the `VisNode` objects into one object. We need to re-adjust the color of the nodes here.
+            if self.theme == 'print':
+                node_actor.GetProperty().SetColor(0, 0, 0) # black
             
             # Add the node actor to the renderer
             renderer.AddActor(node_actor)
@@ -325,7 +330,7 @@ def update(self, reset_camera=True):
 
         # Render the deformed shape if requested
         if self.deformed_shape == True:
-            _DeformedShape(self.model, renderer, self.deformed_scale, self.annotation_size, self.combo_name)
+            _DeformedShape(self.model, renderer, self.deformed_scale, self.annotation_size, self.combo_name, self.render_nodes, self.theme)
 
         # Render the loads if requested
         if (self.combo_name != None or self.case != None) and self.render_loads != False:
@@ -335,15 +340,20 @@ def update(self, reset_camera=True):
         if self.model.Quads or self.model.Plates:
             _RenderContours(self.model, renderer, self.deformed_shape, self.deformed_scale,
                             self.color_map, self.scalar_bar, self.scalar_bar_text_size,
-                            self.combo_name)
+                            self.combo_name, self.theme)
 
-        # Set the window's background to gray
-        renderer.SetBackground(0/255, 0/255, 128/255)
+        # Set the window's background color
+        if self.theme == 'default':
+            renderer.SetBackground(0/255, 0/255, 128/255)  # Blue
+        elif self.theme == 'print':
+            renderer.SetBackground(255, 255, 255)  # White
         
         # Reset the camera
         if reset_camera: renderer.ResetCamera()
 
-# %%
+#%%
+# The code in this section will be deprecated at some point
+
 def RenderModel(model, annotation_size=5, deformed_shape=False, deformed_scale=30,
                 render_loads=True, color_map=None, combo_name='Combo 1', case=None, labels=True,
                 screenshot=None):
@@ -524,7 +534,7 @@ def render_model(model, annotation_size=5, deformed_shape=False, deformed_scale=
 
     # Render the deformed shape if requested
     if deformed_shape == True:
-        _DeformedShape(model, renderer, deformed_scale, annotation_size, combo_name)
+        _DeformedShape(model, renderer, deformed_scale, annotation_size, combo_name, render_nodes=True, theme=theme)
 
     # Render the loads if requested
     if (combo_name != None or case != None) and render_loads != False:
@@ -579,6 +589,7 @@ def render_model(model, annotation_size=5, deformed_shape=False, deformed_scale=
             writer.SetFileName(screenshot)
             writer.Write()
 
+#%%
 # Converts a node object into a node for the viewer
 class VisNode():
 
@@ -847,6 +858,9 @@ def __init__(self, node, annotation_size=5, color=None):
         mapper = vtk.vtkPolyDataMapper()
         mapper.SetInputConnection(self.polydata.GetOutputPort())
         self.actor = vtk.vtkActor()
+        
+        # Set the mapper for the node's actor
+        self.actor.SetMapper(mapper)
       
         # Add color to the actors
         if color == 'red':
@@ -858,9 +872,6 @@ def __init__(self, node, annotation_size=5, color=None):
         elif color == 'black':
             self.actor.GetProperty().SetColor(0, 0, 0)  # Black
             self.lblActor.GetProperty().SetColor(0, 0, 0)  # Black
-        
-        # Set the mapper for the node's actor
-        self.actor.SetMapper(mapper)
 
 class VisSpring():
     
@@ -1327,7 +1338,7 @@ class VisAreaLoad():
     Creates an area load for the viewer
     '''
     
-    def __init__(self, position0, position1, position2, position3, direction, length, label_text, annotation_size=5):
+    def __init__(self, position0, position1, position2, position3, direction, length, label_text, annotation_size=5, theme='default'):
         '''
         Constructor
         '''
@@ -1358,6 +1369,12 @@ def __init__(self, position0, position1, position2, position3, direction, length
         # Add a label
         self.label_actor = ptLoads[0].lblActor
 
+        # Add color to the area load label
+        if theme == 'print':
+            self.label_actor.GetProperty().SetColor(255, 0, 0)  # red
+        elif theme == 'default':
+            self.label_actor.GetProperty().SetColor(0, 255, 0)  # green
+
 def _PerpVector(v):
     '''
     Returns a unit vector perpendicular to v=[i, j, k]
@@ -1473,7 +1490,7 @@ def _PrepContour(model, stress_type='Mx', combo_name='Combo 1'):
             if node.contour != []:
                 node.contour = sum(node.contour)/len(node.contour)
 
-def _DeformedShape(model, renderer, scale_factor, annotation_size, combo_name, color=None):
+def _DeformedShape(model, vtk_renderer, scale_factor, annotation_size, combo_name, render_nodes=True, theme='default'):
     '''
     Renders the deformed shape of a model.
     
@@ -1501,7 +1518,7 @@ def _DeformedShape(model, renderer, scale_factor, annotation_size, combo_name, c
     append_filter = vtk.vtkAppendPolyData()
     
     # Check if nodes are to be rendered
-    if renderer.render_nodes == True:
+    if render_nodes == True:
         
         # Add the deformed nodes to the append filter
         for node in model.Nodes.values():
@@ -1531,10 +1548,10 @@ def _DeformedShape(model, renderer, scale_factor, annotation_size, combo_name, c
     mapper = vtk.vtkPolyDataMapper()
     mapper.SetInputConnection(append_filter.GetOutputPort())
     actor = vtk.vtkActor()
-    if color is None: actor.GetProperty().SetColor(255, 255, 0)  # Yellow
-    elif color == 'black': actor.GetProperty().SetColor(0, 0, 0)  # Black
+    if theme == 'default': actor.GetProperty().SetColor(255, 255, 0)  # Yellow
+    elif theme == 'print': actor.GetProperty().SetColor(0, 0, 0)  # Black
     actor.SetMapper(mapper)
-    renderer.AddActor(actor)
+    vtk_renderer.AddActor(actor)
 
 def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='default'):
 
@@ -1715,7 +1732,7 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
                 position3 = [plate.n_node.X, plate.n_node.Y, plate.n_node.Z]
           
                 # Create an area load and get its data
-                area_load = VisAreaLoad(position0, position1, position2, position3, dir_cos*sign, abs(load_value)/max_area_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size)
+                area_load = VisAreaLoad(position0, position1, position2, position3, dir_cos*sign, abs(load_value)/max_area_load*5*annotation_size, '{:.3g}'.format(load_value), annotation_size, theme)
           
                 # Add the area load's arrows to the overall load polydata
                 polydata.AddInputData(area_load.polydata.GetOutput())
@@ -1757,23 +1774,33 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
     load_mapper = vtk.vtkPolyDataMapper()
     load_mapper.SetInputConnection(polydata.GetOutputPort())
     load_actor = vtk.vtkActor()
-    if theme != 'print': load_actor.GetProperty().SetColor(0, 255, 0)  # Green
-    else: load_actor.GetProperty().SetColor(0, 0, 0)  # Black
     load_actor.SetMapper(load_mapper)
+
+    # Colorize the loads
+    if theme == 'default':
+        load_actor.GetProperty().SetColor(0, 255, 0)  # Green
+    elif theme == 'print':
+        load_actor.GetProperty().SetColor(255, 0, 0)  # Red
+
+    # Add the load actor to the renderer
     renderer.AddActor(load_actor)
     
     # Set up an actor and a mapper for the area load polygons
     polygon_mapper = vtk.vtkPolyDataMapper()
     polygon_mapper.SetInputData(polygon_polydata)
     polygon_actor = vtk.vtkActor()
-    if theme != 'print': polygon_actor.GetProperty().SetColor(0, 255, 0)  # Green
-    else: polygon_actor.GetProperty().SetColor(128, 128, 128)  # Grey
+
     # polygon_actor.GetProperty().SetOpacity(0.5)      # 50% opacity
     polygon_actor.SetMapper(polygon_mapper)
     renderer.AddActor(polygon_actor)
+    
+    # Set the color of the area load polygons
+    if theme == 'default':
+        polygon_actor.GetProperty().SetColor(0, 255, 0)  # Green
+    elif theme == 'print':
+        polygon_actor.GetProperty().SetColor(255, 0, 0)  # Red
 
-def _RenderContours(model, renderer, deformed_shape, deformed_scale, color_map, scalar_bar,
-                    scalar_bar_text_size, combo_name):
+def _RenderContours(model, renderer, deformed_shape, deformed_scale, color_map, scalar_bar, scalar_bar_text_size, combo_name, theme='default'):
   
     # Create a new `vtkCellArray` object to store the elements
     plates = vtk.vtkCellArray()
@@ -1781,7 +1808,7 @@ def _RenderContours(model, renderer, deformed_shape, deformed_scale, color_map,
     # Create a `vtkPoints` object to store the coordinates of the corners of the elements
     plate_points = vtk.vtkPoints()
     
-    # Create 2 lists to store plate results
+    # Create 2 lists to store plate result
     # `results` will store the results in a Python iterable list
     # `plate_results` will store the results in a `vtkDoubleArray` for VTK
     results = []
@@ -1907,6 +1934,11 @@ def _RenderContours(model, renderer, deformed_shape, deformed_scale, color_map,
             scalar_text = vtk.vtkTextProperty()
             scalar_text.SetFontSize(max(int(scalar_bar_text_size), 1))
             scalar_text.SetBold(1)
+
+            # The `vtkTextProperty` object is white by default
+            if theme == 'print':
+                scalar_text.SetColor(255, 255, 255)  # Black
+            
             scalar.SetLabelTextProperty(scalar_text)
 
             scalar.SetMaximumWidthInPixels(100)

From b26bdba229534ea1273801645581769872922900 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Tue, 25 Jul 2023 21:18:29 -0400
Subject: [PATCH 12/20] Updates to load combo tags

---
 Examples/Beam on Elastic Foundation.py | 2 +-
 PyNite/Analysis.py                     | 4 ++--
 PyNite/FEModel3D.py                    | 8 ++++----
 PyNite/LoadCombo.py                    | 8 ++++----
 docs/source/load_combo.rst             | 2 +-
 5 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/Examples/Beam on Elastic Foundation.py b/Examples/Beam on Elastic Foundation.py
index c8bb3dbe..ef02bf72 100644
--- a/Examples/Beam on Elastic Foundation.py	
+++ b/Examples/Beam on Elastic Foundation.py	
@@ -62,7 +62,7 @@
 # Analyze the model. PyNite's standard solver is most appropriate or this model since there are
 # non-linear features (compression-only springs) but no large axial forces that would cause P-Delta
 # effects.
-boef.analyze(check_statics=True)
+boef.analyze(log=True, check_statics=True)
 
 # Render the mdoel with the deformed shape using PyNite's buit-in renderer
 from PyNite.Visualization import Renderer
diff --git a/PyNite/Analysis.py b/PyNite/Analysis.py
index 5dddd3d2..a92f445f 100644
--- a/PyNite/Analysis.py
+++ b/PyNite/Analysis.py
@@ -197,7 +197,7 @@ def _calc_reactions(model, log=False, combo_tags=None):
     if combo_tags is None:
         combo_list = model.LoadCombos.values()
     else:
-        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tag in combo_tags]
+        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tags in combo_tags]
 
     # Calculate the reactions node by node
     for node in model.Nodes.values():
@@ -471,7 +471,7 @@ def _check_statics(model, combo_tags=None):
     if combo_tags is None:
         combo_list = model.LoadCombos.values()
     else:
-        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tag in combo_tags]
+        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tags in combo_tags]
 
     # Step through each load combination
     for combo in combo_list:
diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index 988bdf82..d99c3d62 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -1021,19 +1021,19 @@ def def_releases(self, Member, Dxi=False, Dyi=False, Dzi=False, Rxi=False, Ryi=F
         # Flag the model as unsolved
         self.solution = None
 
-    def add_load_combo(self, name, factors, combo_tag='strength'):
+    def add_load_combo(self, name, factors, combo_tags='strength'):
         """Adds a load combination to the model.
 
         :param name: A unique name for the load combination (e.g. '1.2D+1.6L+0.5S' or 'Gravity Combo').
         :type name: str
         :param factors: A dictionary containing load cases and their corresponding factors (e.g. {'D':1.2, 'L':1.6, 'S':0.5}).
         :type factors: dict
-        :param combo_tag: A description of the type of load combination (e.g. 'strength', 'service'). This has no effect on the analysis. It can be used to mark special combinations for easier filtering through them later on. Defaults to 'service'.
-        :type combo_tag: str, optional
+        :param combo_tags: A description of the type of load combination (e.g. 'strength', 'service'). This has no effect on the analysis. It can be used to mark special combinations for easier filtering through them later on. Defaults to 'service'.
+        :type combo_tags: str, optional
         """            
 
         # Create a new load combination object
-        new_combo = LoadCombo(name, combo_tag, factors)
+        new_combo = LoadCombo(name, combo_tags, factors)
 
         # Add the load combination to the dictionary of load combinations
         self.LoadCombos[name] = new_combo
diff --git a/PyNite/LoadCombo.py b/PyNite/LoadCombo.py
index 07896290..2e2d0ae3 100644
--- a/PyNite/LoadCombo.py
+++ b/PyNite/LoadCombo.py
@@ -2,19 +2,19 @@ class LoadCombo():
     """A class that stores all the information necessary to define a load combination.
     """
 
-    def __init__(self, name, combo_tag=None, factors={}):
+    def __init__(self, name, combo_tags=None, factors={}):
         """Initializes a new load combination.
 
         :param name: A unique name for the load combination.
         :type name: str
-        :param combo_tag: A list of tags for the load combination. This is a list of any strings you would like to use to categorize your load combinations. It is useful for separating load combinations into strength, service, or overstrength combinations as often required by building codes. This parameter has no effect on the analysis, but it can be used to restrict analysis to only the load combinations with the tags you specify.
-        :type combo_tag: list, optional
+        :param combo_tags: A list of tags for the load combination. This is a list of any strings you would like to use to categorize your load combinations. It is useful for separating load combinations into strength, service, or overstrength combinations as often required by building codes. This parameter has no effect on the analysis, but it can be used to restrict analysis to only the load combinations with the tags you specify.
+        :type combo_tags: list, optional
         :param factors: A dictionary of load case names (`keys`) followed by their load factors (`items`). For example, the load combination 1.2D+1.6L would be represented as follows: `{'D': 1.2, 'L': 1.6}`. Defaults to {}.
         :type factors: dict, optional
         """
         
         self.name = name             # A unique user-defined name for the load combination
-        self.combo_tag = combo_tag   # Used to categorize the load combination (e.g. strength or serviceability)
+        self.combo_tags = combo_tags   # Used to categorize the load combination (e.g. strength or serviceability)
         self.factors = factors       # A dictionary containing each load case name and associated load factor
     
     def AddLoadCase(self, case_name, factor):
diff --git a/docs/source/load_combo.rst b/docs/source/load_combo.rst
index bd59e6f7..71e26d6d 100644
--- a/docs/source/load_combo.rst
+++ b/docs/source/load_combo.rst
@@ -22,6 +22,6 @@ must be checked is dead loads factored by 1.2 acting simultaneously with live lo
     # Assume 'D' and 'L' are names previously specified for load cases
     my_model.add_load_combo('Vertical Loads', {'D':1.2, 'L':1.6})
 
-Load combinations can also be passed a 3rd optional argument named ``combo_type``. This has no effect on
+Load combinations can also be passed a 3rd optional argument named ``combo_tags``. This has no effect on
 the analysis, but it can be used to categorize load combinations (e.g. 'strength' or 'service') for
 easier filtering of results later on.
\ No newline at end of file

From cd7b5ba91a2e5ce63aa6a9a8e051cd98720a0145 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 29 Jul 2023 19:22:23 -0400
Subject: [PATCH 13/20] Worked on "themes" for visualization

---
 PyNite/Visualization.py | 99 +++++++++++++++++++++++++----------------
 1 file changed, 61 insertions(+), 38 deletions(-)

diff --git a/PyNite/Visualization.py b/PyNite/Visualization.py
index 40f52fef..a8244e9a 100644
--- a/PyNite/Visualization.py
+++ b/PyNite/Visualization.py
@@ -216,8 +216,10 @@ def update(self, reset_camera=True):
         # Check if nodes are to be rendered
         if self.render_nodes == True:
 
-            if self.theme == 'print': color = 'black'
-            else: color = None
+            if self.theme == 'print':
+                color = 'black'
+            else:
+                color = None
 
             # Create a visual node for each node in the model
             vis_nodes = []
@@ -237,7 +239,7 @@ def update(self, reset_camera=True):
         # Create a visual member for each member in the model
         vis_members = []
         for member in self.model.Members.values():
-            vis_members.append(VisMember(member, self.model.Nodes, self.annotation_size))
+            vis_members.append(VisMember(member, self.model.Nodes, self.annotation_size, self.theme))
         
         # Get the renderer
         renderer = self.renderer
@@ -289,6 +291,9 @@ def update(self, reset_camera=True):
                 node_polydata.AddInputData(vis_node.polydata.GetOutput())
 
                 if self.labels == True:
+
+                    if self.theme == 'print':
+                        vis_node.lblActor.GetProperty().SetColor(0, 0, 0)  # black
                     
                     # Add the actor for the node label
                     renderer.AddActor(vis_node.lblActor)
@@ -306,16 +311,16 @@ def update(self, reset_camera=True):
             node_actor = vtk.vtkActor()
             node_actor.SetMapper(node_mapper)
 
-            # While the `VisNode` object we created earlier does store the color information for the node, it will be lost when we combine all the `VisNode` objects into one object. We need to re-adjust the color of the nodes here.
+            # Adjust the color of the nodes here.
             if self.theme == 'print':
-                node_actor.GetProperty().SetColor(0, 0, 0) # black
+                node_actor.GetProperty().SetColor(0, 0, 0)  # Black
             
             # Add the node actor to the renderer
             renderer.AddActor(node_actor)
 
-            # Add actors for each auxiliary node
+            # Add actors for each auxiliary node. There are usually only a few of these, so they shouldn't bog down interaction if we do them individually rather than with a `vtkAppendFilter`
             for vis_aux_node in vis_aux_nodes:
-            
+
                 # Add the actor for the auxiliary node
                 renderer.AddActor(vis_aux_node.actor)
 
@@ -344,7 +349,7 @@ def update(self, reset_camera=True):
 
         # Set the window's background color
         if self.theme == 'default':
-            renderer.SetBackground(0/255, 0/255, 128/255)  # Blue
+            renderer.SetBackground(0, 0, 128)  # Blue
         elif self.theme == 'print':
             renderer.SetBackground(255, 255, 255)  # White
         
@@ -421,9 +426,11 @@ def render_model(model, annotation_size=5, deformed_shape=False, deformed_scale=
         raise Exception('Unable to render load combination. No load combinations defined.')
     
     # Create a visual node for each node in the model
+    if theme == 'print': color = 'black'
+    else: color = None
     vis_nodes = []
     for node in model.Nodes.values():
-        vis_nodes.append(VisNode(node, annotation_size))
+        vis_nodes.append(VisNode(node, annotation_size, color))
     
     # Create a visual auxiliary node for each auxiliary node in the model
     vis_aux_nodes = []
@@ -438,7 +445,7 @@ def render_model(model, annotation_size=5, deformed_shape=False, deformed_scale=
     # Create a visual member for each member in the model
     vis_members = []
     for member in model.Members.values():
-        vis_members.append(VisMember(member, model.Nodes, annotation_size))
+        vis_members.append(VisMember(member, model.Nodes, annotation_size, theme))
     
     # Create a window
     window = vtk.vtkRenderWindow()
@@ -861,17 +868,21 @@ def __init__(self, node, annotation_size=5, color=None):
         
         # Set the mapper for the node's actor
         self.actor.SetMapper(mapper)
-      
-        # Add color to the actors
-        if color == 'red':
-            self.actor.GetProperty().SetColor(255, 0, 0)  # Red
-            self.lblActor.GetProperty().SetColor(255, 0, 0)  # Red
-        elif color == 'yellow':
-            self.actor.GetProperty().SetColor(255, 255, 0)  # Yellow
-            self.lblActor.GetProperty().SetColor(255, 255, 0)  # Yellow
-        elif color == 'black':
-            self.actor.GetProperty().SetColor(0, 0, 0)  # Black
-            self.lblActor.GetProperty().SetColor(0, 0, 0)  # Black
+
+        # Color will be added to the node actors outside of this class once they are all combined into one `vtkAppendFilter` object.
+
+        # TODO: Delete legacy code below once it's been proven to be unnecessary over time.
+
+        # Add color to the node and label actors if specified
+        # if color == 'red':
+        #     self.actor.GetProperty().SetColor(255, 0, 0)  # Red
+        #     self.lblActor.GetProperty().SetColor(255, 0, 0)  # Red
+        # elif color == 'yellow':
+        #     self.actor.GetProperty().SetColor(255, 255, 0)  # Yellow
+        #     self.lblActor.GetProperty().SetColor(255, 255, 0)  # Yellow
+        # elif color == 'black':
+        #     self.actor.GetProperty().SetColor(0, 0, 0)  # Black
+        #     self.lblActor.GetProperty().SetColor(0, 0, 0)  # Black
 
 class VisSpring():
     
@@ -902,10 +913,8 @@ def __init__(self, spring, nodes, annotation_size=5, color=None):
         mapper = vtk.vtkPolyDataMapper()
         mapper.SetInputConnection(line.GetOutputPort())
 
-        # Set up an actor for the spring
+        # Set up an actor and a mapper for the spring
         self.actor = vtk.vtkActor()
-        if color is None: self.actor.GetProperty().SetColor(255, 0, 255)  # Magenta
-        elif color == 'black': self.actor.GetProperty().SetColor(0, 0, 0)  # Black
         self.actor.SetMapper(mapper)
 
         # Create the text for the spring label
@@ -921,12 +930,20 @@ def __init__(self, spring, nodes, annotation_size=5, color=None):
         self.lblActor.SetMapper(lblMapper)
         self.lblActor.SetScale(annotation_size, annotation_size, annotation_size)
         self.lblActor.SetPosition((Xi+Xj)/2, (Yi+Yj)/2, (Zi+Zj)/2)
+
+        # Add some color
+        if color is None:
+            self.actor.GetProperty().SetColor(255, 0, 255)     # Magenta
+            self.lblActor.GetProperty().SetColor(255, 0, 255)
+        elif color == 'black':
+            self.actor.GetProperty().SetColor(0, 0, 0)         # Black
+            self.lblActor.GetProperty().SetColor(0, 0, 0)
       
 # Converts a member object into a member for the viewer
 class VisMember():
 
     # Constructor
-    def __init__(self, member, nodes, annotation_size=5, color=None):
+    def __init__(self, member, nodes, annotation_size=5, theme='default'):
     
         # Generate a line for the member
         line = vtk.vtkLineSource()
@@ -970,10 +987,10 @@ def __init__(self, member, nodes, annotation_size=5, color=None):
         self.lblActor.SetScale(annotation_size, annotation_size, annotation_size)
         self.lblActor.SetPosition((Xi+Xj)/2, (Yi+Yj)/2, (Zi+Zj)/2)
 
-        # Adjust the color of the member
-        if color == 'black':
-            self.actor.GetProperty().SetColor(0, 0, 0)  # Black
-            self.lblActor.GetProperty().SetColor(0, 0, 0) # Black
+        # Adjust the color of the member if the theme is 'print'
+        if theme == 'print':
+            self.actor.GetProperty().SetColor(0/255, 0/255, 0/255)  # Black
+            self.lblActor.GetProperty().SetColor(0/255, 0/255, 0/255) # Black
 
 # Converts a node object into a node in its deformed position for the viewer
 class VisDeformedNode():
@@ -1371,9 +1388,9 @@ def __init__(self, position0, position1, position2, position3, direction, length
 
         # Add color to the area load label
         if theme == 'print':
-            self.label_actor.GetProperty().SetColor(255, 0, 0)  # red
+            self.label_actor.GetProperty().SetColor(255/255, 0/255, 0/255)  # red
         elif theme == 'default':
-            self.label_actor.GetProperty().SetColor(0, 255, 0)  # green
+            self.label_actor.GetProperty().SetColor(0/255, 255/255, 0/255)  # green
 
 def _PerpVector(v):
     '''
@@ -1548,9 +1565,15 @@ def _DeformedShape(model, vtk_renderer, scale_factor, annotation_size, combo_nam
     mapper = vtk.vtkPolyDataMapper()
     mapper.SetInputConnection(append_filter.GetOutputPort())
     actor = vtk.vtkActor()
-    if theme == 'default': actor.GetProperty().SetColor(255, 255, 0)  # Yellow
-    elif theme == 'print': actor.GetProperty().SetColor(0, 0, 0)  # Black
     actor.SetMapper(mapper)
+
+    # Adjust the color
+    if theme == 'default':
+        actor.GetProperty().SetColor(255/255, 255/255, 0/255)  # Yellow
+    elif theme == 'print':
+        actor.GetProperty().SetColor(26/255, 26/255, 26/255)  # Dark Grey
+    
+    # Add the actor to the renderer
     vtk_renderer.AddActor(actor)
 
 def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='default'):
@@ -1778,9 +1801,9 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
 
     # Colorize the loads
     if theme == 'default':
-        load_actor.GetProperty().SetColor(0, 255, 0)  # Green
+        load_actor.GetProperty().SetColor(0/255, 255/255, 0/255)  # Green
     elif theme == 'print':
-        load_actor.GetProperty().SetColor(255, 0, 0)  # Red
+        load_actor.GetProperty().SetColor(255/255, 0/255, 0/255)  # Red
 
     # Add the load actor to the renderer
     renderer.AddActor(load_actor)
@@ -1796,9 +1819,9 @@ def _RenderLoads(model, renderer, annotation_size, combo_name, case, theme='defa
     
     # Set the color of the area load polygons
     if theme == 'default':
-        polygon_actor.GetProperty().SetColor(0, 255, 0)  # Green
+        polygon_actor.GetProperty().SetColor(0/255, 255/255, 0/255)  # Green
     elif theme == 'print':
-        polygon_actor.GetProperty().SetColor(255, 0, 0)  # Red
+        polygon_actor.GetProperty().SetColor(255/255, 0/255, 0/255)  # Red
 
 def _RenderContours(model, renderer, deformed_shape, deformed_scale, color_map, scalar_bar, scalar_bar_text_size, combo_name, theme='default'):
   
@@ -1937,7 +1960,7 @@ def _RenderContours(model, renderer, deformed_shape, deformed_scale, color_map,
 
             # The `vtkTextProperty` object is white by default
             if theme == 'print':
-                scalar_text.SetColor(255, 255, 255)  # Black
+                scalar_text.SetColor(255/255, 255/255, 255/255)  # Black
             
             scalar.SetLabelTextProperty(scalar_text)
 

From 442c4070ae834e4daaecfb1c6356eae89329fc63 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 29 Jul 2023 19:38:35 -0400
Subject: [PATCH 14/20] Allowed for multiple tags to a single load combo

---
 PyNite/Analysis.py  | 10 ++++++++--
 PyNite/FEModel3D.py | 12 ++++++------
 2 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/PyNite/Analysis.py b/PyNite/Analysis.py
index a92f445f..c3e3343b 100644
--- a/PyNite/Analysis.py
+++ b/PyNite/Analysis.py
@@ -197,7 +197,10 @@ def _calc_reactions(model, log=False, combo_tags=None):
     if combo_tags is None:
         combo_list = model.LoadCombos.values()
     else:
-        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tags in combo_tags]
+        combo_list = []
+        for combo in model.LoadCombos.values():
+            if any(tag in combo.combo_tags for tag in combo_tags):
+                combo_list.append(combo)
 
     # Calculate the reactions node by node
     for node in model.Nodes.values():
@@ -471,7 +474,10 @@ def _check_statics(model, combo_tags=None):
     if combo_tags is None:
         combo_list = model.LoadCombos.values()
     else:
-        combo_list = [combo for combo in model.LoadCombos.values() if combo.combo_tags in combo_tags]
+        combo_list = []
+        for combo in model.LoadCombos.values():
+            if any(tag in combo.combo_tags for tag in combo_tags):
+                combo_list.append(combo)
 
     # Step through each load combination
     for combo in combo_list:
diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index d99c3d62..35c3d4e8 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -1021,15 +1021,15 @@ def def_releases(self, Member, Dxi=False, Dyi=False, Dzi=False, Rxi=False, Ryi=F
         # Flag the model as unsolved
         self.solution = None
 
-    def add_load_combo(self, name, factors, combo_tags='strength'):
+    def add_load_combo(self, name, factors, combo_tags=None):
         """Adds a load combination to the model.
 
         :param name: A unique name for the load combination (e.g. '1.2D+1.6L+0.5S' or 'Gravity Combo').
         :type name: str
         :param factors: A dictionary containing load cases and their corresponding factors (e.g. {'D':1.2, 'L':1.6, 'S':0.5}).
         :type factors: dict
-        :param combo_tags: A description of the type of load combination (e.g. 'strength', 'service'). This has no effect on the analysis. It can be used to mark special combinations for easier filtering through them later on. Defaults to 'service'.
-        :type combo_tags: str, optional
+        :param combo_tags: A list of tags used to categorize load combinations. Default is `None`. This can be useful for filtering results later on, or for limiting analysis to only those combinations with certain tags. This feature is provided for convenience. It is not necessary to use tags.
+        :type combo_tags: list, optional
         """            
 
         # Create a new load combination object
@@ -1964,7 +1964,7 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
         else:
             combo_list = []
             for combo in self.LoadCombos.values():
-                if any([tag in combo.combo_tags for tag in combo_tags]):
+                if any(tag in combo.combo_tags for tag in combo_tags):
                     combo_list.append(combo)
 
         # Step through each load combination
@@ -2159,7 +2159,7 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
         else:
             combo_list = []
             for combo in self.LoadCombos.values():
-                if any([tag in combo.combo_tags for tag in combo_tags]):
+                if any(tag in combo.combo_tags for tag in combo_tags):
                     combo_list.append(combo)
 
         # Step through each load combination
@@ -2320,7 +2320,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         else:
             combo_list = []
             for combo in self.LoadCombos.values():
-                if any([tag in combo.combo_tags for tag in combo_tags]):
+                if any(tag in combo.combo_tags for tag in combo_tags):
                     combo_list.append(combo)
 
         # Step through each load combination

From a0d7adeca434b1d130810d09eaf465809bee4dab Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 29 Jul 2023 19:44:03 -0400
Subject: [PATCH 15/20] Prep for v0.0.80

---
 README.md | 5 +++++
 setup.py  | 2 +-
 2 files changed, 6 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 15e4522e..0c285cb0 100644
--- a/README.md
+++ b/README.md
@@ -61,6 +61,11 @@ Here's a list of projects that use PyNite:
 * Phaenotyp (https://github.com/bewegende-Architektur/Phaenotyp) (https://youtu.be/shloSw9HjVI)
 
 # What's New?
+v0.0.80
+* Refactored/simplified analysis code. Much of it has been moved to a new `Analysis` file that eliminated redundant code.
+* Load combination tags have replaced `combo_type`. You can now use a list of tags to tag your load combinations for easier categorization.
+* You no longer have to run all load combinations. You can now run select combos based on their tags.
+
 v0.0.79
 * Added the option to turn off nodes during visualization.
 * Bug fix for meshing cylinders about the global X or Z axis.
diff --git a/setup.py b/setup.py
index 3c4d7fd3..6c034703 100644
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="PyNiteFEA",
-    version="0.0.78",
+    version="0.0.80",
     author="D. Craig Brinck, PE, SE",
     author_email="Building.Code@outlook.com",
     description="A simple elastic 3D structural finite element library for Python.",

From b9cfb270ad668109ada99d0d3940783b5c0f34c9 Mon Sep 17 00:00:00 2001
From: Kevin Russell <kevin@ottermatics.com>
Date: Thu, 3 Aug 2023 17:24:08 -0700
Subject: [PATCH 16/20] Add fix for 3D cylinder mesh: >allow origin reading for
 input on coordinate >do not adjust number of node spacing in loop as that can
 result in unpredictable behavior >ensure at least one mesh element will be
 made

---
 PyNite/Mesh.py | 24 ++++++++++++++++++------
 1 file changed, 18 insertions(+), 6 deletions(-)

diff --git a/PyNite/Mesh.py b/PyNite/Mesh.py
index e864fb79..1d737b22 100644
--- a/PyNite/Mesh.py
+++ b/PyNite/Mesh.py
@@ -1208,7 +1208,7 @@ class CylinderMesh(Mesh):
     Parameters
     ----------
     mesh_size : number
-        The desired mesh size. This value will only be used to mesh vertically if `num_elements` is
+        The desired mesh element edge size. This value will only be used to mesh vertically if `num_elements` is
         specified. Otherwise it will be used to mesh the circumference too.
     radius : number
         The radius of the cylinder to the element centers
@@ -1273,7 +1273,14 @@ def generate(self):
 
         radius = self.radius
         h = self.h
-        y = self.origin[1]
+
+        if self.axis == 'Y':
+            y = self.origin[1]
+        elif self.axis == 'X':
+            y = self.origin[0]
+        elif self.axis== 'Z':
+            y = self.origin[2]        
+
         n = int(self.start_node[1:])
         q = int(self.start_element[1:])
 
@@ -1283,12 +1290,17 @@ def generate(self):
         if num_elements == None:
             num_elements = int(2*pi/mesh_size)
 
+        #Remaining height to be meshed
+        height = h - y                  
+        
+        #Number of times the plate height fits in the remaining unmeshed height, resulting at least one element
+        n_vert = max(int(abs(height)/mesh_size),1)
+        
+        #Element height in the vertical direction
+        h_y = height/n_vert    
+
         # Mesh the cylinder from the bottom toward the top
         while round(y, 10) < round(h, 10):
-            
-            height = h - y                  # Remaining height to be meshed
-            n_vert = int(height/mesh_size)  # Number of times the plate height fits in the remaining unmeshed height
-            h_y = height/n_vert             # Element height in the vertical direction
         
             # Create a mesh of nodes for the ring
             if self.axis == 'Y':

From 10307013875d410a56c9b35a1e1c79cb68de6009 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 5 Aug 2023 15:05:25 -0400
Subject: [PATCH 17/20] Moved `_renumber` out of `FEModel3D`

---
 PyNite/Analysis.py  | 35 ++++++++++++++++++++++++++++++++++-
 PyNite/FEModel3D.py | 36 +++---------------------------------
 2 files changed, 37 insertions(+), 34 deletions(-)

diff --git a/PyNite/Analysis.py b/PyNite/Analysis.py
index c3e3343b..32d64480 100644
--- a/PyNite/Analysis.py
+++ b/PyNite/Analysis.py
@@ -1,5 +1,8 @@
 from math import isclose
 
+def _prepare_model(model):
+    pass
+
 def _check_stability(model, K):
     """
     Identifies nodal instabilities in a model's stiffness matrix.
@@ -542,4 +545,34 @@ def _check_statics(model, combo_tags=None):
 
     # Print the static check table
     print(statics_table)
-    print('')
\ No newline at end of file
+    print('')
+
+def _renumber(model):
+    """
+    Assigns node and element ID numbers to be used internally by the program. Numbers are
+    assigned according to the order in which they occur in each dictionary.
+    """
+    
+    # Number each node in the model
+    for id, node in enumerate(model.Nodes.values()):
+        node.ID = id
+    
+    # Number each spring in the model
+    for id, spring in enumerate(model.Springs.values()):
+        spring.ID = id
+
+    # Descritize all the physical members and number each member in the model
+    id = 0
+    for phys_member in model.Members.values():
+        phys_member.descritize()
+        for member in phys_member.sub_members.values():
+            member.ID = id
+            id += 1
+    
+    # Number each plate in the model
+    for id, plate in enumerate(model.Plates.values()):
+        plate.ID = id
+    
+    # Number each quadrilateral in the model
+    for id, quad in enumerate(model.Quads.values()):
+        quad.ID = id
\ No newline at end of file
diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index 35c3d4e8..751fce19 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -1950,7 +1950,7 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
                 phys_member.active[combo_name] = True
         
         # Assign an internal ID to all nodes and elements in the model
-        self._renumber()
+        Analysis._renumber(self)
 
         # Get the auxiliary list used to determine how the matrices will be partitioned
         D1_indices, D2_indices, D2 = self._aux_list()
@@ -2137,7 +2137,7 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
                 phys_member.active[combo_name] = True
         
         # Assign an internal ID to all nodes and elements in the model
-        self._renumber()
+        Analysis._renumber(self)
 
         # Get the auxiliary list used to determine how the matrices will be partitioned
         D1_indices, D2_indices, D2 = self._aux_list()
@@ -2306,7 +2306,7 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
                 phys_member.active[combo_name] = True
                
         # Assign an internal ID to all nodes and elements in the model
-        self._renumber()
+        Analysis._renumber(self)
         
         # Get the auxiliary list used to determine how the matrices will be partitioned
         D1_indices, D2_indices, D2 = self._aux_list()
@@ -2530,36 +2530,6 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         
         # Flag the model as solved
         self.solution = 'P-Delta'
-       
-    def _renumber(self):
-        """
-        Assigns node and element ID numbers to be used internally by the program. Numbers are
-        assigned according to the order in which they occur in each dictionary.
-        """
-        
-        # Number each node in the model
-        for id, node in enumerate(self.Nodes.values()):
-            node.ID = id
-        
-        # Number each spring in the model
-        for id, spring in enumerate(self.Springs.values()):
-            spring.ID = id
-
-        # Descritize all the physical members and number each member in the model
-        id = 0
-        for phys_member in self.Members.values():
-            phys_member.descritize()
-            for member in phys_member.sub_members.values():
-                member.ID = id
-                id += 1
-        
-        # Number each plate in the model
-        for id, plate in enumerate(self.Plates.values()):
-            plate.ID = id
-        
-        # Number each quadrilateral in the model
-        for id, quad in enumerate(self.Quads.values()):
-            quad.ID = id
     
     def unique_name(self, dictionary, prefix):
         """Returns the next available unique name for a dictionary of objects.

From fbfa2b45f35f2d0634a73188de28ab0e65aad5d2 Mon Sep 17 00:00:00 2001
From: Craig Brinck <cbrinck@outlook.com>
Date: Sat, 5 Aug 2023 15:23:26 -0400
Subject: [PATCH 18/20] Analysis code simplification/organization

---
 PyNite/Analysis.py  | 30 +++++++++++++++++-
 PyNite/FEModel3D.py | 77 +++++----------------------------------------
 2 files changed, 37 insertions(+), 70 deletions(-)

diff --git a/PyNite/Analysis.py b/PyNite/Analysis.py
index 32d64480..e5fdeb59 100644
--- a/PyNite/Analysis.py
+++ b/PyNite/Analysis.py
@@ -1,7 +1,35 @@
 from math import isclose
+from PyNite.LoadCombo import LoadCombo
 
 def _prepare_model(model):
-    pass
+    """Prepares a model for analysis by ensuring at least one load combination is defined, generating all meshes that have not already been generated, activating all non-linear members, and internally numbering all nodes and elements.
+
+    :param model: The model being prepared for analysis.
+    :type model: FEModel3D
+    """
+
+    # Ensure there is at least 1 load combination to solve if the user didn't define any
+    if model.LoadCombos == {}:
+        # Create and add a default load combination to the dictionary of load combinations
+        model.LoadCombos['Combo 1'] = LoadCombo('Combo 1', factors={'Case 1':1.0})
+    
+    # Generate all meshes
+    for mesh in model.Meshes.values():
+        if mesh.is_generated == False:
+            mesh.generate()
+
+    # Activate all springs and members for all load combinations
+    for spring in model.Springs.values():
+        for combo_name in model.LoadCombos.keys():
+            spring.active[combo_name] = True
+    
+    # Activate all physical members for all load combinations
+    for phys_member in model.Members.values():
+        for combo_name in model.LoadCombos.keys():
+            phys_member.active[combo_name] = True
+    
+    # Assign an internal ID to all nodes and elements in the model. This number is different from the name used by the user to identify nodes and elements.
+    _renumber(model)
 
 def _check_stability(model, K):
     """
diff --git a/PyNite/FEModel3D.py b/PyNite/FEModel3D.py
index 751fce19..4753e983 100644
--- a/PyNite/FEModel3D.py
+++ b/PyNite/FEModel3D.py
@@ -1924,33 +1924,13 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
             print('+-----------+')
             print('| Analyzing |')
             print('+-----------+')
-
+        
         # Import `scipy` features if the sparse solver is being used
         if sparse == True:
             from scipy.sparse.linalg import spsolve
 
-        # Ensure there is at least 1 load combination to solve if the user didn't define any
-        if self.LoadCombos == {}:
-            # Create and add a default load combination to the dictionary of load combinations
-            self.LoadCombos['Combo 1'] = LoadCombo('Combo 1', factors={'Case 1':1.0})
-        
-        # Generate all meshes
-        for mesh in self.Meshes.values():
-            if mesh.is_generated == False:
-                mesh.generate()
-
-        # Activate all springs and members for all load combinations
-        for spring in self.Springs.values():
-            for combo_name in self.LoadCombos.keys():
-                spring.active[combo_name] = True
-        
-        # Activate all physical members for all load combinations
-        for phys_member in self.Members.values():
-            for combo_name in self.LoadCombos.keys():
-                phys_member.active[combo_name] = True
-        
-        # Assign an internal ID to all nodes and elements in the model
-        Analysis._renumber(self)
+        # Prepare the model for analysis
+        Analysis._prepare_model(self)
 
         # Get the auxiliary list used to determine how the matrices will be partitioned
         D1_indices, D2_indices, D2 = self._aux_list()
@@ -2111,33 +2091,13 @@ def analyze_linear(self, log=False, check_stability=True, check_statics=False, s
             print('+-------------------+')
             print('| Analyzing: Linear |')
             print('+-------------------+')
-
+        
         # Import `scipy` features if the sparse solver is being used
         if sparse == True:
             from scipy.sparse.linalg import spsolve
 
-        # Ensure there is at least 1 load combination to solve if the user didn't define any
-        if self.LoadCombos == {}:
-            # Create and add a default load combination to the dictionary of load combinations
-            self.LoadCombos['Combo 1'] = LoadCombo('Combo 1', factors={'Case 1':1.0})
-                
-        # Generate all meshes
-        for mesh in self.Meshes.values():
-            if mesh.is_generated == False:
-                mesh.generate()
-        
-        # Activate all springs for all load combinations
-        for spring in self.Springs.values():
-            for combo_name in self.LoadCombos.keys():
-                spring.active[combo_name] = True
-        
-        # Activate all physical members for all load combinations
-        for phys_member in self.Members.values():
-            for combo_name in self.LoadCombos.keys():
-                phys_member.active[combo_name] = True
-        
-        # Assign an internal ID to all nodes and elements in the model
-        Analysis._renumber(self)
+        # Prepare the model for analysis
+        Analysis._prepare_model(self)
 
         # Get the auxiliary list used to determine how the matrices will be partitioned
         D1_indices, D2_indices, D2 = self._aux_list()
@@ -2284,29 +2244,8 @@ def analyze_PDelta(self, log=False, check_stability=True, max_iter=30, tol=0.01,
         if sparse == True:
             from scipy.sparse.linalg import spsolve
 
-        # Ensure there is at least 1 load combination to solve if the user didn't define any
-        if self.LoadCombos == {}:
-            # Create and add a default load combination to the dictionary of load combinations
-            self.LoadCombos['Combo 1'] = LoadCombo('Combo 1', factors={'Case 1':1.0})
-                
-        # Generate all meshes
-        for mesh in self.Meshes.values():
-            if mesh.is_generated == False:
-                mesh.generate()
-        
-        # Activate all springs for all load combinations. They can be turned inactive
-        # during the course of the tension/compression-only analysis
-        for spring in self.Springs.values():
-            for combo_name in self.LoadCombos.keys():
-                spring.active[combo_name] = True
-                
-        # Activate all physical members for all load combinations
-        for phys_member in self.Members.values():
-            for combo_name in self.LoadCombos.keys():
-                phys_member.active[combo_name] = True
-               
-        # Assign an internal ID to all nodes and elements in the model
-        Analysis._renumber(self)
+        # Prepare the model for analysis
+        Analysis._prepare_model(self)
         
         # Get the auxiliary list used to determine how the matrices will be partitioned
         D1_indices, D2_indices, D2 = self._aux_list()

From 721ee8362d926d278fc0f0e554cf5c28094ae8bd Mon Sep 17 00:00:00 2001
From: Kevin Russell <kevin@ottermatics.com>
Date: Mon, 7 Aug 2023 09:15:52 -0700
Subject: [PATCH 19/20] Revert mesh loop changes, while accomodating at least
 one mesh element.

---
 PyNite/Mesh.py | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/PyNite/Mesh.py b/PyNite/Mesh.py
index 1d737b22..77dfd2e5 100644
--- a/PyNite/Mesh.py
+++ b/PyNite/Mesh.py
@@ -1290,18 +1290,13 @@ def generate(self):
         if num_elements == None:
             num_elements = int(2*pi/mesh_size)
 
-        #Remaining height to be meshed
-        height = h - y                  
-        
-        #Number of times the plate height fits in the remaining unmeshed height, resulting at least one element
-        n_vert = max(int(abs(height)/mesh_size),1)
-        
-        #Element height in the vertical direction
-        h_y = height/n_vert    
-
         # Mesh the cylinder from the bottom toward the top
         while round(y, 10) < round(h, 10):
-        
+
+            height = h - y  #Remaining height to be meshed
+            # Number of times the plate height fits in the remaining unmeshed height, resulting at least one element
+            n_vert = max(int(abs(height)/mesh_size),1)
+            h_y = height/n_vert             # Element height in the vertical direction
             # Create a mesh of nodes for the ring
             if self.axis == 'Y':
                 ring = CylinderRingMesh(radius, h_y, num_elements, thickness, material, self.model, 1, 1, [0, y, 0],

From 23f5d3c328df40f89eea6f935965670b936b1ccf Mon Sep 17 00:00:00 2001
From: SoundsSerious <keavin7@ottermatics.com>
Date: Thu, 10 Aug 2023 21:10:57 -0700
Subject: [PATCH 20/20] uncomment AnnulusMesh.generate()

---
 PyNite/Mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/PyNite/Mesh.py b/PyNite/Mesh.py
index ced6af37..6d10632f 100644
--- a/PyNite/Mesh.py
+++ b/PyNite/Mesh.py
@@ -780,7 +780,7 @@ def __init__(self, mesh_size, outer_radius, inner_radius, thickness, material, m
         self.num_quads_inner = None
         self.num_quads_outer = None
         
-        # self.generate()
+        self.generate()
     
     def generate(self):