Plot stress distribution of plane stress quad model
Example .py file can be downloaded here:
1# import openseespy.opensees as ops
2import opensees as ops # local compilation
3import opsvis as opsv
4
5import matplotlib.pyplot as plt
6
7ops.wipe()
8ops.model('basic', '-ndm', 2, '-ndf', 2)
9ops.node(1, 0., 0.)
10ops.node(2, 0., 1.)
11ops.node(3, 0., 2.)
12ops.node(4, 0., 3.)
13ops.node(5, 0., 4.)
14ops.node(6, 1., 0.)
15ops.node(7, 1., 1.)
16ops.node(8, 1., 2.)
17ops.node(9, 1., 3.)
18ops.node(10, 1., 4.)
19ops.node(11, 2., 0.)
20ops.node(12, 2., 1.)
21ops.node(13, 2., 2.)
22ops.node(14, 2., 3.)
23ops.node(15, 2., 4.)
24ops.node(16, 3., 0.)
25ops.node(17, 3., 1.)
26ops.node(18, 3., 2.)
27ops.node(19, 3., 3.)
28ops.node(20, 3., 4.)
29ops.node(21, 4., 0.)
30ops.node(22, 4., 1.)
31ops.node(23, 4., 2.)
32ops.node(24, 4., 3.)
33ops.node(25, 4., 4.)
34
35ops.nDMaterial('ElasticIsotropic', 1, 1000, 0.3)
36
37ops.element('quad', 1, 1, 6, 7, 2, 1, 'PlaneStress', 1)
38ops.element('quad', 2, 2, 7, 8, 3, 1, 'PlaneStress', 1)
39ops.element('quad', 3, 3, 8, 9, 4, 1, 'PlaneStress', 1)
40ops.element('quad', 4, 4, 9, 10, 5, 1, 'PlaneStress', 1)
41ops.element('quad', 5, 6, 11, 12, 7, 1, 'PlaneStress', 1)
42ops.element('quad', 6, 7, 12, 13, 8, 1, 'PlaneStress', 1)
43ops.element('quad', 7, 8, 13, 14, 9, 1, 'PlaneStress', 1)
44ops.element('quad', 8, 9, 14, 15, 10, 1, 'PlaneStress', 1)
45ops.element('quad', 9, 11, 16, 17, 12, 1, 'PlaneStress', 1)
46ops.element('quad', 10, 12, 17, 18, 13, 1, 'PlaneStress', 1)
47ops.element('quad', 11, 13, 18, 19, 14, 1, 'PlaneStress', 1)
48ops.element('quad', 12, 14, 19, 20, 15, 1, 'PlaneStress', 1)
49ops.element('quad', 13, 16, 21, 22, 17, 1, 'PlaneStress', 1)
50ops.element('quad', 14, 17, 22, 23, 18, 1, 'PlaneStress', 1)
51ops.element('quad', 15, 18, 23, 24, 19, 1, 'PlaneStress', 1)
52ops.element('quad', 16, 19, 24, 25, 20, 1, 'PlaneStress', 1)
53
54ops.fix(1, 1, 1)
55ops.fix(6, 1, 1)
56ops.fix(11, 1, 1)
57ops.fix(16, 1, 1)
58ops.fix(21, 1, 1)
59
60ops.equalDOF(2, 22, 1, 2)
61ops.equalDOF(3, 23, 1, 2)
62ops.equalDOF(4, 24, 1, 2)
63ops.equalDOF(5, 25, 1, 2)
64
65ops.timeSeries('Linear', 1)
66ops.pattern('Plain', 1, 1)
67ops.load(15, 0., -1.)
68
69ops.analysis('Static')
70ops.analyze(1)
71
72# - plot model
73opsv.plot_model()
74plt.axis('equal')
75
76# - plot deformation
77plt.figure()
78opsv.plot_defo()
79plt.axis('equal')
80
81# get values at OpenSees nodes
82sig_out = opsv.sig_out_per_node()
83print(f'sig_out:\n{sig_out}')
84
85# !!! select from sig_out: e.g. vmises
86# j, jstr = 0, 'sxx'
87j, jstr = 1, 'syy'
88# j, jstr = 2, 'sxy'
89# j, jstr = 3, 'vmis'
90# j, jstr = 4, 's1'
91# j, jstr = 5, 's2'
92# j, jstr = 6, 'alfa'
93
94nds_val = sig_out[:, j]
95
96plt.figure()
97opsv.plot_stress_2d(nds_val)
98
99plt.xlabel('x [m]')
100plt.ylabel('y [m]')
101plt.title(f'{jstr}')
102
103# for educational purposes show values at integration points and
104# nodes which can finally be averaged at nodes
105eles_ips_crd, eles_nds_crd, nds_crd, quads_conn = opsv.quad_crds_node_to_ip()
106
107print(f'\neles_ips_crd:\n{eles_ips_crd}')
108print(f'\neles_nds_crd:\n{eles_nds_crd}')
109print(f'\nnds_crd:\n{nds_crd}')
110print(f'\nquads_conn:\n{quads_conn}')
111
112eles_ips_sig_out, eles_nds_sig_out = opsv.sig_out_per_ele_quad()
113
114print(f'\neles_ips_sig_out:\n{eles_ips_sig_out}')
115print(f'\neles_nds_sig_out:\n{eles_nds_sig_out}')
116
117sig_out_indx = j # same as j, jstr
118
119opsv.plot_mesh_with_ips_2d(nds_crd, eles_ips_crd, eles_nds_crd, quads_conn,
120 eles_ips_sig_out, eles_nds_sig_out, sig_out_indx)
121
122plt.xlabel('x [m]')
123plt.ylabel('y [m]')
124
125plt.show()
126
127exit()