Deflection shape fixes for all element types (#77)

src/pycba/analysis.py

@@ -223,7 +223,7 @@ def _forces(self) -> np.ndarray:
 
         for i in range(self._n):
             dof_i = 2 * i
-            fmbr = self._beam.get_cnl(i)
+            fmbr = self._beam.get_ref(i)
             # Cumulatively apply forces in opposite direction
             f[dof_i : dof_i + 4] -= fmbr
         return f

src/pycba/beam.py

@@ -280,9 +280,10 @@ def get_local_span_coords(self, pos: float) -> (int, float):
 
         return ispan, pos_in_span
 
-    def get_cnl(self, i_span: int) -> LoadCNL:
+    def get_ref(self, i_span: int) -> LoadCNL:
         """
-        Returns Consistent Nodal Loads for the member
+        Returns Released End Forces for the member; that is, the Consistent Nodal Loads
+        modified for the element type (i.e. releases)
 
         Parameters
         ----------
@@ -291,18 +292,18 @@ def get_cnl(self, i_span: int) -> LoadCNL:
 
         Returns
         -------
-        cnl : LoadCNL
+        ref : LoadCNL
             The totalled CNL object for the member, considering all loads.
 
         """
-        cnl = np.zeros(4)
+        ref = np.zeros(4)
         L = self.mbr_lengths[i_span]
         eType = self.mbr_eletype[i_span]
 
         for load in self._loads:
             if load.i_span == i_span:
-                cnl += load.get_cnl(L, eType)
-        return cnl
+                ref += load.get_ref(L, eType)
+        return ref
 
     def get_span_k(self, i_span: int) -> np.ndarray:
         """

src/pycba/inf_lines.py

@@ -118,45 +118,67 @@ def get_il(self, poi: float, load_effect: str) -> (np.ndarray, np.ndarray):
         if not self.vResults:
             self.create_ils()
 
-        x = self.vResults[0].results.x
-        dx = x[2] - x[1]
-        idx = np.where(np.abs(x - poi) <= dx * 1e-6)[0][0]
-        #idx = np.abs(x - poi).argmin()
-
+        x = self.vResults[0].results.x   
         npts = len(self.vResults)
         eta = np.zeros(npts)
 
-        #
-        # Getting the correct vertical reaction is tricky
-        # Should be relatively easy to extend to get moment reactions too.
+        # Preparations for reaction ILs
         #
         # Get vector of the node locations
         node_locations = np.cumsum(np.insert(self.ba.beam.mbr_lengths, 0, 0))
-        # The indices of the supported vertical DOFs wrt the node locations vector
-        vert_sup_dof_idx = np.where(np.array(self.ba._beam.restraints)[::2]==-1)[0]
-        # The locations then of these vertical supports
-        vert_sup_locs = node_locations[vert_sup_dof_idx]
-        # The index of the closest vertical support
-        closest_vert_sup_idx = np.abs(vert_sup_locs-poi).argmin()
-        # And its value
-        closest_vert_sup = vert_sup_locs[closest_vert_sup_idx]
-        # And now the indixe of this support in the node locations vector
-        node_idx = np.where(node_locations==closest_vert_sup)[0][0]
-        # And hence its index in the overall DOFs vector
-        dof_idx = 2*node_idx
-
-        # Now we must link the supported DOF to the index in the BeamAnalysis reactions vector
+        # Link the supported DOF to the index in the BeamAnalysis reactions vector
         idx_mask = np.zeros_like(self.ba._beam.restraints)
         idx_mask[np.where(np.array(self.ba._beam.restraints)==-1)] = np.arange(self.ba.beam.no_fixed_restraints)
-        # And finally the index of the vertical support nearest the POI in the reactions vector
-        vert_sup_idx = idx_mask[dof_idx]
 
-        for i, res in enumerate(self.vResults):
-            if load_effect == "V":
+        #idx = np.abs(x - poi).argmin()
+
+        if load_effect.upper() == "V":
+            dx = x[2] - x[1]
+            idx = np.where(np.abs(x - poi) <= dx * 1e-6)[0][0]
+            for i, res in enumerate(self.vResults):
                 eta[i] = res.results.V[idx]
-            elif load_effect == "R":
+
+        elif load_effect.upper() == "R":
+            #
+            # Getting the correct reaction is tricky
+            #
+            # The indices of the supported DOFs wrt the node locations vector
+            vert_sups_dof_idx = np.where(np.array(self.ba._beam.restraints)[::2]==-1)[0]
+            # The locations then of these supports
+            vert_sups_locs = node_locations[vert_sups_dof_idx]
+            # The index of the closest support
+            closest_vert_sup_idx = np.abs(vert_sups_locs-poi).argmin()
+            # And its value
+            closest_vert_sup = vert_sups_locs[closest_vert_sup_idx]
+            # And now the index of this support in the node locations vector
+            vert_sup_node_idx = np.where(node_locations==closest_vert_sup)[0][0]
+            # And hence its index in the overall DOFs vector
+            vert_sup_dof_idx = 2*vert_sup_node_idx
+            # And finally the index of the support nearest the POI in the reactions vector
+            vert_sup_idx = idx_mask[vert_sup_dof_idx]            
+
+            for i, res in enumerate(self.vResults):
                 eta[i] = res.R[vert_sup_idx]
-            else:
+
+        elif load_effect.upper() == "MR":
+            #
+            # Follows the same logic for the vertical reaction
+            #
+            mt_sups_dof_idx = np.where(np.array(self.ba._beam.restraints)[1::2]==-1)[0]
+            mt_sups_locs = node_locations[mt_sups_dof_idx]
+            closest_mt_sup_idx = np.abs(mt_sups_locs-poi).argmin()
+            closest_mt_sup = mt_sups_locs[closest_mt_sup_idx]
+            mt_sup_node_idx = np.where(node_locations==closest_mt_sup)[0][0]
+            mt_sup_dof_idx = 2*mt_sup_node_idx+1
+            mt_sup_idx = idx_mask[mt_sup_dof_idx]
+
+            for i, res in enumerate(self.vResults):
+                eta[i] = res.R[mt_sup_idx]
+
+        else:
+            dx = x[2] - x[1]
+            idx = np.where(np.abs(x - poi) <= dx * 1e-6)[0][0]
+            for i, res in enumerate(self.vResults):
                 eta[i] = res.results.M[idx]
 
         return (np.array(self.pos), eta)

src/pycba/load.py

@@ -182,7 +182,7 @@ def __init__(self, i_span: int):
         """
         self.i_span = i_span
 
-    def get_cnl(self, L, eType):
+    def get_cnl(self, L):
         # Enforce virtual base class
         raise NotImplementedError
 
@@ -215,29 +215,27 @@ def H(self, v: np.ndarray, value: float = 0.0) -> np.ndarray:
         """
         return np.heaviside(v, value)
 
-    def released_end_forces(self, cnl: LoadCNL, L: float, eType: int) -> LoadCNL:
+
+    def get_ref(self, L: float, eType: int) -> LoadCNL:
         """
-        The released end forces for each element type: converts the Consistent Nodal
-        Loads of the applied loading to the correct nodal loading depending on the
-        element type.
+        Returns the Released End Forces for a span of length L of element eType: 
+        converts the Consistent Nodal Loads of the applied loading to the correct nodal 
+        loading depending on the element type.
 
         Parameters
         ----------
-        cnl : LoadCNL
-            The nodal loading statically consistent with the externally applied loads.
         L : float
-            The length of the member.
+            The length of the member
         eType : int
-            The element type.
+            The member element type
 
         Returns
         -------
         LoadCNL
-            The nodal loads to be applied in the analysis, consistent with the element
-            type.
-
+            Released End Forces for this load type: the nodal loads to be applied in 
+            the analysis, consistent with the element type.
         """
-
+        cnl = self.get_cnl(L, eType)
         ref = np.zeros(4)
         fm = 6 / (4 * L)  # flexibility coeff for moment
 
@@ -252,10 +250,10 @@ def released_end_forces(self, cnl: LoadCNL, L: float, eType: int) -> LoadCNL:
             ref[2] = -fm * cnl.Ma
             ref[3] = 0.5 * cnl.Ma
         elif eType == 4:  # keep only vertical, remove moments
-            ref[0] = 0
-            ref[1] = 1.0 * cnl.Va
-            ref[2] = 0
-            ref[3] = 1.0 * cnl.Va
+            ref[0] = -(cnl.Ma+cnl.Mb)/L
+            ref[1] = 1.0 * cnl.Ma
+            ref[2] = (cnl.Ma+cnl.Mb)/L
+            ref[3] = 1.0 * cnl.Mb
         else:
             # no nothing if it is FF
             pass
@@ -319,8 +317,7 @@ def get_cnl(self, L: float, eType: int) -> LoadCNL:
             Ma=w * L**2 / 12.0,
             Mb=-w * L**2 / 12.0,
         )
-
-        return self.released_end_forces(cnl, L, eType)
+        return cnl
 
     def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults:
         """
@@ -416,8 +413,7 @@ def get_cnl(self, L, eType) -> LoadCNL:
             Ma=P * a * b**2 / L**2,
             Mb=-P * a**2 * b / L**2,
         )
-        # implicit conversion to tuple in correct order
-        return self.released_end_forces(cnl, L, eType)
+        return cnl
 
     def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults:
         """
@@ -512,7 +508,7 @@ def get_cnl(self, L, eType) -> LoadCNL:
             Mb=-(w * c / L**2) * (t * s**2 + (t - 2 * s) * c**2 / 12),
         )
         # implicit conversion to tuple in correct order
-        return self.released_end_forces(cnl, L, eType)
+        return cnl
 
     def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults:
         """
@@ -608,8 +604,7 @@ def get_cnl(self, L, eType) -> LoadCNL:
             Ma=Ma,
             Mb=Mb,
         )
-        # implicit conversion to tuple in correct order
-        return self.released_end_forces(cnl, L, eType)
+        return cnl
 
     def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults:
         """
@@ -690,8 +685,7 @@ def get_cnl(self, L, eType) -> LoadCNL:
             Ma=(m * b / L**2) * (2 * a - b),
             Mb=(m * a / L**2) * (2 * b - a),
         )
-        # implicit conversion to tuple in correct order
-        return self.released_end_forces(cnl, L, eType)
+        return cnl
 
     def get_mbr_results(self, x: np.ndarray, L: float) -> MemberResults:
         """

src/pycba/results.py

@@ -8,7 +8,7 @@
 import matplotlib.pyplot as plt
 from scipy import integrate
 from .beam import Beam
-from .load import MemberResults, LoadMaMb
+from .load import MemberResults, LoadMaMb, LoadCNL
 from copy import deepcopy
 
 
@@ -105,7 +105,7 @@ def _member_analysis(self, beam: Beam, d: np.ndarray) -> List[MemberResults]:
             fmbr = np.zeros(4)
             for j in range(4):
                 fmbr[j] = np.sum(kb[j][:] * dmbr[:])
-            fmbr += beam.get_cnl(i)
+            fmbr += beam.get_ref(i)
             res = self._member_values(beam, i, fmbr, dmbr)
             # Shift x vals by location of mbr starting point
             res.x += sumL
@@ -138,6 +138,9 @@ def _member_values(
         """
 
         L = beam.mbr_lengths[i_span]
+        EI = beam.mbr_EIs[i_span]
+        etype = beam.mbr_eletype[i_span]
+
         dx = L / self.npts
         x = np.zeros(self.npts + 3)
         x[1 : self.npts + 2] = dx * np.arange(0, self.npts + 1)
@@ -148,16 +151,31 @@ def _member_values(
         res = MaMb.get_mbr_results(x, L)
 
         # Now get the results for all the applied loads on a simple span
+        Ma=0
+        Mb=0
         for load in beam._loads:
             if load.i_span != i_span:
                 continue
             res += load.get_mbr_results(x, L)
+            cnl = load.get_cnl(L, etype)
+            Ma += cnl.Ma
+            Mb += cnl.Mb
+
+        # Check element type for any released displacements        
+        R0 = d[1]
+        theta = 0
+        phi_i = 0
+        # Account for end release if joint not already rotating
+        if etype != 1 and abs(R0) < 1e-6 :
+            theta = (d[2]-d[0])/L
+            phi_i = (L/(3*EI))*(-(f[1]-0.5*f[3])+(Ma-0.5*Mb))
+
+        R0 -= phi_i - theta
 
         # And superimpose end displacements using Moment-Area
-        h = L / self.npts
-        EI = beam.mbr_EIs[i_span]
+        h = L / self.npts        
 
-        R = integrate.cumtrapz(res.M[1:-1], dx=h, initial=0) / EI + d[1]
+        R = integrate.cumtrapz(res.M[1:-1], dx=h, initial=0) / EI + R0
         D = integrate.cumtrapz(R, dx=h, initial=0) + d[0]
 
         res.R[1:-1] = R