1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

# ===----------------------------------------------------------------------===//
# 
#         OpenSees - Open System for Earthquake Engineering Simulation    
#                Structural Artificial Intelligence Laboratory
#                          stairlab.berkeley.edu
# 
# ===----------------------------------------------------------------------===//
#
# Static analysis of a wrench in the plane.
#
# The mesh is created from 7 blocks:
#                 ______
#                /     /
#               /     /
#              /     /
#             /     /
#            /     /
#           /     /
#          /     /
#         /     /
#        /  7  /
#       /_____/      ____
#     /         \
#    /     6     \    70
#   /_____________\  ____
#  | \           / |
#  |  \    3    /  |  20
#  | 4 \_______/ 2 | _________
#  |___/       \___|      ____ 20
#   \5 |       | 1/   70
#    \_|       |_/        __
#      .       o     ____ __ 10
#              ^
#              #
#
# Chrystal Chern and Claudio Perez
#
import sees
import opensees.openseespy as ops
import numpy as np

angwrench = np.arctan(1/5)

# Quadrilateral blocks that comprise the wrench:
blocks = {
    1:    {
       1:  [   0,         0],
       2:  [  20,        10],
       3:  [  40, 160-70-40],
       4:  [   0, 160-70-40],
       6:  [  35,        30]},
#      6:  [  25,        30]},
    2:    {
       1:  [   0, 160-70-40],
       2:  [  40, 160-70-40],
#      6:  [  45, 160-70-20],
       3:  [  40, 160-70   ],
       4:  [ -20, 160-70-20]},
    3:    {
       1:  [ -20, 160-70-20],
       2:  [  40, 160-70   ],
       3:  [-115, 160-70   ],
       4:  [ -60,        70]},
    4:    {
       1:  [ -60,        70],
       2:  [-115,        90],
       3:  [-115,        50],
       4:  [ -75,        50]},
    5:    {
       1:  [ -75,        50],
       2:  [-115,        50],
       3:  [ -95,        10],
       4:  [ -75,         0]},
    6:    {
       1:  [  40,        90],
       5:  [  33,       112],
       2:  [   0,       160],
       3:  [-50*np.cos(angwrench), 160 + 50*np.sin(angwrench)],
       4:  [-115, 160-70   ]},
    7:    {
       1:  [   0,       160],
       2:  [250*np.sin(angwrench),                      160+250*np.cos(angwrench)],
       3:  [250*np.sin(angwrench)-50*np.cos(angwrench), 160+250*np.cos(angwrench)],
       4:  [-50*np.cos(angwrench),                      160+ 50*np.sin(angwrench)]}
}

# Subdivisions to create within each block:
divs = {
    1: (3,3),
    2: (3,3),
    3: (3,4),
    4: (3,3),
    5: (3,3),
    6: (4,4),
    7: (6,4)
}


def create_quads():
    model = ops.Model(ndm=2, ndf=2)
    model.nDMaterial("ElasticIsotropic", 1, 200e3, 0.25)

    for num,block in blocks.items():
        model.surface(divs[num],
                      element="quad", args=(1, "PlaneStrain", 1),
                      points = block)

    return model


def create_tris():
    model = ops.Model(ndm=2, ndf=2)
    model.nDMaterial("ElasticIsotropic", 1, 200e3, 0.25)


    elem = 1
    for num,block in blocks.items():
        # Because no element argument is passed, only nodes are created.
        # Next we will go back over the newly created cells and manually
        # create triangles.
        mesh = model.surface(divs[num], points = block)

        # For each new 4-node cell, create two triangles
        for cell in mesh.cells:
            nodes = mesh.cells[cell]
            model.element("tri31",   elem, (nodes[0], nodes[1], nodes[2]), 10, "PlaneStrain", 1)
            model.element("tri31", elem+1, (nodes[0], nodes[2], nodes[3]), 10, "PlaneStrain", 1)
            elem += 2

    return model


def create_boundary(model):
    # Load magnitude
    P = 700

    # Fix the first node, which is at (0.0, 0.0)
    model.fix(1, 1, 1)

    # Create a load pattern
    model.pattern("Plain", 1, "Linear")

    for node in model.getNodeTags():
        coord = model.nodeCoord(node)
        # Fix corner of block 4
        if np.linalg.norm(np.array(coord) - blocks[4][4]) < 1e-12:
            model.fix(node, 1,1)

        # Add load to the corner of block 7
        elif np.linalg.norm(np.array(coord) - blocks[7][3]) < 1e-12:
            model.load(node, (P, 0), pattern=1)



#model = create_quads()
model = create_tris()
create_boundary(model)

model.analysis("Static")
model.integrator("LoadControl", 1)
model.analyze(1)


# Render the deformed shape
sees.serve(sees.render(model, lambda i: [500*u for u in model.nodeDisp(i)], canvas="gltf"))