clean up

Author: claudioperez

content/en/benchmarks/frame-0002/cantilever.py

@@ -3,7 +3,7 @@
 import xara.units.iks as units
 import veux
 import os
-from xsection.library import from_aisc
+from xsection.library import from_aisc, Rectangle
 import numpy as np
 
 def create_cantilever(shape, length, element, material,
@@ -70,16 +70,19 @@ def analyze(model, P):
     model.pattern("Plain", 1, "Linear", load={nn-1: [0,0,1, 0,0,0]})
     model.integrator("LoadControl", P)
     model.analysis("Static")#, pattern=1, integrator=1)
+    model.test("Residual", 1e-10, 3)
 
-    model.analyze(1)
+    assert model.analyze(1) == 0
 
 
 
 
 if __name__ == "__main__":
-    shape = from_aisc("W14x48", units=units)
+#   shape = Rectangle(d=14, b=10, mesh_scale=1/100) #
+    shape = from_aisc("W14x48", units=units, mesh_scale=1/100)
+
     E  = 29e3*units.ksi
-    G  = 11.2e3*units.ksi * 4/5
+    G  = 11.2e3*units.ksi # * 4/5
 
     A  = shape.cnn()[0,0]
     GA = A*G
@@ -90,10 +93,10 @@ def analyze(model, P):
         model = create_cantilever(shape,
                                 length=L,
                                 material={"name": "ElasticIsotropic", "E": E, "G": G},
-                                element=os.environ.get("Element", "ForceFrame"), 
-                                section="ShearFiber", 
+                                element=os.environ.get("Element", "ForceFrame"),
+                                section="ShearFiber",
                                 shear=shear,
-                                ne=6, 
+                                ne=1,
                                 nen=2
                 )
 
@@ -105,9 +108,9 @@ def analyze(model, P):
         u_euler = P*L**3/(3*EI)
         u_shear = P*L/GA * shear
 
-        print(f"Uz = {uz:.3f}, Uz theory = {u_euler+u_shear:.3f} ({u_euler:.3f} + {u_shear:.6f})")
+        print(f"Uz = {uz:.6f}, Uz theory = {u_euler+u_shear:.6f} ({u_euler:.6f} + {u_shear:.6f})")
 
-        model.eval(f"verify value [nodeDisp {end} 3] {u_euler+u_shear:.6f} 1e-4")
+        model.eval(f"verify value [nodeDisp {end} 3] {u_euler+u_shear:.12f} 1e-6")
 
     # a = veux.create_artist(model)
     # a.draw_sections()

content/en/benchmarks/frame-0011/e0011.py

@@ -1,11 +1,13 @@
 # Test of the warping DOF
 # Linear 7-DOF analysis of a cantilever subjected to a torque
+#
 import os
 import sys
 
 import veux
 from veux.motion import Motion
 from xsection.library import WideFlange, HollowRectangle, Channel, Rectangle, Circle
+
 import opensees.openseespy as ops
 
 # External libraries
@@ -31,7 +33,7 @@ def create_cantilever(aspect,
     nen = 2
     nn = ne*(nen-1)+1
 
-    model = ops.Model(ndm=3, ndf=7) #, echo_file=open(f"test-{case}-{element[:5]}.tcl", "w+"))
+    model = ops.Model(ndm=3, ndf=7)
 
     model.eval(f"set E {E}")
     model.eval(f"set G {G}")
@@ -43,11 +45,11 @@ def create_cantilever(aspect,
 
 
     if section == "Elastic":
-        cmm = shape.torsion.cmm()
-        cnn = shape.torsion.cnn()
-        cnv = shape.torsion.cnv()
-        cnm = shape.torsion.cnm()
-        cmw = shape.torsion.cmw()
+        cmm = shape.cmm()
+        cnn = shape.cnn()
+        cnv = shape.cnv()
+        cnm = shape.cnm()
+        cmw = shape.cmw()
         A = cnn[0,0]
         model.section("ElasticFrame", sec,
                         E=E,
@@ -66,20 +68,27 @@ def create_cantilever(aspect,
                         Sz=-cmw[2,0]
         )
     else:
+
+
         model.section("ShearFiber", 1, GJ=0)
+
+
         for fiber in shape.create_fibers():
+
             model.fiber(**fiber, material=mat, section=1)
 
 
+
     model.geomTransf("Linear", 1, (0,0,1))
 
+
     for i,x in enumerate(np.linspace(0, L, nn)):
         model.node(i, (x,0,0))
 
     for i in range(ne):
         start = i * (nen - 1)
         nodes = list(range(start, start + nen))
-        model.element(element, i+1, nodes, section=sec, transform=1, shear=1)
+        model.element(element, i+1, nodes, section=1, transform=1, shear=1)
 
     wi = int(case in "cb")
     wj = int(case in "c")
@@ -124,6 +133,9 @@ def create_cantilever(aspect,
                 mesh_scale=1/200
             )
 
+
+
+
     print(shape.summary())
 
     fig,  ax  = plt.subplots()
@@ -183,7 +195,7 @@ def create_cantilever(aspect,
             P.append(model.getTime())
             if model.analyze(1) != 0:
                 print(f"Failed at time = {model.getTime()}")
-                break 
+                break
 
         model.reactions()
 
@@ -206,14 +218,14 @@ def create_cantilever(aspect,
 
         xs = [(x+model.nodeCoord(model.eleNodes(tag)[0], 1))/L for tag in model.getEleTags() for x in model.eleResponse(tag, "integrationPoints")]
         ws = [
-            model.eleResponse(tag, "section", i+1, "resultant")[6]/(Mmax*L) 
-            for tag in model.getEleTags() 
+            model.eleResponse(tag, "section", i+1, "resultant")[6]/(Mmax*L)
+            for tag in model.getEleTags()
             for i in range(len(model.eleResponse(tag, "integrationPoints")))
         ]
         ax4.plot(xs, ws, "-", color=color, label=case)
         vs = [
             model.eleResponse(tag, "section", i+1, "resultant")[9]/Mmax
-            for tag in model.getEleTags() 
+            for tag in model.getEleTags()
             for i in range(len(model.eleResponse(tag, "integrationPoints")))
         ]
         ax5.plot(xs, vs, "-", color=color)
@@ -245,5 +257,3 @@ def create_cantilever(aspect,
     ax4.figure.savefig("img/e0011-resultants.png")
     plt.show()
 
-
-

content/en/benchmarks/frame-0011/e0011_6dof.py

@@ -39,7 +39,7 @@ def create_cantilever(case="a", element="ExactFrame", section="Elastic"):
                 d  = 21.62,
                 b  = 8.42,
                 # saint_venant="linear"
-            ).create_shape()
+            )
 
     if section == "Elastic":
         cmm = shape.torsion.cmm()
@@ -66,9 +66,8 @@ def create_cantilever(case="a", element="ExactFrame", section="Elastic"):
         )
     else:
         model.section("ShearFiber", 1, GJ=0)
-        for fiber in shape.fibers():
-            y, z = fiber.location
-            model.fiber(y, z, fiber.area, warp=[fiber.warp[0], [0,0,0], [0,0,0]], material=mat, section=1)
+        for fiber in shape.create_fibers():
+            model.fiber(**fiber, material=mat, section=1)
 
 
     model.geomTransf("Linear", 1, (0,0,1))

content/en/benchmarks/frame-0013/a.py

@@ -0,0 +1,19 @@
+import opensees.openseespy as ops
+
+m = ops.Model(ndm=3, ndf=6)
+
+m.node(1, (0,0,0))
+m.node(2, (1,2,3))
+m.fix(1, (1,1,1, 1,1,1))
+
+m.geomTransf("Corotational", 1, (0,0,1))
+m.section("ElasticFrame", 1, A=1, Iy=1, Iz=1, J=1, G=1, E=1, Ay=1, Az=1)
+m.element("ForceFrame", 1, (1,2), transform=1, section=1)
+
+m.pattern("Plain", 1, "Linear", loads={2: (0.2,1,1, 0, 0.3,0)})
+m.integrator("LoadControl", 0.001)
+m.test("FixedNumIter", 1)
+m.analysis("Static")
+m.analyze(1, operation="increment")
+m.analyze(1, operation="iteration")
+m.analyze(1, operation="iteration")

content/en/benchmarks/frame-0013/e0013.py

@@ -0,0 +1,177 @@
+# Cantilever with channel section and eccentric loading
+#
+import veux
+from veux.motion import Motion
+from shps.shapes import Channel
+import opensees.openseespy as ops
+
+# External libraries
+import numpy as np
+import matplotlib.pyplot as plt
+try:
+    plt.style.use("veux-web")
+except:
+    pass
+
+def create_cantilever(ne, offset, element, section, nen=2):
+    model = ops.Model(ndm=3, ndf=6)
+
+    E = 2.1e4 # MPa, or 210 GPa
+    v = 0.30 #0.5*E/G - 1
+    G = 0.5*E/(1+v) # 8076.92
+
+    nmn = ne*(nen-1)+1
+    L  = 900
+
+
+    mat = 1
+    sec = 1
+    model.material('ElasticIsotropic', mat, E, v) #G=G)
+
+    shape = Channel(d=30, b=10, tf=1.6, tw=1.0).create_shape()
+
+    shape = shape.translate(offset)
+
+
+    # _m = shape.mesh
+    # veux.serve(veux.render((_m.nodes, _m.cells())))
+
+    if "fiber" in section.lower():
+        warp = os.environ.get("Warp", None)
+        print(f"Section = Fiber; ", warp)
+        model.section("ShearFiber", sec, GJ=0)
+
+        for fiber in shape.fibers(warp=warp):
+            model.fiber(y, z, fiber.area, mat, fiber.warp[0], section=sec) #fiber.warp[1], section=sec)
+
+    else:
+        print("Section = Elastic")
+        cmm = shape.torsion.cmm()
+        cnn = shape.torsion.cnn()
+        cnv = shape.torsion.cnv()
+        cnm = shape.torsion.cnm()
+        cmw = shape.torsion.cmw()
+        A = cnn[0,0]
+        model.section("ElasticFrame", sec,
+                        E=E,
+                        G=G,
+                        A=A,
+                        Ay=1*A,
+                        Az=1*A,
+                        Qy=cnm[0,1],
+                        Qz=cnm[2,0],
+                        Iy=cmm[1,1],
+                        Iz=cmm[2,2],
+                        J =shape.torsion.torsion_constant(),
+                        Ry= cnv[1,0],
+                        Rz=-cnv[2,0],
+                        Sy= cmw[1,0],
+                        Sz=-cmw[2,0]
+        )
+
+
+#   model.geomTransf("Corotational", 1, (0,1,0))
+    model.geomTransf("Corotational", 1, (0,0,1))
+
+    for i,x in enumerate(np.linspace(0, L, nmn)):
+        model.node(i, (x,0,0))
+
+    model.fix(0,  (1,1,1,  1,1,1))
+
+    for i in range(ne):
+        start = i * (nen - 1)
+        nodes = list(range(start, start + nen))
+        model.element(element, i+1, nodes, section=sec, transform=1)
+
+
+    return model, shape
+
+def analyze(element, section, pattern="node", nen=2):
+    ne = 20
+    en = ne*(nen-1)
+    if pattern == "node":
+        offset = (0, 15)
+    else:
+        offset = (0,  0)
+
+    model,shape = create_cantilever(ne, offset, element=element, section=section, nen=nen)
+    if False:
+        artist = veux.create_artist(model, model_config=dict(extrude_outline=shape))
+        artist.draw_nodes(size=10)
+        # artist.draw_sections()
+        # veux.serve(artist)
+        motion = Motion(artist)
+
+    #
+    # Apply vertical load
+    #
+    speed  = 1 # animation frames
+    Pmax   = 20 # kN
+    model.pattern("Plain", 1, "Linear")
+    if pattern == "node":
+        print("Pattern = node")
+        model.load(en, (0,0,-1,  0,0,0), pattern=1)
+#       model.load(ne, (0,-1,0,  0,0,0), pattern=1)
+
+    else:
+        print("Pattern = element")
+        model.eleLoad("Frame", "Dirac",
+                      force = [0, 0, -1],
+                      basis = "global",
+                      offset=[1.0,0, -15],
+                      pattern=1,
+                      elements=[ne]
+        )
+
+    model.system('Umfpack')
+    model.integrator("LoadControl", Pmax/1000)#, 8, Pmax/500, Pmax/2)
+    model.test("NormDispIncr", 1e-8, 100, 0)
+#   model.test('NormUnbalance',1e-6,10,1)
+    model.algorithm("Newton")
+    model.analysis("Static")
+#   input()
+
+    u = []
+    v = []
+    w = []
+    P = []
+#   for i in range(50):
+    while model.getTime() < Pmax:
+#       motion.advance(time=model.getTime()*speed)
+#       motion.draw_sections(rotation=model.nodeRotation,
+#                            position=model.nodeDisp)
+#       u.append(-model.nodeDisp(ne, 1))
+#       v.append( model.nodeDisp(ne, 2))
+#       w.append(-model.nodeDisp(ne, 3))
+#       P.append( model.getTime())
+        if model.analyze(1) != 0:
+            print(f"Failed at time = {model.getTime()} with v = {v[-1]}")
+            break
+
+
+    if False:
+        fig, ax = plt.subplots()
+        ax.set_xlabel(r"Displ, $v$")
+        ax.set_ylabel("Load, $P$")
+        # ax.set_xlim([0,    300])
+        ax.set_ylim([0,   Pmax])
+        ax.axvline(0, color='black', linestyle='-', linewidth=1)
+        ax.axhline(0, color='black', linestyle='-', linewidth=1)
+        ax.plot(u, P, label="$u$")
+        ax.plot(v, P, label="$v$")
+        ax.plot(w, P, label="$w$")
+        ax.legend()
+        fig.savefig("img/e0013.png")
+        plt.show()
+
+        motion.add_to(artist.canvas)
+        veux.serve(artist)
+
+if __name__ == "__main__":
+    import os
+    analyze(pattern = os.environ.get("Pattern", "node"),
+            element = os.environ.get("Element", "ExactFrame"),
+            section = os.environ.get("Section", "ShearFiber"),
+            nen=2
+            )
+

content/en/benchmarks/frame-0013/index.md

@@ -25,6 +25,7 @@ $$
 
 Beams:
 
+- Rong et al 2020, Section 7.2
 - Hsiao KM, Lin W Y (2000) A co-rotational formulation for thin-walled beams with monosymmetric open section. DOI: 10.1016/S0045-7825(99)00471-5
 - D Manta, R Goncalves (2016) A geometrically exact Kirchhoff beam model including torsion warping.