Merge pull request #215 from JWock82/main

Preparing for a new release

PyNite/Analysis.py

@@ -392,47 +392,16 @@ def _unpartition_disp(model, D1, D2, D1_indices, D2_indices):
     # Step through each node in the model
     for node in model.nodes.values():
 
-        if node.ID*6 + 0 in D2_indices:
-            # Get the enforced displacement
-            D[(node.ID*6 + 0, 0)] = D2[D2_indices.index(node.ID*6 + 0), 0]
-        else:
-            # Get the calculated displacement
-            D[(node.ID*6 + 0, 0)] = D1[D1_indices.index(node.ID*6 + 0), 0]
-
-        if node.ID*6 + 1 in D2_indices:
-            # Get the enforced displacement
-            D[(node.ID*6 + 1, 0)] = D2[D2_indices.index(node.ID*6 + 1), 0]
-        else:
-            # Get the calculated displacement
-            D[(node.ID*6 + 1, 0)] = D1[D1_indices.index(node.ID*6 + 1), 0]
-
-        if node.ID*6 + 2 in D2_indices:
-            # Get the enforced displacement
-            D[(node.ID*6 + 2, 0)] = D2[D2_indices.index(node.ID*6 + 2), 0]
-        else:
-            # Get the calculated displacement
-            D[(node.ID*6 + 2, 0)] = D1[D1_indices.index(node.ID*6 + 2), 0]
-
-        if node.ID*6 + 3 in D2_indices:
-            # Get the enforced rotation
-            D[(node.ID*6 + 3, 0)] = D2[D2_indices.index(node.ID*6 + 3), 0]
-        else:
-            # Get the calculated rotation
-            D[(node.ID*6 + 3, 0)] = D1[D1_indices.index(node.ID*6 + 3), 0]
+        # Step through each degree of freedom at the node
+        for i in range(6):
 
-        if node.ID*6 + 4 in D2_indices:
-            # Get the enforced rotation
-            D[(node.ID*6 + 4, 0)] = D2[D2_indices.index(node.ID*6 + 4), 0]
-        else:
-            # Get the calculated rotation
-            D[(node.ID*6 + 4, 0)] = D1[D1_indices.index(node.ID*6 + 4), 0]
-
-        if node.ID*6 + 5 in D2_indices:
-            # Get the enforced rotation
-            D[(node.ID*6 + 5, 0)] = D2[D2_indices.index(node.ID*6 + 5), 0]
-        else:
-            # Get the calculated rotation
-            D[(node.ID*6 + 5, 0)] = D1[D1_indices.index(node.ID*6 + 5), 0]
+            # Check if the dof is in the list of enforced displacements
+            if node.ID*6 + i in D2_indices:
+                # Get the enforced displacement
+                D[(node.ID*6 + i, 0)] = D2[D2_indices.index(node.ID*6 + i), 0]
+            else:
+                # Get the calculated displacement
+                D[(node.ID*6 + i, 0)] = D1[D1_indices.index(node.ID*6 + i), 0]
 
     # Return the displacement vector
     return D
@@ -503,110 +472,52 @@ def _sum_displacements(model, Delta_D1, Delta_D2, D1_indices, D2_indices, combo)
         node.RY[combo.name] += Delta_D[node.ID*6 + 4, 0]
         node.RZ[combo.name] += Delta_D[node.ID*6 + 5, 0]
 
-def _check_TC_convergence(model, combo_name='Combo 1', log=True):
-    
+def _check_TC_convergence(model, combo_name="Combo 1", log=True, spring_tolerance=0, member_tolerance=0):
+
     # 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')
+    if log:
+        print("- Checking for tension/compression-only support spring convergence")
+    # Loop through each node and each directional spring to check and update their active status
     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
-
+        for direction in ["DX", "DY", "DZ", "RX", "RY", "RZ"]:
+            spring = getattr(node, f"spring_{direction}")
+            displacement = getattr(node, direction)[combo_name]
+
+            if spring[1] is not None:
+                # Determine if the spring should be active based on its designation ('+' or '-') and the displacement
+                should_be_active = (
+                    spring[1] == "-" and displacement <= -spring_tolerance
+                ) or (spring[1] == "+" and displacement >= spring_tolerance)
+
+                # Check if there's a need to switch the active state of the spring
+                if spring[2] != should_be_active:
+                    spring[2] = should_be_active
+                    convergence = False
+
     # 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:
+            if (
+                spring.tension_only == True
+                and spring.axial(combo_name) > spring_tolerance
+            ):
                 spring.active[combo_name] = False
                 convergence = False
-            
+
             # Check if compression-only conditions exist
-            elif spring.comp_only == True and spring.axial(combo_name) < 0:
+            elif (
+                spring.comp_only == True
+                and spring.axial(combo_name) < -spring_tolerance
+            ):
                 spring.active[combo_name] = False
                 convergence = False
 
@@ -617,19 +528,43 @@ def _check_TC_convergence(model, combo_name='Combo 1', log=True):
         # 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:
+            # Check if a tension-only conditions exist
+            if (
+                phys_member.tension_only == True
+                and phys_member.max_axial(combo_name) > member_tolerance
+            ):
+                # Deactivate the physical member
                 phys_member.active[combo_name] = False
+
+                # Deactivate all the sub-members
+                for sub_member in phys_member.sub_members.values():
+                    sub_member.active[combo_name] = False
+
+                # Flag the analysis as not converged
                 convergence = False
 
-            # Check if compression-only conditions exist
-            elif phys_member.comp_only == True and phys_member.min_axial(combo_name) < 0:
+            # Check if a compression-only conditions exist
+            elif (
+                phys_member.comp_only == True
+                and phys_member.min_axial(combo_name) < -member_tolerance
+            ):
+                # Deactivate the physical member
                 phys_member.active[combo_name] = False
+
+                # Deactivate all the sub-members
+                for sub_member in phys_member.sub_members.values():
+                    sub_member.active[combo_name] = False
+
+                # Flag the analysis as not converged
                 convergence = False
 
+        # Reset the sub-member's flag to unsolved. This will allow it to resolve for the same load combination after subsequent iterations have made further changes.
+        for sub_member in phys_member.sub_members.values():
+            sub_member._solved_combo = None
+
     # Return whether the TC analysis has converged
     return convergence
-
+  
 def _calc_reactions(model, log=False, combo_tags=None):
     """
     Calculates reactions internally once the model is solved.
@@ -1001,7 +936,7 @@ def _partition_D(model):
 
     D1_indices = [] # A list of the indices for the unknown nodal displacements
     D2_indices = [] # A list of the indices for the known nodal displacements
-    D2 = []         # A list of the values of the known nodal displacements (D != None)
+    D2 = []         # A list of the values of the known nodal displacements
 
     # Create the auxiliary table
     for node in model.nodes.values():

PyNite/BeamSegY.py

@@ -20,14 +20,14 @@ def moment(self, x, P_delta=False):
         P1 = self.P1
         w1 = self.w1
         w2 = self.w2
-        delta1 = self.delta1
-        delta = self.deflection(x)
         L = self.Length()
 
         # M = M1 + V1*x + w1*x**2/2 + x**3*(-w1 + w2)/(6*L)
         M = -M1 - V1*x - w1*x**2/2 - x**3*(-w1 + w2)/(6*L)
 
         if P_delta == True:
+            delta1 = self.delta1
+            delta = self.deflection(x)
             M += P1*(delta - delta1)
 
         return M
@@ -49,12 +49,13 @@ def slope(self, x, P_delta=False):
         w1 = self.w1
         w2 = self.w2
         theta_1 = self.theta1
-        delta_x = self.deflection(x, P_delta)
-        delta_1 = self.delta1
+
         L = self.Length()
         EI = self.EI
 
         if P_delta == True:
+            delta_x = self.deflection(x, P_delta)
+            delta_1 = self.delta1
             return theta_1 + (-V1*x**2/2 - w1*x**3/6 + x*(-M1 - P1*delta_1 + P1*delta_x) + x**4*(w1 - w2)/(24*L))/EI
         else:
             return theta_1 + (-V1*x**2/2 - w1*x**3/6 + x*(-M1) + x**4*(w1 - w2)/(24*L))/EI

PyNite/BeamSegZ.py

@@ -4,6 +4,8 @@
 
 @author: D. Craig Brinck, SE
 """
+from numpy import full, array
+
 # %%
 # A mathematically continuous beam segment
 class BeamSegZ():
@@ -103,17 +105,16 @@ def moment(self, x, P_delta=False):
         V1 = self.V1
         M1 = self.M1
         P1 = self.P1
-        Px = self.axial(x)
         w1 = self.w1
         w2 = self.w2
-        delta_1 = self.delta1
-        delta_x = self.deflection(x)
         L = self.Length()
 
         M = M1 - V1*x - w1*x**2/2 - x**3*(-w1 + w2)/(6*L)
 
         # Include the P-little-delta moment if a P-Delta analysis was run
         if P_delta == True:
+            delta_1 = self.delta1
+            delta_x = self.deflection(x)
             M += P1*(delta_x - delta_1)
 
         # Return the computed moment
@@ -129,14 +130,20 @@ def axial(self, x):
 
         return P1 + (p2 - p1)/(2*L)*x**2 + p1*x
 
-    def Torsion(self):
+    def Torsion(self, x = 0):
         """
         Returns the torsional moment in the segment.
         """
 
         # The torsional moment is constant across the segment
         # This can be updated in the future for distributed torsional forces
-        return self.T1
+
+        #As the return value is not calculated as a function of x (for now), we need to check
+        #whether x is an array, and if so, return a results array of the same length
+        if isinstance(x, (int, float)):
+            return self.T1
+        else:
+            return full(len(x), self.T1)
 
     def slope(self, x, P_delta=False):
         """Returns the slope of the elastic curve at any point `x` along the segment.
@@ -154,8 +161,6 @@ def slope(self, x, P_delta=False):
         P1 = self.P1
         w1 = self.w1
         w2 = self.w2
-        delta_1 = self.delta1
-        delta_x = self.deflection(x, P_delta)
         theta_1 = self.theta1
         L = self.Length()
         EI = self.EI
@@ -164,6 +169,8 @@ def slope(self, x, P_delta=False):
         theta_x = theta_1 - (-V1*x**2/2 - w1*x**3/6 + x*M1 + x**4*(w1 - w2)/(24*L))/EI
 
         if P_delta == True:
+            delta_1 = self.delta1
+            delta_x = self.deflection(x, P_delta)
             theta_x -= x*(-P1*delta_1 + P1*delta_x)/EI
 
         # TODO: This is an old equation left for reference. Delete it after the new equation has been proved over time.

PyNite/FEModel3D.py

@@ -1914,7 +1914,7 @@ def D(self, combo_name='Combo 1'):
         # Return the global displacement vector
         return self._D[combo_name]
 
-    def analyze(self, log=False, check_stability=True, check_statics=False, max_iter=30, sparse=True, combo_tags=None):
+    def analyze(self, log=False, check_stability=True, check_statics=False, max_iter=30, sparse=True, combo_tags=None, spring_tolerance=0, member_tolerance=0):
         """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.
@@ -2006,9 +2006,10 @@ def analyze(self, log=False, check_stability=True, check_statics=False, max_iter
                 Analysis._store_displacements(self, D1, D2, D1_indices, D2_indices, combo)
 
                 # Check for tension/compression-only convergence
-                convergence = Analysis._check_TC_convergence(self, combo.name, log=log)
+                convergence = Analysis._check_TC_convergence(self, combo.name, log=log, spring_tolerance=spring_tolerance, member_tolerance=member_tolerance)
 
                 if convergence == False:
+
                     if log: print('- Tension/compression-only analysis did not converge. Adjusting stiffness matrix and reanalyzing.')
                 else:
                     if log: print('- Tension/compression-only analysis converged after ' + str(iter_count) + ' iteration(s)')
@@ -2027,7 +2028,7 @@ 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, combo_tags)
-        
+
         # Flag the model as solved
         self.solution = 'Linear TC'