cmp - update cyclic example

Author: claudioperez

content/en/examples/FrameCycling/column_cycles.py

@@ -5,6 +5,7 @@
 # Written: Vesna Terzic (vesna@berkeley.edu)
 # Created: 12/2011
 #
+import os
 from math import sqrt, pi
 import opensees.openseespy as ops
 import matplotlib.pyplot as plt
@@ -170,19 +171,9 @@ def lehman_section(model, D, clearCover):
     # Aggregate shear to the RC section
     secTag = secnTag + 1
     model.section("Aggregator", secTag, steelTag + 6, "Vy", "-section", secnTag)
-    return secTag 
+    return secTag
 
-if __name__ == "__main__":
-    # Set element type (force-based = 1, displacement-based = 2)
-    eleType = 1
-
-    # Number of finite elements and integration points per element
-    if eleType == 1:
-        ne = 1
-        nIP = 5
-    elif eleType == 2:
-        ne = 4
-        nIP = 3
+def create_column(ne, nIP, element):
 
     # Create a 2D model with 3 DOFs per node
     model = ops.Model(ndm=2, ndf=3)
@@ -211,12 +202,15 @@ def lehman_section(model, D, clearCover):
     model.geomTransf("Corotational", transfTag)
 
     for i in range(ne):
-        if eleType == 1:
-            model.element("forceBeamColumn", i + 1, (i + 1, i + 2), nIP, secTag, transfTag)
-        elif eleType == 2:
-            model.element("dispBeamColumn", i + 1, (i + 1, i + 2), nIP, secTag, transfTag)
+        model.element(element, i + 1, (i + 1, i + 2), nIP, secTag, transfTag)
+
+    return model
 
 
+def analyze_column(ne, nIP, eleType):
+
+    model = create_column(ne, nIP, eleType)
+
     # Define Gravity Load
     IDctrlNode = ne + 1
 
@@ -268,8 +262,8 @@ def lehman_section(model, D, clearCover):
     model.sp(IDctrlNode, IDctrlDOF, 1.0, pattern=2)
 
     # Define recorders for displacement and force
-    model.recorder("Node", "disp", "-file", "out/Disp.out", "-time", "-node", ne + 1, dof=1)
-    model.recorder("Node", "reaction", "-file", "out/Force.out", "-time", "-node", 1, dof=1)
+    model.recorder("Node", "disp", "-time", file="out/_disp.out", node=ne+1, dof=1)
+    model.recorder("Node", "reaction", "-time", file="out/_force.out", node=1, dof=1)
 
     # Cyclic analysis
     model.constraints("Penalty", 1.0e14, 1.0e14)
@@ -281,5 +275,52 @@ def lehman_section(model, D, clearCover):
     model.analysis("Static")
 
     ok = model.analyze(anpts)
+    if ok != 0:
+        raise Exception("Analysis failed")
+
+    model.wipe()
 
+    with open("out/_disp.out", "r") as f:
+        lines = f.readlines()
+        disp = [float(line.split()[1]) for line in lines]
+        time = [float(line.split()[0]) for line in lines]
+    with open("out/_force.out", "r") as f:
+        lines = f.readlines()
+        force = [float(line.split()[1]) for line in lines]
 
+    os.remove("out/_disp.out")
+    os.remove("out/_force.out")
+
+
+    return time, disp, force
+
+if __name__ == "__main__":
+
+    fig, ax = plt.subplots()
+
+    # Set element type (force-based = 1, displacement-based = 2)
+    for element in "dispBeamColumn", "forceBeamColumn":
+        # Number of finite elements and integration points per element
+        if "force" in element:
+            ne = 1
+            nIP = 5
+        else:
+            # Displacement formulation
+            ne = 4
+            nIP = 3
+
+        time, disp, force = analyze_column(ne, nIP, element)
+
+        ax.plot(disp, force, label=element)
+
+    ax.set_xlabel("Displacement [in]")
+    ax.set_ylabel("Force [kips]")
+    fig.legend()
+    plt.show()
+
+
+    fig, ax = plt.subplots()
+    ax.plot(time, disp)
+    ax.set_xlabel("Time [s]")
+    ax.set_ylabel("Displacement [in]")
+    plt.show()

content/en/examples/FrameCycling/index.md

@@ -37,18 +37,15 @@ The intent of this example is to show how to properly model inelastic beam-colum
 ## Model Generation
 
 - The column is modeled using either force-based or displacement-based elements.
-- The column height <code>H</code> is 96 inches, and the diameter <code>D</code> is 24 inches.
-- The clear cover of concrete <code>clearCover</code> is 0.75 inches.
+- The column height <code>H</code> is 96 inches
 
 ```python
 import xara
+
 # Create a 2D model with 3 DOFs per node
 model = xara.Model(ndm=2, ndf=3)
 
-# Input parameters
 H = 96.0 * inch  # Column height
-D = 24.0 * inch  # Column diameter
-clearCover = 0.75 * inch  # Clear cover of concrete
 
 # Define nodes
 model.node(1, (0.0, 0.0))
@@ -108,7 +105,13 @@ $$
 
 ### Section
 
+- The diameter <code>D</code> is 24 inches.
+- The clear cover of concrete <code>clearCover</code> is 0.75 inches.
+
 ```python
+D = 24.0 * inch  # Column diameter
+clearCover = 0.75 * inch  # Clear cover of concrete
+
 theta = 360.0 / numBars
 model.section("fiberSec", secnTag, GJ=1e8)
 # Core patch
@@ -141,6 +144,12 @@ for i in range(ne):
 
 ### Cycling
 
+```python
+# Define recorders for displacement and force
+model.recorder("Node", "disp", "-time", file="out/Disp.out", node=ne+1, dof=1)
+model.recorder("Node", "reaction", "-time", file="out/Force.out", node=1, dof=1)
+```
+
 ## Resources
 
 The slides from the original presentation can be downloaded <a href="https://opensees.berkeley.edu/wiki/images/c/c5/FFvsDF.pdf">here</a>

content/en/examples/FrameCycling/plot.py

@@ -0,0 +1,33 @@
+import matplotlib.pyplot as plt
+plt.style.use("veux-web")
+import numpy as np
+
+def main(argv):
+    x = None 
+    y_files = []
+    argi = iter(argv)
+    for arg in argi:
+        if arg == "--x":
+            x = next(argi)
+        else:
+            if x is None:
+                x = np.loadtxt(arg)
+            else:
+                y_files.append(arg)
+
+    print(y_files)
+    ys = [
+        np.loadtxt(f)[:,-1]
+        for f in y_files
+    ]
+    if len(ys) > 0:
+        for y in ys:
+            plt.plot(x, y)
+    else:
+        plt.plot(x)
+
+    plt.show()
+
+if __name__ == "__main__":
+    import sys
+    main(sys.argv[1:])

content/en/examples/FrameOffsets/diamond.py

@@ -14,7 +14,7 @@ def diamond_solution(EI, EA, L, off=0):
                    [0, 0,              1],
                    [0, 0, off/np.sqrt(2)]])
 
-    A  = Ag@Ao #Ag[:,:3]@Ao[:3,:]
+    A  = Ag@Ao
     Ke = np.array([[EA/Le,       0,       0],
                    [    0, 4*EI/Le, 2*EI/Le],
                    [    0, 2*EI/Le, 4*EI/Le]])

content/en/examples/FrameOffsets/index.md

@@ -1,6 +1,35 @@
 ---
 title: Rigid Offsets
-draft: true
+draft: false
 ---
 
+In this example we investigate the `offset` argument for frame geometric transformations.
+
+$$
+\boldsymbol{A}_v = \begin{bmatrix}
+-dX/L_e   &  dY/L_e   & 0 &  dX/L_e   \\
+-dY/L_e^2 & -dX/L_e^2 & 0 &  dY/L_e^2 \\
+-dY/L_e^2 & -dX/L_e^2 & 1 &  dY/L_e^2 \\
+\end{bmatrix}
+$$
+
+$$
+\boldsymbol{A}_o = 
+\begin{bmatrix}
+0 & 1 &               0 \\
+1 & 0 & -L_o/\sqrt{2}   \\
+0 & 0 &               1 \\
+0 & 0 &  L_o/\sqrt{2}   \\
+\end{bmatrix}
+$$
+
+$$
+\boldsymbol{K}_e = \begin{bmatrix}
+EA/L_e &       0 &       0 \\
+    0 & 4 EI/L_e & 2 EI/L_e \\
+    0 & 2 EI/L_e & 4 EI/L_e \\
+\end{bmatrix}
+$$
+
+{{< fold diamond.py >}}