v0.8.7

doc/source/CHANGELOG.rst

@@ -1,6 +1,10 @@
 Changelog
 =============
 
+v0.8.7
+---------
+- Updated the class :class:`opstool.analysis.MomentCurvature` to moment-curvature analysis of fiber section.
+
 v0.8.6
 ---------
 - Fixed ``get_frame_props`` bugs in :py:class:`opstool.preprocessing.SecMesh`

doc/source/src/api/ana.MomentCurvature.rst

@@ -122,12 +122,12 @@ Extract limit state points based on fiber strain thresholds or other criteria.
    # The concrete fiber in the confined area reaches the ultimate compressive strain 0.0144
    phiu, Mu = mc.get_limit_state(matTag=matTagCCore,
                                  threshold=-0.0144,
-                                 use_peak_drop20=False
+                                 peak_drop=False
                                  )
-   # or use_peak_drop20
+   # or use peak_drop 
    # phiu, Mu = mc.get_limit_state(matTag=matTagCCore,
    #                               threshold=-0.0144,
-   #                               use_peak_drop20=True
+   #                               peak_drop=0.2
    #                              )
 
 Equivalent linearization according to area:

src/opstool/__about__.py

@@ -1 +1 @@
-__version__ = "0.8.6"
+__version__ = "0.8.7"

src/opstool/analysis/_sec_analysis.py

@@ -1,6 +1,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import openseespy.opensees as ops
+from typing import Union
 
 from ._smart_analyze import SmartAnalyze
 
@@ -22,12 +23,13 @@ def __init__(self, sec_tag: int, axial_force: float = 0) -> None:
         self.phi, self.M, self.FiberData = None, None, None
 
     def analyze(
-        self,
-        axis: str = "y",
-        max_phi: float = 1.0,
-        incr_phi: float = 1e-4,
-        post_peak_ratio: float = 0.8,
-        smart_analyze: bool = True,
+            self,
+            axis: str = "y",
+            max_phi: float = 0.25,
+            incr_phi: float = 1e-4,
+            limit_peak_ratio: float = 0.8,
+            smart_analyze: bool = True,
+            debug: bool = False
     ):
         """Performing Moment-Curvature Analysis.
 
@@ -39,12 +41,13 @@ def analyze(
             The maximum curvature to analyze, by default 1.0.
         incr_phi : float, optional
             Curvature analysis increment, by default 1e-4.
-        post_peak_ratio : float, optional
-            A ratio of the moment intensity after the peak used to stop the analysis., by default 0.7,
-            i.e. a 30% drop after peak.
+        limit_peak_ratio : float, optional
+            A ratio of the moment intensity after the peak used to stop the analysis., by default 0.8,
+            i.e. a 20% drop after peak.
         smart_analyze : bool, optional
             Whether to use smart analysis options, by default True.
-
+        debug: bool, optional
+            Whether to use debug mode when smart analysis is True, by default False.
         .. note::
             The termination of the analysis depends on whichever reaches `max_phi` or `post_peak_ratio` first.
         """
@@ -54,31 +57,52 @@ def analyze(
             axis=axis,
             max_phi=max_phi,
             incr_phi=incr_phi,
-            stop_ratio=post_peak_ratio,
+            stop_ratio=limit_peak_ratio,
             smart_analyze=smart_analyze,
+            debug=debug
         )
 
-    def plot_M_phi(self):
-        """Plot the moment-curvature relationship."""
+    def plot_M_phi(self, return_ax: bool = False):
+        """Plot the moment-curvature relationship.
+
+        Parameters
+        ------------
+        return_ax: bool, default=False
+            If True, return the axes for the plot of matplotlib.
+
+        """
         _, ax = plt.subplots(1, 1, figsize=(10, 10 * 0.618))
         ax.plot(self.phi, self.M, lw=3, c="blue")
         ax.set_title(
             "$M-\\phi$",
             fontsize=28,
         )
-        ax.set_xlabel("$\phi$", fontsize=25)
+        ax.set_xlabel("$\\phi$", fontsize=25)
         ax.set_ylabel("$M$", fontsize=25)
         plt.xticks(fontsize=15)
         plt.yticks(fontsize=15)
-        ax.spines["bottom"].set_linewidth(1.2)
-        ax.spines["left"].set_linewidth(1.2)
-        ax.spines["right"].set_linewidth(1.2)
-        ax.spines["top"].set_linewidth(1.2)
+        for loc in ["bottom", "left", "right", "top"]:
+            ax.spines[loc].set_linewidth(1.0)
         plt.gcf().subplots_adjust(bottom=0.15)
-        plt.show()
+        if return_ax:
+            return ax
+        else:
+            plt.show()
+
+    def plot_fiber_responses(
+            self,
+            return_ax: bool = False
+
+    ):
+        """Plot the stress-strain histories of fiber by loc and matTag.
+        Plot the fiber response of the material Tag ``mat`` closest to the ``loc``.
+
+        Parameters
+        -----------
+        return_ax: bool, default=False
+            If True, return the axes for the plot of matplotlib.
+        """
 
-    def plot_fiber_responses(self):
-        """Plot the stress-strain histories for all fibers of each material."""
         fiber_data = self.FiberData
         matTags = np.unique(fiber_data[-1][:, 3])
 
@@ -96,10 +120,13 @@ def plot_fiber_responses(self):
             ax.set_title(f"matTag = {mat:.0f}", fontsize=15)
             ax.tick_params(labelsize=12)
             ax.set_ylabel("stress", fontsize=16)
-        ax.set_xlabel("strain", fontsize=16)
+            ax.set_xlabel("strain", fontsize=16)
         plt.subplots_adjust(wspace=0.02, hspace=0.4)
         plt.tight_layout()
-        plt.show()
+        if return_ax:
+            return axs
+        else:
+            plt.show()
 
     def get_phi(self):
         """Get the curvature array.
@@ -157,13 +184,13 @@ def get_fiber_data(self):
         Returns
         -------
         Shape-(n,m,6) Array.
-            n is the length of moment and curvature array, m is the fiber number,
+            n is the length of analysis steps, m is the fiber number,
             6 contain ("yCoord", "zCoord", "area", 'mat', "stress", "strain")
         """
         return self.FiberData
 
     def get_limit_state(
-        self, matTag: int = 1, threshold: float = 0, use_peak_drop20: bool = False
+            self, matTag: int = 1, threshold: float = 0, peak_drop: Union[float, bool] = False
     ):
         """Get the curvature and moment corresponding to a certain limit state.
 
@@ -173,9 +200,11 @@ def get_limit_state(
             The OpenSeesPy material Tag used to determine the limit state., by default 1
         threshold : float, optional
             The strain threshold used to determine the limit state by material `matTag`, by default 0
-        use_peak_drop20 : bool, optional
-            If True, A 20% drop from the peak value of the moment will be used as the limit state,
-            and `matTag` and `threshold` are not needed, by default False.
+        peak_drop : Union[float, bool], optional
+            If True, A default 20% drop from the peak value of the moment will be used as the limit state;
+            If float in [0, 1], this means that the ratio of ultimate strength to peak value is
+            specified by this value, for example, peak_drop = 0.3, the ultimate strength = 0.7 * peak.
+            `matTag` and `threshold` are not needed, by default False.
 
         Returns
         -------
@@ -185,12 +214,16 @@ def get_limit_state(
         phi = self.phi
         M = self.M
         fiber_data = self.FiberData
-        if use_peak_drop20:
+        if peak_drop:
+            if peak_drop is True:
+                ratio_ = 0.8
+            else:
+                ratio_ = 1 - peak_drop
             idx = np.argmax(M)
-            au = np.argwhere(M[idx:] <= np.max(M) * 0.80)
+            au = np.argwhere(M[idx:] <= np.max(M) * ratio_)
             if au.size == 0:
                 raise RuntimeError(
-                    "Peak strength does not drop 20%, please increase target ductility ratio!"
+                    f"Peak strength does not drop {1 - ratio_}, please increase target ductility ratio!"
                 )
             else:
                 bu = np.min(au) + idx - 1
@@ -213,7 +246,7 @@ def get_limit_state(
         M_u = M[bu]
         return Phi_u, M_u
 
-    def bilinearize(self, phiy: float, My: float, phiu: float, plot: bool = False):
+    def bilinearize(self, phiy: float, My: float, phiu: float, plot: bool = False, return_ax: bool = False):
         """Bilinear Approximation of Moment-Curvature Relation.
 
         Parameters
@@ -226,6 +259,8 @@ def bilinearize(self, phiy: float, My: float, phiu: float, plot: bool = False):
             The limit curvature.
         plot : bool, optional
             If True, plot the bilinear approximation, by default False.
+        return_ax: bool, default=False
+            If True, return the axes for the plot of matplotlib.
 
         Returns
         -------
@@ -241,8 +276,8 @@ def bilinearize(self, phiy: float, My: float, phiu: float, plot: bool = False):
         Phi_eq = (k * phiu - np.sqrt((k * phiu) ** 2 - 2 * k * Q)) / k
         M_eq = k * Phi_eq
 
-        M_new = np.insert(M[0 : bu + 1], 0, 0, axis=None)
-        Phi_new = np.insert(phi[0 : bu + 1], 0, 0, axis=None)
+        M_new = np.insert(M[0: bu + 1], 0, 0, axis=None)
+        Phi_new = np.insert(phi[0: bu + 1], 0, 0, axis=None)
 
         if plot:
             _, ax = plt.subplots(1, 1, figsize=(10, 10 * 0.618))
@@ -255,7 +290,7 @@ def bilinearize(self, phiy: float, My: float, phiu: float, plot: bool = False):
                 ms=12,
                 mec="black",
                 mfc="#0099e5",
-                label="Initial Yield ($\phi_y$,$M_y$)",
+                label="Initial Yield ($\\phi_y$,$M_y$)",
             )
             ax.plot(
                 Phi_eq,
@@ -264,14 +299,14 @@ def bilinearize(self, phiy: float, My: float, phiu: float, plot: bool = False):
                 ms=15,
                 mec="black",
                 mfc="#ff4c4c",
-                label="Equivalent Yield ($\phi_{{eq}}$,$M_{{eq}}$)",
+                label="Equivalent Yield ($\\phi_{{eq}}$,$M_{{eq}}$)",
             )
             maxy = np.max(ax.get_yticks())
             ax.vlines(phiu, 0, maxy, colors="#34bf49", linestyles="dashed", lw=0.75)
             txt = (
-                f"$\phi_y$={phiy:.3E}, $M_y$={My:.3E}\n"
-                f"$\phi_{{eq}}$={Phi_eq:.3E}, $M_{{eq}}$={M_eq:.3E}\n"
-                f"$\phi_{{u}}$={phiu:.3E}, $M_{{u}}$={M[bu]:.3E}"
+                f"$\\phi_y$={phiy:.3E}, $M_y$={My:.3E}\n"
+                f"$\\phi_{{eq}}$={Phi_eq:.3E}, $M_{{eq}}$={M_eq:.3E}\n"
+                f"$\\phi_{{u}}$={phiu:.3E}, $M_{{u}}$={M[bu]:.3E}"
             )
             ax.text(
                 0.5,
@@ -286,19 +321,20 @@ def bilinearize(self, phiy: float, My: float, phiu: float, plot: bool = False):
                 "Moment-Curvature",
                 fontsize=22,
             )
-            ax.set_xlabel("$\phi$", fontsize=20)
+            ax.set_xlabel("$\\phi$", fontsize=20)
             ax.set_ylabel("$M$", fontsize=20)
             plt.xticks(fontsize=15)
             plt.yticks(fontsize=15)
             # ax.set_xlim(0, np.max(ax.get_xticks()))
             ax.set_ylim(0, maxy)
-            ax.spines["bottom"].set_linewidth(1.2)
-            ax.spines["left"].set_linewidth(1.2)
-            ax.spines["right"].set_linewidth(1.2)
-            ax.spines["top"].set_linewidth(1.2)
+            for loc in ["bottom", "left", "right", "top"]:
+                ax.spines[loc].set_linewidth(1.0)
             ax.legend(loc="lower center", fontsize=15)
             plt.gcf().subplots_adjust(bottom=0.15)
-            plt.show()
+            if return_ax:
+                return ax
+            else:
+                plt.show()
         return Phi_eq, M_eq
 
 
@@ -328,7 +364,7 @@ def _create_axial_resp(p):
 
 
 def _analyze(
-    sec_tag, P=0, axis="y", max_phi=1, incr_phi=1e-5, stop_ratio=0.7, smart_analyze=True
+        sec_tag, P=0.0, axis="y", max_phi=0.5, incr_phi=1e-5, stop_ratio=0.8, smart_analyze=True, debug: bool = False
 ):
     _create_model(sec_tag=sec_tag)
     if P != 0:
@@ -351,44 +387,47 @@ def _analyze(
         userControl = {
             "analysis": "Static",
             "testType": "NormDispIncr",
-            "testTol": 1.0e-8,
+            "testTol": 1.0e-10,
             "tryAlterAlgoTypes": True,
-            "algoTypes": [10, 40, 30, 20],
+            "algoTypes": [40, 30, 20],
             "relaxation": 0.5,
             "minStep": 1.0e-12,
             "printPer": 10000000000,
-            "debugMode": False,
+            "debugMode": debug,
         }
         analysis = SmartAnalyze(analysis_type="Static", **userControl)
         segs = analysis.static_split(protocol, maxStep=incr_phi)
         ii = 0
         while True:
             seg = segs[ii]
             ii += 1
-            analysis.StaticAnalyze(2, dof, seg)
+            ok = analysis.StaticAnalyze(2, dof, seg)
             curr_M = ops.getLoadFactor(2)
             curr_Phi = ops.nodeDisp(2, dof)
             cond1 = np.abs(curr_M) < np.max(np.abs(M)) * (stop_ratio - 0.02)
-            cond2 = np.abs(curr_Phi) > max_phi
-            if cond1 or cond2:
-                break
+            cond2 = np.abs(curr_Phi) >= max_phi
             PHI.append(ops.nodeDisp(2, dof))
             M.append(curr_M)
             FIBER_RESPONSES.append(_get_fiber_sec_data(ele_tag=1))
-
+            if cond1 or cond2:
+                break
+            if ok < 0:
+                raise RuntimeError("Analysis failed!")
     else:
         ops.integrator("DisplacementControl", 2, dof, incr_phi)
         while True:
-            ops.analyze(1)
+            ok = ops.analyze(1)
             curr_M = ops.getLoadFactor(2)
             curr_Phi = ops.nodeDisp(2, dof)
             cond1 = np.abs(curr_M) < np.max(np.abs(M)) * (stop_ratio - 0.02)
             cond2 = np.abs(curr_Phi) > max_phi
-            if cond1 or cond2:
-                break
             PHI.append(ops.nodeDisp(2, dof))
             M.append(curr_M)
             FIBER_RESPONSES.append(_get_fiber_sec_data(ele_tag=1))
+            if cond1 or cond2:
+                break
+            if ok < 0:
+                raise RuntimeError("Analysis failed!")
     FIBER_RESPONSES[0] = FIBER_RESPONSES[1] * 0
     return np.abs(PHI), np.abs(M), np.array(FIBER_RESPONSES)
 
@@ -398,3 +437,9 @@ def _get_fiber_sec_data(ele_tag: int):
     # From column 1 to 6: "yCoord", "zCoord", "area", 'mat', "stress", "strain"
     fiber_data = np.array(fiber_data).reshape((-1, 6))
     return fiber_data
+
+
+def _get_center(ys, zs, areas):
+    yo = np.sum(ys * areas) / np.sum(areas)
+    zo = np.sum(zs * areas) / np.sum(areas)
+    return yo, zo