import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from matplotlib.collections import PolyCollection
from matplotlib.patches import Circle, Polygon, Wedge
from matplotlib.animation import FuncAnimation
from matplotlib.path import Path
import matplotlib.tri as tri
from math import isclose
from .settings import *
from .defo import *
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d.proj3d import proj_transform
# from https://gist.github.com/WetHat/1d6cd0f7309535311a539b42cccca89c
class Arrow3D(FancyArrowPatch):
def __init__(self, x, y, z, dx, dy, dz, *args, **kwargs):
super().__init__((0, 0), (0, 0), *args, **kwargs)
self._xyz = (x, y, z)
self._dxdydz = (dx, dy, dz)
def draw(self, renderer):
x1, y1, z1 = self._xyz
dx, dy, dz = self._dxdydz
x2, y2, z2 = (x1 + dx, y1 + dy, z1 + dz)
xs, ys, zs = proj_transform((x1, x2), (y1, y2), (z1, z2), self.axes.M)
self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
super().draw(renderer)
def do_3d_projection(self, renderer=None):
x1, y1, z1 = self._xyz
dx, dy, dz = self._dxdydz
x2, y2, z2 = (x1 + dx, y1 + dy, z1 + dz)
xs, ys, zs = proj_transform((x1, x2), (y1, y2), (z1, z2), self.axes.M)
self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
return np.min(zs)
def _arrow3D(ax, x, y, z, dx, dy, dz, *args, **kwargs):
'''Add an 3d arrow to an `Axes3D` instance.'''
arrow = Arrow3D(x, y, z, dx, dy, dz, *args, **kwargs)
ax.add_artist(arrow)
setattr(Axes3D, 'arrow3D', _arrow3D)
def _plot_model_2d(node_labels, element_labels, offset_nd_label, axis_off,
fig_wi_he, fig_lbrt, nodes_only, fmt_model,
fmt_model_nodes_only, node_supports, gauss_points, fmt_gauss_points,
fmt_model_truss, truss_node_offset, ax):
if not ax:
if fig_wi_he:
fig_wi, fig_he = fig_wi_he
fig, ax = plt.subplots(figsize=(fig_wi / 2.54, fig_he / 2.54))
else:
fig, ax = plt.subplots()
if fig_lbrt:
fleft, fbottom, fright, ftop = fig_lbrt
fig.subplots_adjust(left=fleft, bottom=fbottom, right=fright, top=ftop)
max_x_crd, max_y_crd, max_crd = -np.inf, -np.inf, -np.inf
node_tags = ops.getNodeTags()
ele_tags = ops.getEleTags()
# nodes only in the OpenSees domain
if len(ele_tags) == 0 or nodes_only:
for node_tag in node_tags:
ax.plot(ops.nodeCoord(node_tag)[0],
ops.nodeCoord(node_tag)[1],
**fmt_model_nodes_only)
if node_labels:
ax.text(ops.nodeCoord(node_tag)[0],
ops.nodeCoord(node_tag)[1],
f'{node_tag}', va='bottom', ha='left', color='blue')
ax.axis('equal')
ax.grid(False)
return ax
for node_tag in node_tags:
x_crd = ops.nodeCoord(node_tag)[0]
y_crd = ops.nodeCoord(node_tag)[1]
if x_crd > max_x_crd:
max_x_crd = x_crd
if y_crd > max_y_crd:
max_y_crd = y_crd
max_crd = np.amax([max_x_crd, max_y_crd])
_offset = 0.005 * max_crd
_offnl = 0.003 * max_crd
# 0. first print node labels
if node_labels:
for j, node_tag in enumerate(node_tags):
# xycrd = ops.nodeCoord(node_tag)
if not offset_nd_label == 'above':
offset_nd_label_x, offset_nd_label_y = _offnl, _offnl
va = 'bottom'
# va = 'center'
ha = 'left'
else:
offset_nd_label_x, offset_nd_label_y = 0.0, _offnl
va = 'bottom'
ha = 'center'
ax.text(ops.nodeCoord(node_tag)[0]+offset_nd_label_x,
ops.nodeCoord(node_tag)[1]+offset_nd_label_y,
f'{node_tag}', va=va, ha=ha, color='blue')
for i, ele_tag in enumerate(ele_tags):
ele_classtag = ops.getEleClassTags(ele_tag)[0]
if (ele_classtag == EleClassTag.ElasticBeam2d
or ele_classtag == EleClassTag.ForceBeamColumn2d
or ele_classtag == EleClassTag.DispBeamColumn2d
or ele_classtag == EleClassTag.TimoshenkoBeamColumn2d
or ele_classtag == EleClassTag.ElasticTimoshenkoBeam2d
or ele_classtag == EleClassTag.MVLEM
or ele_classtag == EleClassTag.SFI_MVLEM
or ele_classtag == EleClassTag.DispBeamColumn2dThermal
or ele_classtag == EleClassTag.ForceBeamColumn2dThermal):
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd_nodes = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd_nodes[i, :] = ops.nodeCoord(ele_node_tag)
ecrd_eles = np.copy(ecrd_nodes)
ele_offsets = np.array(ops.eleResponse(ele_tag, 'offsets'))
nz_offsets = np.nonzero(ele_offsets)[0] # tuple of arrays
if np.any(nz_offsets):
# modify ecrd_eles
ecrd_eles[:, 0] += ele_offsets[[0, 3]]
ecrd_eles[:, 1] += ele_offsets[[1, 4]]
ax.plot([ecrd_nodes[0, 0], ecrd_eles[0, 0]],
[ecrd_nodes[0, 1], ecrd_eles[0, 1]],
**fmt_model_rigid_offset)
ax.plot([ecrd_nodes[1, 0], ecrd_eles[1, 0]],
[ecrd_nodes[1, 1], ecrd_eles[1, 1]],
**fmt_model_rigid_offset)
else:
pass
ax.plot(ecrd_eles[:, 0], ecrd_eles[:, 1], **fmt_model)
if gauss_points:
Lgps = ops.eleResponse(ele_tag, 'integrationPoints')
for Lgpi in Lgps:
dxi = (ecrd_eles[1, 0] - ecrd_eles[0, 0]) * Lgpi / bar_length(ecrd_eles)
dyi = (ecrd_eles[1, 1] - ecrd_eles[0, 1]) * Lgpi / bar_length(ecrd_eles)
ax.plot(ecrd_eles[0, 0] + dxi, ecrd_eles[0, 1] + dyi, **fmt_gauss_points)
# location of label
xt = sum(ecrd_eles[:, 0]) / nen
yt = sum(ecrd_eles[:, 1]) / nen
if element_labels:
if ecrd_eles[1, 0] - ecrd_eles[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y = _offset, 0.0
elif ecrd_eles[1, 1] - ecrd_eles[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y = 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y = 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, f'{ele_tag}', va=va, ha=ha,
color='red')
if ele_classtag in [EleClassTag.truss,EleClassTag.trussSection,EleClassTag.CorotTruss]:
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
if truss_node_offset:
Lx = ecrd[1, 0] - ecrd[0, 0]
Ly = ecrd[1, 1] - ecrd[0, 1]
dLx = truss_node_offset * Lx
dLy = truss_node_offset * Ly
ecrd[0, 0] += dLx
ecrd[1, 0] -= dLx
ecrd[0, 1] += dLy
ecrd[1, 1] -= dLy
ax.plot(ecrd[:, 0], ecrd[:, 1], **fmt_model_truss)
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y = _offset, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y = 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y = 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, f'{ele_tag}', va=va, ha=ha,
color='red')
# 2d triangular (tri31) elements plot_model
elif (ele_classtag == EleClassTag.tri3n):
nen = 3
nodes_geo_order = [0, 1, 2, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1], **fmt_model)
if element_labels:
va = 'center'
ha = 'center'
ax.text(xt, yt, f'{ele_tag}', va=va, ha=ha, color='red')
# 2d quadrilateral (quad) elements plot_model
elif (ele_classtag == EleClassTag.quad4n):
nen = 4
nodes_geo_order = [0, 1, 2, 3, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1], **fmt_model)
if element_labels:
va = 'center'
ha = 'center'
ax.text(xt, yt, f'{ele_tag}', va=va, ha=ha, color='red')
# 2d quadrilateral (quad8n) elements plot_model
elif (ele_classtag == EleClassTag.quad8n):
nen = 8
nodes_geo_order = [0, 4, 1, 5, 2, 6, 3, 7, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1], **fmt_model)
if element_labels:
va = 'center'
ha = 'center'
ax.text(xt, yt, f'{ele_tag}', va=va, ha=ha, color='red')
# 2d quadrilateral (quad9n) elements plot_model
elif (ele_classtag == EleClassTag.quad9n):
nen = 9
nodes_geo_order = [0, 4, 1, 5, 2, 6, 3, 7, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1], **fmt_model)
if element_labels:
va = 'center'
ha = 'center'
ax.text(xt, yt, f'{ele_tag}', va=va, ha=ha, color='red')
ax.scatter([ecrd[8, 0]], [ecrd[8, 1]], s=2, color='g')
# 2d triangle (tri6n) elements plot_model
elif (ele_classtag == EleClassTag.tri6n):
nen = 6
nodes_geo_order = [0, 3, 1, 4, 2, 5, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1], **fmt_model)
if element_labels:
va = 'center'
ha = 'center'
ax.text(xt, yt, f'{ele_tag}', va=va, ha=ha, color='red')
# 2d joint2d element plot_model
elif (ele_classtag == EleClassTag.Joint2D):
nen = 5
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label - central node
xt = ecrd[4, 0]
yt = ecrd[4, 1]
ax.plot([ecrd[0, 0], ecrd[2, 0], ecrd[2, 0], ecrd[0, 0],
ecrd[0, 0]],
[ecrd[1, 1], ecrd[1, 1], ecrd[3, 1], ecrd[3, 1],
ecrd[1, 1]], **fmt_model_joint2d)
if element_labels:
ax.text(xt, yt, f'{ele_tag}', va='center', ha='center', color='red')
# 2d zeroLength, two node link element plot_model
elif (ele_classtag == EleClassTag.ZeroLength or
ele_classtag == EleClassTag.TwoNodeLink):
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = ecrd[0, 0]
yt = ecrd[0, 1]
# ax.plot([ecrd[0, 0], ecrd[2, 0], ecrd[2, 0], ecrd[0, 0],
# ecrd[0, 0]],
# [ecrd[1, 1], ecrd[1, 1], ecrd[3, 1], ecrd[3, 1],
# ecrd[1, 1]], **fmt_model_joint2d)
if element_labels:
ax.text(xt, yt, f'{ele_tag}', va='top', ha='left',
color='red')
## MASONRY PANEL 12 NODES
elif (ele_classtag == EleClassTag.MasonPan12 ):
nen= 12
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen,2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i,:]=ops.nodeCoord(ele_node_tag)
xt=(ecrd[3,0]-ecrd[0,0])/2
yt=(ecrd[9,1]-ecrd[0,1])/2
ax.plot([ecrd[9, 0], ecrd[3, 0]],
[ecrd[9,1],ecrd[3,1]], **fmt_model_truss)
ax.plot([ecrd[10, 0], ecrd[2, 0]],
[ecrd[10,1],ecrd[2,1]], **fmt_model_truss)
ax.plot([ecrd[8, 0], ecrd[4, 0]],
[ecrd[8,1],ecrd[4,1]], **fmt_model_truss)
ax.plot([ecrd[0, 0], ecrd[6, 0]],
[ecrd[0,1],ecrd[6,1]], **fmt_model_truss)
ax.plot([ecrd[1, 0], ecrd[5, 0]],
[ecrd[1,1],ecrd[5 ,1]], **fmt_model_truss)
ax.plot([ecrd[11, 0], ecrd[7, 0]],
[ecrd[11,1],ecrd[7,1]], **fmt_model_truss)
if element_labels:
va = 'center'
ha = 'center'
ax.text(xt, yt, f'Panel: {ele_tag}', va=va, ha=ha, color='red')
if node_supports:
_plot_supports(node_tags, ax)
ax.axis('equal')
ax.grid(False)
return ax
def _plot_supports(node_tags, ax):
ndim = ops.getNDM()[0]
# fix 2d support: square
verts = [(-1., 0.),
(1., 0.),
(1., -2.),
(-1., -2.),
(-1., 0.)]
codes = [Path.MOVETO,
Path.LINETO,
Path.LINETO,
Path.LINETO,
Path.CLOSEPOLY]
path_fix = Path(verts, codes)
# hinge/pin 2d support: ^
verts = [(0., 0.),
(-0.2, -0.25),
(0.2, -0.25),
(0., 0.)]
codes = [Path.MOVETO,
Path.LINETO,
Path.LINETO,
Path.CLOSEPOLY]
path_pin = Path(verts, codes)
# roller horiz 2d support: o_
verts = [(0., 0.),
(1., 0.),
(1., -1.),
(1., -2.),
(0., -2.),
(-1., -2.),
(-1., -1.),
(-1., 0.),
(0., 0.),
(-1.5, -2.),
(1.5, -2.)]
codes = [Path.MOVETO,
Path.CURVE3,
Path.LINETO,
Path.CURVE3,
Path.LINETO,
Path.CURVE3,
Path.LINETO,
Path.CURVE3,
Path.LINETO,
Path.MOVETO,
Path.LINETO]
path_roller_horiz = Path(verts, codes)
# roller vert 2d support: o|
verts = [(0., 0.),
(0., 1.),
(1., 1.),
(2., 1.),
(2., 0.),
(2., -1.),
(1., -1.),
(0., -1.),
(0., 0.),
(2, -1.5),
(2, 1.5)]
codes = [Path.MOVETO,
Path.CURVE3,
Path.LINETO,
Path.CURVE3,
Path.LINETO,
Path.CURVE3,
Path.LINETO,
Path.CURVE3,
Path.LINETO,
Path.MOVETO,
Path.LINETO]
path_roller_vert = Path(verts, codes)
m_color = 'm'
m_fstyle = 'full'
m_size = 16
fixed_nodes = ops.getFixedNodes()
for node_tag in fixed_nodes:
nd_crd = ops.nodeCoord(node_tag)
fixed_dofs = ops.getFixedDOFs(node_tag)
if ndim == 2:
if (fixed_dofs == [1, 2, 3]):
m_type = path_fix
elif (fixed_dofs == [1, 2] or fixed_dofs == [2, 1]):
m_type = path_pin
elif (fixed_dofs == [1]):
m_type = path_roller_vert
elif (fixed_dofs == [2]):
m_type = path_roller_horiz
ax.plot(nd_crd[0], nd_crd[1], marker=m_type, markersize=m_size,
color=m_color, fillstyle=m_fstyle)
elif ndim == 3:
if (fixed_dofs == [1, 2, 3, 4, 5, 6]):
m_type = path_fix
elif (fixed_dofs == [1, 2, 3]):
m_type = path_pin
else:
m_type = None
ax.plot(nd_crd[0], nd_crd[1], nd_crd[2], marker=m_type, markersize=m_size,
color=m_color, fillstyle=m_fstyle)
return ax
def _plot_model_3d(node_labels, element_labels, offset_nd_label,
axis_off, az_el, fig_wi_he, fig_lbrt, local_axes,
nodes_only, fmt_model, node_supports, gauss_points,
fmt_gauss_points, fmt_model_truss,
truss_node_offset, ax):
node_tags = ops.getNodeTags()
ele_tags = ops.getEleTags()
azim, elev = az_el
if not ax:
if fig_wi_he:
fig_wi, fig_he = fig_wi_he
fig = plt.figure(figsize=(fig_wi / 2.54, fig_he / 2.54))
# fig.subplots_adjust(left=.08, bottom=.08, right=.985, top=.94)
else:
fig = plt.figure()
if fig_lbrt:
fleft, fbottom, fright, ftop = fig_lbrt
fig.subplots_adjust(left=fleft, bottom=fbottom, right=fright, top=ftop)
ax = fig.add_subplot(111, projection=Axes3D.name)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.view_init(azim=azim, elev=elev)
max_x_crd, max_y_crd, max_z_crd, max_crd = -np.inf, -np.inf, \
-np.inf, -np.inf
# nodes only in the OpenSees domain
if len(ele_tags) == 0 or nodes_only:
for node_tag in node_tags:
ax.plot(ops.nodeCoord(node_tag)[0],
ops.nodeCoord(node_tag)[1],
ops.nodeCoord(node_tag)[2],
color='green', marker='o')
if node_labels:
ax.text(ops.nodeCoord(node_tag)[0],
ops.nodeCoord(node_tag)[1],
ops.nodeCoord(node_tag)[2],
f'{node_tag}', va='bottom', ha='left', color='blue')
return ax
for node_tag in node_tags:
x_crd = ops.nodeCoord(node_tag)[0]
y_crd = ops.nodeCoord(node_tag)[1]
z_crd = ops.nodeCoord(node_tag)[2]
if x_crd > max_x_crd:
max_x_crd = x_crd
if y_crd > max_y_crd:
max_y_crd = y_crd
if z_crd > max_z_crd:
max_z_crd = z_crd
if offset_nd_label == 0 or offset_nd_label == 0.:
_offset = 0.
else:
max_crd = np.amax([max_x_crd, max_y_crd, max_z_crd])
_offset = 0.005 * max_crd
if node_labels:
for node_tag in node_tags:
ax.text(ops.nodeCoord(node_tag)[0] + _offset,
ops.nodeCoord(node_tag)[1] + _offset,
ops.nodeCoord(node_tag)[2] + _offset,
f'{node_tag}', va='bottom', ha='left', color='blue')
for i, ele_tag in enumerate(ele_tags):
ele_classtag = ops.getEleClassTags(ele_tag)[0]
nen = np.shape(ops.eleNodes(ele_tag))[0]
if (ele_classtag == EleClassTag.ElasticBeam3d or
ele_classtag == EleClassTag.ForceBeamColumn3d or
ele_classtag == EleClassTag.DispBeamColumn3d or
ele_classtag == EleClassTag.ElasticTimoshenkoBeam3d or
ele_classtag == EleClassTag.Pipe or
ele_classtag == EleClassTag.CurvedPipe or
ele_classtag == EleClassTag.DispBeamColumn3dThermal or
ele_classtag == EleClassTag.ForceBeamColumn3dThermal):
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd_nodes = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd_nodes[i, :] = ops.nodeCoord(ele_node_tag)
ecrd_eles = np.copy(ecrd_nodes)
ele_offsets = np.array(ops.eleResponse(ele_tag, 'offsets'))
nz_offsets = np.nonzero(ele_offsets)[0] # tuple of arrays
if np.any(nz_offsets):
# modify ecrd_eles
ecrd_eles[:, 0] += ele_offsets[[0, 3]]
ecrd_eles[:, 1] += ele_offsets[[1, 4]]
ecrd_eles[:, 2] += ele_offsets[[2, 5]]
ax.plot([ecrd_nodes[0, 0], ecrd_eles[0, 0]],
[ecrd_nodes[0, 1], ecrd_eles[0, 1]],
[ecrd_nodes[0, 2], ecrd_eles[0, 2]],
**fmt_model_rigid_offset)
ax.plot([ecrd_nodes[1, 0], ecrd_eles[1, 0]],
[ecrd_nodes[1, 1], ecrd_eles[1, 1]],
[ecrd_nodes[1, 2], ecrd_eles[1, 2]],
**fmt_model_rigid_offset)
else:
pass
ax.plot(ecrd_eles[:, 0],
ecrd_eles[:, 1],
ecrd_eles[:, 2], **fmt_model)
# ax.plot(ecrd[:, 0], ecrd[:, 1], ecrd[:, 2], **fmt_model)
# location of label
xt = sum(ecrd_eles[:, 0]) / nen
yt = sum(ecrd_eles[:, 1]) / nen
zt = sum(ecrd_eles[:, 2]) / nen
# fixme: placement of node_tag labels
if element_labels:
if ecrd_eles[1, 0] - ecrd_eles[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd_eles[1, 1] - ecrd_eles[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd_eles[1, 2] - ecrd_eles[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt + offset_x,
yt + offset_y,
zt + offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
if local_axes:
# eo = Eo[i, :]
xloc = ops.eleResponse(ele_tag, 'xlocal')
yloc = ops.eleResponse(ele_tag, 'ylocal')
zloc = ops.eleResponse(ele_tag, 'zlocal')
g = np.vstack((xloc, yloc, zloc))
L = bar_length(ecrd_eles)
alen = 0.1 * L
# plot local axis directional vectors: workaround quiver = arrow
ax.quiver(xt, yt, zt, g[0, 0], g[0, 1], g[0, 2], color='g',
lw=2, length=alen, alpha=.8, normalize=True)
ax.quiver(xt, yt, zt, g[1, 0], g[1, 1], g[1, 2], color='b',
lw=2, length=alen, alpha=.8, normalize=True)
ax.quiver(xt, yt, zt, g[2, 0], g[2, 1], g[2, 2], color='r',
lw=2, length=alen, alpha=.8, normalize=True)
if gauss_points:
Lgps = ops.eleResponse(ele_tag, 'integrationPoints')
for Lgpi in Lgps:
dxi = (ecrd_eles[1, 0] - ecrd_eles[0, 0]) * Lgpi / bar_length(ecrd_eles)
dyi = (ecrd_eles[1, 1] - ecrd_eles[0, 1]) * Lgpi / bar_length(ecrd_eles)
dzi = (ecrd_eles[1, 2] - ecrd_eles[0, 2]) * Lgpi / bar_length(ecrd_eles)
ax.plot(ecrd_eles[0, 0] + dxi, ecrd_eles[0, 1] + dyi, ecrd_eles[0, 2] + dzi, **fmt_gauss_points)
elif ele_classtag in [EleClassTag.truss,EleClassTag.trussSection,EleClassTag.CorotTruss]:
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
if truss_node_offset:
Lx = ecrd[1, 0] - ecrd[0, 0]
Ly = ecrd[1, 1] - ecrd[0, 1]
Lz = ecrd[1, 2] - ecrd[0, 2]
dLx = truss_node_offset * Lx
dLy = truss_node_offset * Ly
dLz = truss_node_offset * Lz
ecrd[0, 0] += dLx
ecrd[1, 0] -= dLx
ecrd[0, 1] += dLy
ecrd[1, 1] -= dLy
ecrd[0, 2] += dLz
ecrd[1, 2] -= dLz
ax.plot(ecrd[:, 0], ecrd[:, 1], ecrd[:, 2], **fmt_model_truss)
# fixme: placement of node_tag labels
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# quad in 3d
# elif nen == 4:
elif ele_classtag == EleClassTag.FourNodeTetrahedron:
nen = 4
nodes_geo_order = [0, 1, 2, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1],
ecrd[nodes_geo_order, 2],
**fmt_model)
for j in range(3):
ax.plot(ecrd[[j, 3], 0],
ecrd[[j, 3], 1],
ecrd[[j, 3], 2], **fmt_model)
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# if node_labels:
# for node_tag in node_tags:
# ax.text(ops.nodeCoord(node_tag)[0]+_offset,
# ops.nodeCoord(node_tag)[1]+_offset,
# ops.nodeCoord(node_tag)[2]+_offset,
# f'{node_tag}', va='bottom', ha='left', color='blue')
elif (ele_classtag in {EleClassTag.TenNodeTetrahedron,
EleClassTag.TenNodeTetrahedronSK}):
nen = 10
nodes_geo_order = [0, 4, 1, 5, 2, 6, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1],
ecrd[nodes_geo_order, 2],
'b.-')
for j in range(3):
ax.plot(ecrd[[j, 7+j, 3], 0],
ecrd[[j, 7+j, 3], 1],
ecrd[[j, 7+j, 3], 2], 'b.-')
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# if node_labels:
# for node_tag in node_tags:
# ax.text(ops.nodeCoord(node_tag)[0]+_offset,
# ops.nodeCoord(node_tag)[1]+_offset,
# ops.nodeCoord(node_tag)[2]+_offset,
# f'{node_tag}', va='bottom', ha='left', color='blue')
elif (ele_classtag == EleClassTag.ShellMITC4 or
ele_classtag == EleClassTag.ASDShellQ4 or
ele_classtag == EleClassTag.ShellNLDKGQ or
ele_classtag == EleClassTag.ShellDKGQ or
ele_classtag == EleClassTag.quad4n):
nen = 4
nodes_geo_order = [0, 1, 2, 3, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1],
ecrd[nodes_geo_order, 2],
**fmt_model)
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# if node_labels:
# for node_tag in node_tags:
# ax.text(ops.nodeCoord(node_tag)[0]+_offset,
# ops.nodeCoord(node_tag)[1]+_offset,
# ops.nodeCoord(node_tag)[2]+_offset,
# f'{node_tag}', va='bottom', ha='left', color='blue')
elif (ele_classtag == EleClassTag.ShellNLDKGT or
ele_classtag == EleClassTag.ShellDKGT):
nen = 3
nodes_geo_order = [0, 1, 2, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1],
ecrd[nodes_geo_order, 2],
**fmt_model)
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# if node_labels:
# for node_tag in node_tags:
# ax.text(ops.nodeCoord(node_tag)[0]+_offset,
# ops.nodeCoord(node_tag)[1]+_offset,
# ops.nodeCoord(node_tag)[2]+_offset,
# f'{node_tag}', va='bottom', ha='left', color='blue')
# mvlem 3d plot model
elif (ele_classtag == EleClassTag.MVLEM_3D
or ele_classtag == EleClassTag.SFI_MVLEM_3D):
nen = 4
nodes_geo_order = [0, 1, 3, 2, 0]
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1],
ecrd[nodes_geo_order, 2],
**fmt_model)
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# if node_labels:
# for node_tag in node_tags:
# ax.text(ops.nodeCoord(node_tag)[0]+_offset,
# ops.nodeCoord(node_tag)[1]+_offset,
# ops.nodeCoord(node_tag)[2]+_offset,
# f'{node_tag}', va='bottom', ha='left', color='blue')
# 8-node brick, 3d model
# elif nen == 8:
elif (ele_classtag == EleClassTag.brick8n or
ele_classtag == EleClassTag.SSPbrick):
nen = 8
nodes_geo_order = np.array([0, 1, 2, 3, 0], dtype=int) # bottom face nodes loop
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1],
ecrd[nodes_geo_order, 2],
**fmt_model)
ax.plot(ecrd[nodes_geo_order + 4, 0],
ecrd[nodes_geo_order + 4, 1],
ecrd[nodes_geo_order + 4, 2],
**fmt_model)
for j in range(4):
ax.plot(ecrd[[j, j+4], 0],
ecrd[[j, j+4], 1],
ecrd[[j, j+4], 2], **fmt_model)
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# if node_labels:
# for node_tag in node_tags:
# ax.text(ops.nodeCoord(node_tag)[0]+_offset,
# ops.nodeCoord(node_tag)[1]+_offset,
# ops.nodeCoord(node_tag)[2]+_offset,
# f'{node_tag}', va='bottom', ha='left', color='blue')
# 20-node brick, 3d model
# elif nen == 20:
elif (ele_classtag == EleClassTag.brick20n):
nen = 20
nodes_geo_order = np.array([0, 8, 1, 9, 2, 10, 3, 11, 0], dtype=int) # bottom face nodes loop
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
ax.plot(ecrd[nodes_geo_order, 0],
ecrd[nodes_geo_order, 1],
ecrd[nodes_geo_order, 2],
**fmt_model, mfc='g', mec='g')
ax.plot(ecrd[nodes_geo_order + 4, 0],
ecrd[nodes_geo_order + 4, 1],
ecrd[nodes_geo_order + 4, 2],
**fmt_model, mfc='g', mec='g')
# for j in range(2):
# ax.plot(ecrd[[0+j*4, 8+j*4, 1+j*4, 9+j*4, 2+j*4, 10+j*4, 3+j*4, 11+j*4, 0+j*4], 0],
# ecrd[[0+j*4, 8+j*4, 1+j*4, 9+j*4, 2+j*4, 10+j*4, 3+j*4, 11+j*4, 0+j*4], 1],
# ecrd[[0+j*4, 8+j*4, 1+j*4, 9+j*4, 2+j*4, 10+j*4, 3+j*4, 11+j*4, 0+j*4], 2],
# **fmt_model, mfc='g', mec='g')
for j in range(4):
ax.plot(ecrd[[j, 16+j, 4+j], 0],
ecrd[[j, 16+j, 4+j], 1],
ecrd[[j, 16+j, 4+j], 2], **fmt_model,
mfc='g', mec='g')
# location of label
xt = sum(ecrd[:, 0]) / nen
yt = sum(ecrd[:, 1]) / nen
zt = sum(ecrd[:, 2]) / nen
if element_labels:
if ecrd[1, 0] - ecrd[0, 0] == 0:
va = 'center'
ha = 'left'
offset_x, offset_y, offset_z = _offset, 0.0, 0.0
elif ecrd[1, 1] - ecrd[0, 1] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, _offset, 0.0
elif ecrd[1, 2] - ecrd[0, 2] == 0:
va = 'bottom'
ha = 'center'
offset_x, offset_y, offset_z = 0.0, 0.0, _offset
else:
va = 'bottom'
ha = 'left'
offset_x, offset_y, offset_z = 0.03, 0.03, 0.03
ax.text(xt+offset_x, yt+offset_y, zt+offset_z, f'{ele_tag}',
va=va, ha=ha, color='red')
# if node_labels:
# for node_tag in node_tags:
# ax.text(ops.nodeCoord(node_tag)[0]+_offset,
# ops.nodeCoord(node_tag)[1]+_offset,
# ops.nodeCoord(node_tag)[2]+_offset,
# f'{node_tag}', va='bottom', ha='left', color='blue')
# 3d zeroLength, two node link element plot_model
elif (ele_classtag == EleClassTag.ZeroLength or
ele_classtag == EleClassTag.TwoNodeLink):
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# location of label
xt = ecrd[0, 0]
yt = ecrd[0, 1]
zt = ecrd[0, 2]
# ax.plot([ecrd[0, 0], ecrd[2, 0], ecrd[2, 0], ecrd[0, 0],
# ecrd[0, 0]],
# [ecrd[1, 1], ecrd[1, 1], ecrd[3, 1], ecrd[3, 1],
# ecrd[1, 1]], **fmt_model_joint2d)
if element_labels:
ax.text(xt, yt, zt, f'{ele_tag}', va='top', ha='left',
color='red')
if node_supports:
_plot_supports(node_tags, ax)
ax.set_box_aspect((np.ptp(ax.get_xlim3d()),
np.ptp(ax.get_ylim3d()),
np.ptp(ax.get_zlim3d())))
return ax
[docs]
def plot_model(node_labels=1, element_labels=1, offset_nd_label=False,
axis_off=0, az_el=az_el, fig_wi_he=False,
fig_lbrt=False, local_axes=True, nodes_only=False,
fmt_model=fmt_model,
fmt_model_nodes_only=fmt_model_nodes_only,
node_supports=True, gauss_points=True,
fmt_gauss_points=fmt_gauss_points,
fmt_model_truss=fmt_model_truss,
truss_node_offset=0.96, ax=False):
"""Plot defined model of the structure.
Args:
node_labels (int): 1 - plot node labels, 0 - do not plot them;
(default: 1)
element_labels (int): 1 - plot element labels, 0 - do not plot
them; (default: 1)
offset_nd_label (bool): False - do not offset node labels from the
actual node location. This option can enhance visibility.
axis_off (int): 0 - turn off axes, 1 - display axes; (default: 0)
az_el (tuple): contains azimuth and elevation for 3d plots. For 2d
plots this parameter is neglected.
fig_wi_he (tuple): contains width and height of the figure
fig_lbrt (tuple): a tuple contating left, bottom, right and top offsets
local_axes (bool): True - show cross section local axes or False.
The green, red and blue arrows denote the element axis direction,
the z-local axis and the y-local axis.
nodes_only (bool): True - show the nodes only, although the elements
are defined. Default: False.
fmt_model (dict): A dictionary containing formatting the line and markers
of the model elements. The formatting options can be: linewidth, color,
marker, markersize. See matplotlib.plot documentation for more details,
if necessary.
node_supports (bool): True - show the supports. Default: True.
gauss_points (bool): True - show the integration (Gauss) points. Default: True.
fmt_gauss_points (dict): A dictionary containing formatting the marker
of the gauss point.
truss_node_offset (float): If non-zero, the nodes are offset to show pin markers.
The number should preferably be between 0.90-0.97.
Zero (0) or False means no offset. Default: 0.96.
ax: axis object.
Usage:
``plot_model()`` - plot model with node and element labels.
``plot_model(node_labels=0, element_labels=0)`` - plot model without node
element labels
``plot_model(fig_wi_he=(20., 14.))`` - plot model in a window 20 cm long,
and 14 cm high.
``plot_model(nodes_only=True)`` - plot only the nodes even thought the elements
are defined.
``plot_model(node_supports=False)`` - plot the model without the supports.
"""
# az_el - azimut, elevation used for 3d plots only
node_tags = ops.getNodeTags()
ndim = ops.getNDM()[0]
if ndim == 2:
ax = _plot_model_2d(node_labels, element_labels,
offset_nd_label, axis_off, fig_wi_he,
fig_lbrt, nodes_only, fmt_model,
fmt_model_nodes_only, node_supports,
gauss_points, fmt_gauss_points,
fmt_model_truss, truss_node_offset, ax)
if axis_off:
ax.axis('off')
elif ndim == 3:
ax = _plot_model_3d(node_labels, element_labels,
offset_nd_label, axis_off, az_el,
fig_wi_he, fig_lbrt, local_axes,
nodes_only, fmt_model, node_supports,
gauss_points, fmt_gauss_points,
fmt_model_truss, truss_node_offset, ax)
if axis_off:
ax.axis('off')
else:
print(f'\nWarning! ndim: {ndim} not supported yet.')
# plt.show() # call this from main py file for more control
return ax
def plot_supports_and_loads_2d(nep=17):
"""This function is removed. Use plot_model() and
plot_loads_2d() for showing support conditions and loads,
respectively.
"""
print('\nWarning! plot_supports_and_loads_2d has been removed.') # noqa: E501
print('\nWarning! Use plot_model(node_supports=True) for showing supports and') # noqa: E501
print('\nWarning! plot_loads_2d() for showing loads') # noqa: E501
def plot_loads_2d(nep, sfac, fig_wi_he, fig_lbrt, fmt_model_loads,
node_supports, truss_node_offset, ax):
"""Display the nodal and element loads applied to the 2d models.
This is a local function: use plot_load()
"""
if not ax:
if fig_wi_he:
fig_wi, fig_he = fig_wi_he
fig, ax = plt.subplots(figsize=(fig_wi / 2.54, fig_he / 2.54))
else:
fig, ax = plt.subplots()
if fig_lbrt:
fleft, fbottom, fright, ftop = fig_lbrt
fig.subplots_adjust(left=fleft, bottom=fbottom, right=fright, top=ftop)
ax = plot_model(node_labels=0, element_labels=0, fmt_model=fmt_model_loads,
node_supports=node_supports,
truss_node_offset=truss_node_offset, ax=ax)
if not sfac:
ratio = 0.1
min_x, max_x = ax.get_xlim()
min_y, max_y = ax.get_ylim()
xsfac = ratio * abs(max_x - min_x)
ysfac = ratio * abs(max_y - min_y)
sfac = max(xsfac, ysfac)
node_tags = ops.getNodeTags()
ele_tags = ops.getEleTags()
Wx = False
waa = False
wab = False
Ew = get_Ew_data_from_ops_domain()
for ele_tag in ele_tags:
ele_class_tag = ops.getEleClassTags(ele_tag)[0]
if (ele_class_tag == EleClassTag.ElasticBeam2d
or ele_class_tag == EleClassTag.ForceBeamColumn2d
or ele_class_tag == EleClassTag.DispBeamColumn2d
or ele_class_tag == EleClassTag.TimoshenkoBeamColumn2d
or ele_class_tag == EleClassTag.ElasticTimoshenkoBeam2d
or ele_class_tag == EleClassTag.DispBeamColumn2dThermal
or ele_class_tag == EleClassTag.ForceBeamColumn2dThermal):
# by default no element load
ele_load_data = []
if ele_tag in Ew:
ele_load_data = Ew[ele_tag]
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd_nodes = np.zeros((nen, 2))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd_nodes[i, :] = ops.nodeCoord(ele_node_tag)
ecrd_eles = np.copy(ecrd_nodes)
ele_offsets = np.array(ops.eleResponse(ele_tag, 'offsets'))
nz_offsets = np.nonzero(ele_offsets)[0] # tuple of arrays
if np.any(nz_offsets):
# modify ecrd_eles
ecrd_eles[:, 0] += ele_offsets[[0, 3]]
ecrd_eles[:, 1] += ele_offsets[[1, 4]]
else:
pass
Lxy = ecrd_eles[1, :] - ecrd_eles[0, :]
L = np.sqrt(Lxy @ Lxy)
cosa, cosb = Lxy / L
xl = np.linspace(0., L, nep)
xl2 = np.linspace(0., L, 5)
dxs = np.array([0.1, 0.2, 0.3, 0.4])
xl3 = np.append(0., np.cumsum(dxs)) * L
one = np.ones(nep)
for ele_load_data_i in ele_load_data:
ele_load_type = ele_load_data_i[0]
if ele_load_type == '-beamPoint':
Pt, aL, Pa = ele_load_data_i[1:4]
a = aL * L
s = sfac * np.sign(Pt)
s_0 = np.zeros(2)
s_0[0] = ecrd_eles[0, 0] + a * cosa
s_0[1] = ecrd_eles[0, 1] + a * cosb
s_p = np.copy(s_0)
s_p[0] -= s * cosb
s_p[1] += s * cosa
ax.arrow(s_0[0], s_0[1],
s_p[0] - s_0[0], s_p[1] - s_0[1],
head_starts_at_zero=True, # default False
head_width=0.1 * sfac, head_length=0.2 * sfac,
fc='b', ec='b',
length_includes_head=False, shape='full')
ax.text(sum(ecrd_eles[:, 0])/2, sum(ecrd_eles[:, 1])/2, f'{abs(Pt)}', color='b')
elif ele_load_type == '-beamUniform':
n_ele_load_data = len(ele_load_data_i)
if n_ele_load_data == 3:
Wy, Wx = ele_load_data_i[1], ele_load_data_i[2]
text_string = f'q = {abs(Wy):.3g}, {abs(Wx):.3g}'
s = sfac * one * np.sign(Wy)
elif n_ele_load_data == 7:
wta, waa, aL, bL, wtb, wab = ele_load_data_i[1:7]
text_string = f'q = {abs(wta)}, {abs(waa)}, {aL}, {bL}, {abs(wtb)}, {abs(wab)}'
a, b = aL * L, bL * L
bma = b - a
s = np.zeros(nep)
indx = 0
for x in np.nditer(xl):
xma = x - a
wtx = wta + (wtb - wta) * xma / bma
if x < a:
pass
elif x >= a and x <= b:
s[indx] = wtx
elif x > b:
pass
indx += 1
s_0 = np.zeros((nep, 2))
s_0[0, :] = [ecrd_eles[0, 0], ecrd_eles[0, 1]]
s_0[1:, 0] = s_0[0, 0] + xl[1:] * cosa
s_0[1:, 1] = s_0[0, 1] + xl[1:] * cosb
s_p = np.copy(s_0)
s_p[:, 0] += s * cosb
s_p[:, 1] -= s * cosa
# for i in np.arange(nep):
# ax.arrow(s_0[i, 0], s_0[i, 1],
# s_p[i, 0] - s_0[i, 0], s_p[i, 1] - s_0[i, 1],
# head_width=0.1*sfac, head_length=0.2*sfac,
# head_starts_at_zero=True, # default False
# fc='r', ec='r',
# length_includes_head=True, shape='full')
for i in np.arange(nep):
ax.arrow(s_p[i, 0], s_p[i, 1],
s_0[i, 0] - s_p[i, 0], s_0[i, 1] - s_p[i, 1],
# width = 0.005,
lw=1,
head_width=0.1 * sfac, head_length=0.2 * sfac,
# head_starts_at_zero=True, # default False
# overhang=0.5,
fc='r', ec='r',
length_includes_head=True, shape='full',
joinstyle='round')
plt.plot([s_p[0, 0], s_p[-1, 0]],
[s_p[0, 1], s_p[-1, 1]], 'r', lw=1)
ax.text(s_p[int((nep - 1) / 3), 0],
s_p[int((nep - 1) / 3), 1],
text_string, va='bottom', ha='center', color='r')
if Wx != 0:
if Wx < 0:
ax.quiver(s_0[:-1, 0], s_0[:-1, 1],
-1 * (s_0[1:, 0] - s_0[:-1, 0]),
-1 * (s_0[1:, 1] - s_0[:-1, 1]),
scale_units='xy', angles='xy', scale=0.8,
color='m')
else:
ax.quiver(s_0[:-1, 0], s_0[:-1, 1],
s_0[1:, 0] - s_0[:-1, 0],
s_0[1:, 1] - s_0[:-1, 1],
scale_units='xy', angles='xy', scale=0.8,
color='m')
if waa != 0 or wab != 0:
sa = np.zeros((5, 2))
sa[0, :] = [ecrd_eles[0, 0], ecrd_eles[0, 1]]
sa[1:, 0] = sa[0, 0] + xl3[1:] * cosa
sa[1:, 1] = sa[0, 1] + xl3[1:] * cosb
for i, xl in enumerate(xl3[:-1]):
ax.arrow(sa[i, 0], sa[i, 1],
sa[i+1, 0]-sa[i, 0], sa[i+1, 1]-sa[i, 1],
width=0.01*L,
fc='m', ec='m',
length_includes_head=True, shape='full')
for node_tag in node_tags:
nd_crd = ops.nodeCoord(node_tag)
nodal_loads = ops.nodeUnbalance(node_tag)
nodal_loads_idx = np.nonzero(nodal_loads)
if nodal_loads_idx[0].size:
for kier in nodal_loads_idx[0]:
if kier == 0 or kier == 1:
kolor = 'b'
if kier == 0:
kier2 = np.sign(nodal_loads[kier])
dx, dy = sfac * np.sign(kier2), 0.
elif kier == 1:
kier2 = np.sign(nodal_loads[kier])
dx, dy = 0., sfac * np.sign(kier2)
ax.arrow(nd_crd[0] - dx, nd_crd[1] - dy,
dx, dy,
lw=3,
head_width=0.1 * sfac, head_length=0.2 * sfac,
fc=kolor, ec=kolor,
length_includes_head=True, shape='full',
joinstyle='round')
ax.text(nd_crd[0] - dx, nd_crd[1] - dy,
f' {abs(nodal_loads[kier]):.5g}', color=kolor)
elif kier == 2:
kolor = 'r'
kier2 = np.sign(nodal_loads[kier])
nodal_load_str = f'{abs(nodal_loads[kier]):.5g}'
if kier2 > 0:
marker_type = r'$\curvearrowleft$'
ax.text(nd_crd[0], nd_crd[1], f'\n {nodal_load_str}',
color=kolor, va='top', ha='right')
elif kier2 < 0:
marker_type = r'$\curvearrowright$'
ax.text(nd_crd[0], nd_crd[1], f'\n {nodal_load_str}',
color=kolor, va='top', ha='left')
ax.plot(nd_crd[0], nd_crd[1], marker=marker_type, markersize=30,
color=kolor)
ax.axis('equal')
ax.grid(False)
return ax
def plot_loads_3d(nep, sfac, fig_wi_he, fig_lbrt, fmt_model_loads,
node_supports, truss_node_offset, local_axes, ax):
"""Display the nodal and element loads applied to the 3d models.
This is a local function: use plot_load()
"""
if not ax:
if fig_wi_he:
fig_wi, fig_he = fig_wi_he
fig = plt.figure(figsize=(fig_wi / 2.54, fig_he / 2.54))
# fig.subplots_adjust(left=.08, bottom=.08, right=.985, top=.94)
else:
fig = plt.figure()
if fig_lbrt:
fleft, fbottom, fright, ftop = fig_lbrt
fig.subplots_adjust(left=fleft, bottom=fbottom, right=fright, top=ftop)
ax = fig.add_subplot(111, projection=Axes3D.name)
ax = plot_model(node_labels=0, element_labels=0, fmt_model=fmt_model_loads,
node_supports=node_supports, local_axes=local_axes, ax=ax)
if not sfac:
ratio = 0.15 # initial scale 10% of the max dimension
min_x, max_x = ax.get_xlim()
min_y, max_y = ax.get_ylim()
min_z, max_z = ax.get_zlim()
xsfac = ratio * abs(max_x - min_x)
ysfac = ratio * abs(max_y - min_y)
zsfac = ratio * abs(max_z - min_z)
sfac = max(xsfac, ysfac, zsfac)
node_tags = ops.getNodeTags()
ele_tags = ops.getEleTags()
# ndf = ops.getNDF()[0]
Wx = False
waa = False
wab = False
Ew = get_Ew_data_from_ops_domain_3d()
for ele_tag in ele_tags: # plot load 3d
ele_class_tag = ops.getEleClassTags(ele_tag)[0]
if (ele_class_tag == EleClassTag.ElasticBeam3d
or ele_class_tag == EleClassTag.ForceBeamColumn3d
or ele_class_tag == EleClassTag.DispBeamColumn3d
or ele_class_tag == EleClassTag.DispBeamColumn3dThermal
or ele_class_tag == EleClassTag.ForceBeamColumn3dThermal):
# by default no element load
ele_load_data = []
if ele_tag in Ew:
ele_load_data = Ew[ele_tag]
nen = 2
ele_node_tags = ops.eleNodes(ele_tag)
ecrd_nodes = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd_nodes[i, :] = ops.nodeCoord(ele_node_tag)
ecrd_eles = np.copy(ecrd_nodes)
ele_offsets = np.array(ops.eleResponse(ele_tag, 'offsets'))
nz_offsets = np.nonzero(ele_offsets)[0] # tuple of arrays
if np.any(nz_offsets): # plot load 3d
# modify ecrd_eles
ecrd_eles[:, 0] += ele_offsets[[0, 3]]
ecrd_eles[:, 1] += ele_offsets[[1, 4]]
ecrd_eles[:, 2] += ele_offsets[[2, 5]]
else:
pass
# step 2
xloc = ops.eleResponse(ele_tag, 'xlocal')
yloc = ops.eleResponse(ele_tag, 'ylocal')
zloc = ops.eleResponse(ele_tag, 'zlocal')
g = np.vstack((xloc, yloc, zloc))
G, L = rot_transf_3d(ecrd_eles, g)
xl = np.linspace(0., L, nep)
xl2 = np.linspace(0., L, 5)
# xl2 = np.logspace(0., 1., nep) * L / nep
dxs = np.array([0.1, 0.2, 0.3, 0.4])
# dxs = np.array([0.15, 0.2, 0.275, 0.375])
xl3 = np.append(0., np.cumsum(dxs)) * L
one = np.ones(nep)
for ele_load_data_i in ele_load_data: # plot load 3d
ele_load_type = ele_load_data_i[0]
# fixme uncomment when regular uniformLoad works! (1)
if ele_load_type == '-beamPoint': # plot load 3d
Py, Pz, aL, Px = ele_load_data_i[1:5]
Pxyz = [Px, Py, Pz]
xa = aL * L
text_string = f'Pyzx:{Py},{Pz},{Px}'
s_0 = np.zeros(3)
s_0[0] = ecrd_eles[0, 0] + xa * g[0, 0]
s_0[1] = ecrd_eles[0, 1] + xa * g[0, 1]
s_0[2] = ecrd_eles[0, 2] + xa * g[0, 2]
for k, Pi in enumerate(Pxyz):
if Pi != 0:
dir_plt = k
s = sfac * np.sign(Pi)
dx = s * g[dir_plt, 0]
dy = s * g[dir_plt, 1]
dz = s * g[dir_plt, 2]
# plot arrows: 3d beamPoint
nn = 1.0 # additional scaling for arrows
ax.arrow3D(s_0[0], s_0[1], s_0[2],
nn * dx,
nn * dy,
nn * dz, mutation_scale=8,
ec='b', fc='c')
ax.text(s_0[0] + dx, s_0[1] + dy, s_0[2] + dz,
f'{abs(Pi)}', va='bottom', ha='center', color='r')
elif ele_load_type == '-beamUniform': # plot load 3d
n_ele_load_data = len(ele_load_data_i)
# constant uniform element load
if n_ele_load_data == 4:
# eload_type, Wy, Wx = ele_load_data[0], ele_load_data[1], ele_load_data[2]
Wy, Wz, Wx = ele_load_data_i[1:4]
text_string = f'q = {Wy}, {Wz}, {Wx}'
Wxyz = [Wx, Wy, Wz]
s_0 = np.zeros((nep, 3))
s_0[0, :] = [ecrd_eles[0, 0], ecrd_eles[0, 1], ecrd_eles[0, 2]]
s_0[1:, 0] = s_0[0, 0] + xl[1:] * g[0, 0]
s_0[1:, 1] = s_0[0, 1] + xl[1:] * g[0, 1]
s_0[1:, 2] = s_0[0, 2] + xl[1:] * g[0, 2]
for k, Wi in enumerate(Wxyz):
if Wi != 0:
dir_plt = k
s = sfac * one * np.sign(Wi) # fixme: 2d is different
s_p = np.copy(s_0)
s_p[:, 0] += s * g[dir_plt, 0]
s_p[:, 1] += s * g[dir_plt, 1]
s_p[:, 2] += s * g[dir_plt, 2]
if k == 0:
nn = 0.3
else:
nn = 1
# plot arrows 3D beamUniform
for i in np.arange(nep):
ax.arrow3D(s_0[i, 0], s_0[i, 1], s_0[i, 2],
nn * (s_p[i, 0] - s_0[i, 0]),
nn * (s_p[i, 1] - s_0[i, 1]),
nn * (s_p[i, 2] - s_0[i, 2]), mutation_scale=8)
ax.text(s_p[5, 0], s_p[5, 1], s_p[5, 2],
f'{abs(Wi)}', va='bottom', ha='center', color='r')
# triangular or trapezoidal element load
elif n_ele_load_data == 7:
wta, waa, aL, bL, wtb, wab = ele_load_data_i[1:7]
text_string = f'q = {wta}, {waa}, {aL}, {bL}, {wtb}, {wab}'
a, b = aL * L, bL * L
bma = b - a
s = np.zeros(nep)
indx = 0
for x in np.nditer(xl):
xma = x - a
wtx = wta + (wtb - wta) * xma / bma
if x < a:
pass
elif x >= a and x <= b:
s[indx] = wtx
elif x > b:
pass
indx += 1
for node_tag in node_tags:
nd_crd = ops.nodeCoord(node_tag)
nodal_loads = ops.nodeUnbalance(node_tag) # plot load 3d
nodal_loads_idx = np.nonzero(nodal_loads)
if nodal_loads_idx[0].size:
for kier in nodal_loads_idx[0]:
dx, dy, dz = 0., 0., 0.
# horizontal or vertical nodal force (load)
if kier == 0 or kier == 1 or kier == 2:
kolor = 'b'
if kier == 0:
kier2 = np.sign(nodal_loads[kier])
if kier2 > 0:
pos_or_neg = '+'
elif kier2 < 0:
pos_or_neg = '-'
dx = sfac * np.sign(kier2)
elif kier == 1:
kier2 = np.sign(nodal_loads[kier])
if kier2 > 0:
pos_or_neg = '+'
elif kier2 < 0:
pos_or_neg = '-'
dy = sfac * np.sign(kier2)
elif kier == 2:
kier2 = np.sign(nodal_loads[kier])
if kier2 > 0:
pos_or_neg = '+'
elif kier2 < 0:
pos_or_neg = '-'
dz = sfac * np.sign(kier2)
elif kier == 3 or kier == 4 or kier == 5:
kolor = 'g'
if kier == 3:
kier2 = np.sign(nodal_loads[kier])
if kier2 > 0:
pos_or_neg = '+'
elif kier2 < 0:
pos_or_neg = '-'
dx = sfac * np.sign(kier2)
elif kier == 4:
kier2 = np.sign(nodal_loads[kier])
if kier2 > 0:
pos_or_neg = '+'
elif kier2 < 0:
pos_or_neg = '-'
dy = sfac * np.sign(kier2)
elif kier == 5:
kier2 = np.sign(nodal_loads[kier])
if kier2 > 0:
pos_or_neg = '+'
elif kier2 < 0:
pos_or_neg = '-'
dz = sfac * np.sign(kier2)
nn = 1.0 # additional scaling for arrows
ax.arrow3D(*nd_crd,
nn * dx,
nn * dy,
nn * dz,
mutation_scale=10, ec=kolor, fc=kolor)
if kier == 3 or kier == 4 or kier == 5:
nn = 0.5 # additional scaling for arrows
ax.arrow3D(nd_crd[0] + dx/5, nd_crd[1] + dy/5, nd_crd[2] + dz/5,
nn * dx,
nn * dy,
nn * dz,
mutation_scale=10, ec=kolor, fc=kolor)
ax.text(nd_crd[0] + dx, nd_crd[1] + dy, nd_crd[2] + dz,
f'{abs(nodal_loads[kier])}', va='bottom', ha='center', color='r')
return ax
[docs]
def plot_load(nep=11, sfac=False, fig_wi_he=False,
fig_lbrt=False, fmt_model_loads=fmt_model_loads,
node_supports=True, truss_node_offset=0, local_axes=False, ax=False):
"""Display the nodal and element loads applied to the 2d and 3d models.
Args:
nep (int): number of arrows when displacing element distributed loads
(default: 11)
node_supports (bool): True - show the node support conditions.
Default: False.
"""
ndim = ops.getNDM()[0]
if ndim == 2:
ax = plot_loads_2d(nep, sfac, fig_wi_he,
fig_lbrt, fmt_model_loads,
node_supports, truss_node_offset, ax)
elif ndim == 3:
ax = plot_loads_3d(nep, sfac, fig_wi_he,
fig_lbrt, fmt_model_loads,
node_supports, truss_node_offset,
local_axes, ax)
# elif ndim == 1:
# ax = plot_loads_1d()
else:
print(f'\nWarning! ndim: {ndim} not supported yet.')
return ax
def get_nodal_loads_from_ops_domain():
load_at_nodes = {}
ibeg, iend = 0, 0
node_load_tags_all_patterns = ops.getNodeLoadTags()
node_load_data_all_patterns = ops.getNodeLoadData()
for node_load_tag in node_load_tags_all_patterns:
load_at_nodes[node_load_tag] = []
for node_load_tag in node_load_tags_all_patterns:
ndf = ops.getNDF(node_load_tag)[0]
iend = ibeg+ndf
node_load_data = node_load_data_all_patterns[ibeg:iend]
ibeg = iend
load_at_nodes[node_load_tag].append(node_load_data)
return load_at_nodes
def get_Ew_data_from_ops_domain():
Ew = {}
ibeg = 0
iend = 0
ele_load_tags_all_patterns = ops.getEleLoadTags()
ele_load_types_all_patterns = ops.getEleLoadClassTags()
ele_load_data_all_patterns = ops.getEleLoadData()
for ele_load_tag in ele_load_tags_all_patterns:
Ew[ele_load_tag] = []
# iterate over ele_load_tags/classtags
for ele_load_type, ele_load_tag in zip(ele_load_types_all_patterns,
ele_load_tags_all_patterns):
# -beamUniform
if ele_load_type == LoadTag.Beam2dUniformLoad:
iend = ibeg+LoadTag.Beam2dUniformLoad_ndata
ele_load_data = ele_load_data_all_patterns[ibeg:iend]
ibeg = iend
wy, wx = ele_load_data[0], ele_load_data[1]
Ew[ele_load_tag].append(['-beamUniform', wy, wx])
# -beamUniform partial
elif ele_load_type == LoadTag.Beam2dPartialUniformLoad:
iend = ibeg+LoadTag.Beam2dPartialUniformLoad_ndata
ele_load_data = ele_load_data_all_patterns[ibeg:iend]
ibeg = iend
wta, wtb, waa, wab, aL, bL = ele_load_data
Ew[ele_load_tag].append(['-beamUniform', wta, waa, aL, bL, wtb, wab])
# -beamPoint
elif ele_load_type == LoadTag.Beam2dPointLoad:
iend = ibeg+LoadTag.Beam2dPointLoad_ndata
ele_load_data = ele_load_data_all_patterns[ibeg:iend]
ibeg = iend
Pt, Pa, aL = ele_load_data
Ew[ele_load_tag].append(['-beamPoint', Pt, aL, Pa])
else:
print(f'\nWarning! ele_load_type:\n{ele_load_type} - Unknown element load Error')
return Ew
def get_Ew_data_from_ops_domain_3d():
Ew = {}
ibeg = 0
iend = 0
ele_load_tags_all_patterns = ops.getEleLoadTags()
ele_load_types_all_patterns = ops.getEleLoadClassTags()
ele_load_data_all_patterns = ops.getEleLoadData()
for ele_load_tag in ele_load_tags_all_patterns:
Ew[ele_load_tag] = []
# iterate over ele_load_tags/classtags
for ele_load_type, ele_load_tag in zip(ele_load_types_all_patterns,
ele_load_tags_all_patterns):
# -beamUniform
if ele_load_type == LoadTag.Beam3dUniformLoad:
iend = ibeg+LoadTag.Beam3dUniformLoad_ndata
ele_load_data = ele_load_data_all_patterns[ibeg:iend]
ibeg = iend
wy, wz, wx = ele_load_data[0], ele_load_data[1], ele_load_data[2]
Ew[ele_load_tag].append(['-beamUniform', wy, wz, wx])
# -beamUniform partial
elif ele_load_type == LoadTag.Beam3dPartialUniformLoad:
iend = ibeg+LoadTag.Beam3dPartialUniformLoad_ndata
ele_load_data = ele_load_data_all_patterns[ibeg:iend]
ibeg = iend
wy, wz, wa, aL, bL = ele_load_data
Ew[ele_load_tag].append(['-beamUniform', wy, wz, wa, aL, bL])
# -beamPoint
elif ele_load_type == LoadTag.Beam3dPointLoad:
iend = ibeg+LoadTag.Beam3dPointLoad_ndata
ele_load_data = ele_load_data_all_patterns[ibeg:iend]
ibeg = iend
Py, Pz, Px, aL = ele_load_data
Ew[ele_load_tag].append(['-beamPoint', Py, Pz, aL, Px])
else:
print(f'\nWarning! ele_load_type:\n{ele_load_type} - Unknown element load Error')
return Ew
def plot_reactions_2d(sfac, fig_wi_he,
fig_lbrt, fmt_model_loads,
truss_node_offset, ax):
"""Display support reactions of the 2d models.
This is a local function: use plot_reactions()
"""
if not ax:
if fig_wi_he:
fig_wi, fig_he = fig_wi_he
fig, ax = plt.subplots(figsize=(fig_wi / 2.54, fig_he / 2.54))
else:
fig, ax = plt.subplots()
if fig_lbrt:
fleft, fbottom, fright, ftop = fig_lbrt
fig.subplots_adjust(left=fleft, bottom=fbottom, right=fright, top=ftop)
ax = plot_model(node_labels=0, element_labels=0, fmt_model=fmt_model_loads,
node_supports=False,
truss_node_offset=truss_node_offset, ax=ax)
if not sfac:
ratio = 0.1
min_x, max_x = ax.get_xlim()
min_y, max_y = ax.get_ylim()
xsfac = ratio * abs(max_x - min_x)
ysfac = ratio * abs(max_y - min_y)
sfac = max(xsfac, ysfac)
fixed_nodes = ops.getFixedNodes()
for node_tag in fixed_nodes:
nd_crd = ops.nodeCoord(node_tag)
fixed_dofs = ops.getFixedDOFs(node_tag)
for dof in fixed_dofs:
reaction = ops.nodeReaction(node_tag, dof)
if not isclose(reaction, 0., abs_tol=1e-9):
if dof in [1, 2]:
kolor = 'b'
reac_sign = np.sign(reaction)
if dof == 1:
dx = sfac * reac_sign
dy = 0.
ha_my = 'center'
elif dof == 2:
reac_sign = np.sign(reaction)
dx = 0.
dy = sfac * reac_sign
ha_my = 'left'
if reac_sign > 0:
x0, y0 = nd_crd[0] - dx, nd_crd[1] - dy
x0t, y0t = x0, y0
else:
x0, y0 = nd_crd[0], nd_crd[1]
x0t, y0t = x0 + dx, y0 + dy
ax.arrow(x0, y0,
dx, dy,
lw=3,
head_width=0.1 * sfac, head_length=0.2 * sfac,
fc=kolor, ec=kolor,
length_includes_head=True, shape='full',
joinstyle='round')
ax.text(x0t, y0t,
f' {abs(reaction):.5g}', color=kolor, va='bottom', ha=ha_my)
elif dof == 3:
kolor = 'r'
kier2 = np.sign(reaction)
reaction_str = f'{abs(reaction):.5g}'
if kier2 > 0:
marker_type = r'$\curvearrowleft$'
ax.text(nd_crd[0], nd_crd[1], f'\n {reaction_str}',
color=kolor, va='top', ha='right')
elif kier2 < 0:
marker_type = r'$\curvearrowright$'
ax.text(nd_crd[0], nd_crd[1], f'\n {reaction_str}',
color=kolor, va='top', ha='left')
ax.plot(nd_crd[0], nd_crd[1], marker=marker_type, markersize=30,
color=kolor)
ax.axis('equal')
ax.grid(False)
return ax
def plot_reactions_3d(sfac, fig_wi_he, fig_lbrt, fmt_model_loads,
truss_node_offset, local_axes, ax):
"""Display support reactions of the 3d models.
This is a local function: use plot_reactions()
"""
if not ax:
if fig_wi_he:
fig_wi, fig_he = fig_wi_he
fig = plt.figure(figsize=(fig_wi / 2.54, fig_he / 2.54))
else:
fig = plt.figure()
if fig_lbrt:
fleft, fbottom, fright, ftop = fig_lbrt
fig.subplots_adjust(left=fleft, bottom=fbottom, right=fright, top=ftop)
ax = fig.add_subplot(111, projection=Axes3D.name)
ax = plot_model(node_labels=0, element_labels=0, fmt_model=fmt_model_loads,
node_supports=False, local_axes=local_axes, ax=ax)
if not sfac:
ratio = 0.15 # initial scale 10% of the max dimension
min_x, max_x = ax.get_xlim()
min_y, max_y = ax.get_ylim()
min_z, max_z = ax.get_zlim()
xsfac = ratio * abs(max_x - min_x)
ysfac = ratio * abs(max_y - min_y)
zsfac = ratio * abs(max_z - min_z)
sfac = max(xsfac, ysfac, zsfac)
fixed_nodes = ops.getFixedNodes()
for node_tag in fixed_nodes:
nd_crd = ops.nodeCoord(node_tag)
fixed_dofs = ops.getFixedDOFs(node_tag)
for dof in fixed_dofs:
reaction = ops.nodeReaction(node_tag, dof)
dx, dy, dz = 0., 0., 0.
if not isclose(reaction, 0., abs_tol=1e-9):
if dof in [1, 2, 3]:
kolor = 'b'
if dof == 1:
dx = sfac * np.sign(reaction)
# xcomp = 1
elif dof == 2:
dy = sfac * np.sign(reaction)
# ycomp = 1
elif dof == 3:
dz = sfac * np.sign(reaction)
# zcomp = 1
nn = 1.0 # additional scaling for arrows
ax.arrow3D(*nd_crd, nn * dx, nn * dy, nn * dz,
mutation_scale=10, ec=kolor, fc=kolor)
ax.text(nd_crd[0] + dx,
nd_crd[1] + dy,
nd_crd[2] + dz,
f'{abs(reaction):.5g}', va='bottom', ha='center', color=kolor)
if dof in [4, 5, 6]:
kolor = 'g'
if dof == 4:
dx = sfac * np.sign(reaction)
elif dof == 5:
dy = sfac * np.sign(reaction)
elif dof == 6:
dz = sfac * np.sign(reaction)
nn = 0.5 # additional scaling for arrows
ax.arrow3D(nd_crd[0] + dx,
nd_crd[1] + dy,
nd_crd[2] + dz,
nn * 2 * dx,
nn * 2 * dy,
nn * 2 * dz,
mutation_scale=10, ec=kolor, fc=kolor)
ax.arrow3D(nd_crd[0] + 1.5 * dx,
nd_crd[1] + 1.5 * dy,
nd_crd[2] + 1.5 * dz,
nn * 2 * dx,
nn * 2 * dy,
nn * 2 * dz,
mutation_scale=10, ec=kolor, fc=kolor)
ax.text(nd_crd[0] + 3 * dx,
nd_crd[1] + 3 * dy,
nd_crd[2] + 3 * dz,
f'{abs(reaction):.5g}', va='bottom', ha='center', color=kolor)
return ax
[docs]
def plot_reactions(sfac=False, fig_wi_he=False,
fig_lbrt=False, fmt_model_loads=fmt_model_loads,
truss_node_offset=0, local_axes=False, ax=False):
"""Display support reactions of 2d and 3d models.
"""
ndim = ops.getNDM()[0]
ops.reactions()
if ndim == 2:
ax = plot_reactions_2d(sfac, fig_wi_he,
fig_lbrt, fmt_model_loads,
truss_node_offset, ax)
elif ndim == 3:
ax = plot_reactions_3d(sfac, fig_wi_he,
fig_lbrt, fmt_model_loads,
truss_node_offset,
local_axes, ax)
else:
print(f'\nWarning! ndim: {ndim} not supported yet.')
return ax
[docs]
def plot_extruded_shapes_3d(ele_shapes, az_el=az_el,
fig_wi_he=False,
fig_lbrt=False, ax=False):
"""Plot an extruded 3d model based on cross-section dimenions.
Three arrows present local section axes: green - local x-axis,
red - local z-axis, blue - local y-axis.
Args:
ele_shapes (dict): keys are ele_tags and values are lists of:
shape_type (str): 'rect' - rectangular shape, 'I' - double T shape
and shape_args (list): list of floats, which necessary section
dimensions.
For 'rect' the list is [b d]; width and depth,
for 'I' shape - [bf d tw tf]; flange width, section depth, web
and flange thicknesses
Example: ele_shapes = {1: ['rect', [b, d]],
2: ['I', [bf, d, tw, tf]]}
az_el (tuple): azimuth and elevation
fig_wi_he: figure width and height in centimeters
fig_lbrt: figure left, bottom, right, top boundaries
Usage:
::
ele_shapes = {1: ['circ', [b]],
2: ['rect', [b, h]],
3: ['I', [b, h, b/10., h/6.]]}
opsv.plot_extruded_shapes_3d(ele_shapes)
Notes:
- For now only rectangular, circular and double T sections are supported.
- This function can be a source of inconsistency because OpenSees lacks
functions to return section dimensions. A workaround is to have own
Python helper functions to reuse data specified once
"""
ele_tags = ops.getEleTags()
azim, elev = az_el
if not ax:
if fig_wi_he:
fig_wi, fig_he = fig_wi_he
fig = plt.figure(figsize=(fig_wi / 2.54, fig_he / 2.54))
# fig.subplots_adjust(left=.08, bottom=.08, right=.985, top=.94)
else:
fig = plt.figure()
if fig_lbrt:
fleft, fbottom, fright, ftop = fig_lbrt
fig.subplots_adjust(left=fleft, bottom=fbottom, right=fright, top=ftop)
ax = fig.add_subplot(111, projection=Axes3D.name)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.view_init(azim=azim, elev=elev)
for i, ele_tag in enumerate(ele_tags):
ele_node_tags = ops.eleNodes(ele_tag)
nen = 2
ecrd = np.zeros((nen, 3))
for i, ele_node_tag in enumerate(ele_node_tags):
ecrd[i, :] = ops.nodeCoord(ele_node_tag)
# mesh outline
ax.plot(ecrd[:, 0], ecrd[:, 1], ecrd[:, 2], 'k--', solid_capstyle='round',
solid_joinstyle='round', dash_capstyle='butt',
dash_joinstyle='round')
# eo = Eo[i, :]
xloc = ops.eleResponse(ele_tag, 'xlocal')
yloc = ops.eleResponse(ele_tag, 'ylocal')
zloc = ops.eleResponse(ele_tag, 'zlocal')
g = np.vstack((xloc, yloc, zloc))
G, L = rot_transf_3d(ecrd, g)
# by default empty
shape_type, shape_args = ['circ', [0.]]
if ele_tag in ele_shapes:
shape_type, shape_args = ele_shapes[ele_tag]
if shape_type == 'rect' or shape_type == 'I':
if shape_type == 'rect':
verts = _plot_extruded_shapes_3d_rect(ecrd, g,
shape_args)
elif shape_type == 'I':
verts = _plot_extruded_shapes_3d_double_T(ecrd, g,
shape_args)
# plot 3d sides
ax.add_collection3d(Poly3DCollection(verts, linewidths=0.6,
edgecolors='k', alpha=.25))
elif shape_type == 'circ':
X, Y, Z = _plot_extruded_shapes_3d_circ(ecrd, g,
shape_args)
ax.plot_surface(X, Y, Z, linewidths=0.4, edgecolors='k',
alpha=.25)
# common for all members
Xm, Ym, Zm = sum(ecrd[:, 0]) / 2, sum(ecrd[:, 1]) / 2, sum(ecrd[:, 2]) / 2
alen = 0.1*L
# plot local axis directional vectors: workaround quiver = arrow
ax.quiver(Xm, Ym, Zm, g[0, 0], g[0, 1], g[0, 2], color='g',
lw=2, length=alen, alpha=.8, normalize=True)
ax.quiver(Xm, Ym, Zm, g[1, 0], g[1, 1], g[1, 2], color='b',
lw=2, length=alen, alpha=.8, normalize=True)
ax.quiver(Xm, Ym, Zm, g[2, 0], g[2, 1], g[2, 2], color='r',
lw=2, length=alen, alpha=.8, normalize=True)
ax.set_box_aspect((np.ptp(ax.get_xlim3d()),
np.ptp(ax.get_ylim3d()),
np.ptp(ax.get_zlim3d())))
def _plot_extruded_shapes_3d_double_T(ecrd, g, shape_args):
bf, d, tw, tf = shape_args
za, zb = tw / 2, bf / 2
ya, yb = d / 2 - tf, d / 2
g10a, g11a, g12a = g[1, 0] * ya, g[1, 1] * ya, g[1, 2] * ya
g10b, g11b, g12b = g[1, 0] * yb, g[1, 1] * yb, g[1, 2] * yb
g20a, g21a, g22a = g[2, 0] * za, g[2, 1] * za, g[2, 2] * za
g20b, g21b, g22b = g[2, 0] * zb, g[2, 1] * zb, g[2, 2] * zb
pts = np.zeros((24, 3))
# beg node (I) cross-section vertices crds, counter-clockwise order
pts[0] = [ecrd[0, 0] - g10b - g20b, ecrd[0, 1] - g11b - g21b, ecrd[0, 2] - g12b - g22b]
pts[1] = [ecrd[0, 0] - g10a - g20b, ecrd[0, 1] - g11a - g21b, ecrd[0, 2] - g12a - g22b]
pts[2] = [ecrd[0, 0] - g10a - g20a, ecrd[0, 1] - g11a - g21a, ecrd[0, 2] - g12a - g22a]
pts[3] = [ecrd[0, 0] + g10a - g20a, ecrd[0, 1] + g11a - g21a, ecrd[0, 2] + g12a - g22a]
pts[4] = [ecrd[0, 0] + g10a - g20b, ecrd[0, 1] + g11a - g21b, ecrd[0, 2] + g12a - g22b]
pts[5] = [ecrd[0, 0] + g10b - g20b, ecrd[0, 1] + g11b - g21b, ecrd[0, 2] + g12b - g22b]
pts[6] = [ecrd[0, 0] + g10b + g20b, ecrd[0, 1] + g11b + g21b, ecrd[0, 2] + g12b + g22b]
pts[7] = [ecrd[0, 0] + g10a + g20b, ecrd[0, 1] + g11a + g21b, ecrd[0, 2] + g12a + g22b]
pts[8] = [ecrd[0, 0] + g10a + g20a, ecrd[0, 1] + g11a + g21a, ecrd[0, 2] + g12a + g22a]
pts[9] = [ecrd[0, 0] - g10a + g20a, ecrd[0, 1] - g11a + g21a, ecrd[0, 2] - g12a + g22a]
pts[10] = [ecrd[0, 0] - g10a + g20b, ecrd[0, 1] - g11a + g21b, ecrd[0, 2] - g12a + g22b]
pts[11] = [ecrd[0, 0] - g10b + g20b, ecrd[0, 1] - g11b + g21b, ecrd[0, 2] - g12b + g22b]
# end node (J) cross-section vertices
pts[12] = [ecrd[1, 0] - g10b - g20b, ecrd[1, 1] - g11b - g21b, ecrd[1, 2] - g12b - g22b]
pts[13] = [ecrd[1, 0] - g10a - g20b, ecrd[1, 1] - g11a - g21b, ecrd[1, 2] - g12a - g22b]
pts[14] = [ecrd[1, 0] - g10a - g20a, ecrd[1, 1] - g11a - g21a, ecrd[1, 2] - g12a - g22a]
pts[15] = [ecrd[1, 0] + g10a - g20a, ecrd[1, 1] + g11a - g21a, ecrd[1, 2] + g12a - g22a]
pts[16] = [ecrd[1, 0] + g10a - g20b, ecrd[1, 1] + g11a - g21b, ecrd[1, 2] + g12a - g22b]
pts[17] = [ecrd[1, 0] + g10b - g20b, ecrd[1, 1] + g11b - g21b, ecrd[1, 2] + g12b - g22b]
pts[18] = [ecrd[1, 0] + g10b + g20b, ecrd[1, 1] + g11b + g21b, ecrd[1, 2] + g12b + g22b]
pts[19] = [ecrd[1, 0] + g10a + g20b, ecrd[1, 1] + g11a + g21b, ecrd[1, 2] + g12a + g22b]
pts[20] = [ecrd[1, 0] + g10a + g20a, ecrd[1, 1] + g11a + g21a, ecrd[1, 2] + g12a + g22a]
pts[21] = [ecrd[1, 0] - g10a + g20a, ecrd[1, 1] - g11a + g21a, ecrd[1, 2] - g12a + g22a]
pts[22] = [ecrd[1, 0] - g10a + g20b, ecrd[1, 1] - g11a + g21b, ecrd[1, 2] - g12a + g22b]
pts[23] = [ecrd[1, 0] - g10b + g20b, ecrd[1, 1] - g11b + g21b, ecrd[1, 2] - g12b + g22b]
# list of sides
verts = [[pts[0], pts[1], pts[2], pts[3], pts[4], pts[5], pts[6],
pts[7], pts[8], pts[9], pts[10], pts[11]], # beg
[pts[12], pts[13], pts[14], pts[15], pts[16], pts[17],
pts[18], pts[19], pts[20], pts[21], pts[22], pts[23]], # end
[pts[0], pts[12], pts[23], pts[11]], # bot 1
[pts[5], pts[17], pts[18], pts[6]], # top 2
[pts[9], pts[8], pts[20], pts[21]], # 3
[pts[2], pts[3], pts[15], pts[14]], # 4
[pts[11], pts[10], pts[22], pts[23]], # 5
[pts[0], pts[1], pts[13], pts[12]], # 6
[pts[7], pts[6], pts[18], pts[19]], # 7
[pts[4], pts[5], pts[17], pts[16]], # 8
[pts[10], pts[9], pts[21], pts[22]], # 9
[pts[2], pts[1], pts[13], pts[14]], # 10
[pts[7], pts[8], pts[20], pts[19]], # 11
[pts[3], pts[4], pts[16], pts[15]]] # 12
return verts
def _plot_extruded_shapes_3d_rect(ecrd, g, shape_args):
b, d = shape_args
b2, d2 = b / 2, d / 2
g10, g11, g12 = g[1, 0] * d2, g[1, 1] * d2, g[1, 2] * d2
g20, g21, g22 = g[2, 0] * b2, g[2, 1] * b2, g[2, 2] * b2
pts = np.zeros((8, 3))
# beg node cross-section vertices
# collected i-beg, j-end node coordinates, counter-clockwise order
pts[0] = [ecrd[0, 0] - g10 - g20, ecrd[0, 1] - g11 - g21, ecrd[0, 2] - g12 - g22]
pts[1] = [ecrd[0, 0] + g10 - g20, ecrd[0, 1] + g11 - g21, ecrd[0, 2] + g12 - g22]
pts[2] = [ecrd[0, 0] + g10 + g20, ecrd[0, 1] + g11 + g21, ecrd[0, 2] + g12 + g22]
pts[3] = [ecrd[0, 0] - g10 + g20, ecrd[0, 1] - g11 + g21, ecrd[0, 2] - g12 + g22]
# end node cross-section vertices
pts[4] = [ecrd[1, 0] - g10 - g20, ecrd[1, 1] - g11 - g21, ecrd[1, 2] - g12 - g22]
pts[5] = [ecrd[1, 0] + g10 - g20, ecrd[1, 1] + g11 - g21, ecrd[1, 2] + g12 - g22]
pts[6] = [ecrd[1, 0] + g10 + g20, ecrd[1, 1] + g11 + g21, ecrd[1, 2] + g12 + g22]
pts[7] = [ecrd[1, 0] - g10 + g20, ecrd[1, 1] - g11 + g21, ecrd[1, 2] - g12 + g22]
# list of 4-node sides
verts = [[pts[0], pts[1], pts[2], pts[3]], # beg
[pts[4], pts[5], pts[6], pts[7]], # end
[pts[0], pts[4], pts[5], pts[1]], # bottom
[pts[3], pts[7], pts[6], pts[2]], # top
[pts[0], pts[4], pts[7], pts[3]], # front
[pts[1], pts[5], pts[6], pts[2]]] # back
return verts
def _plot_extruded_shapes_3d_circ(ecrd, g, shape_args):
d = shape_args[0]
r = d / 2
Lxyz = ecrd[1, :] - ecrd[0, :]
L = np.sqrt(Lxyz @ Lxyz)
xl = np.linspace(0, L, 3) # subdivde member length
alpha = np.linspace(0, 2 * np.pi, 10) # subdivide circle
xl, alpha = np.meshgrid(xl, alpha)
X = (ecrd[0, 0] + g[0, 0] * xl + r * np.sin(alpha) * g[1, 0]
+ r * np.cos(alpha) * g[2, 0])
Y = (ecrd[0, 1] + g[0, 1] * xl + r * np.sin(alpha) * g[1, 1]
+ r * np.cos(alpha) * g[2, 1])
Z = (ecrd[0, 2] + g[0, 2] * xl + r * np.sin(alpha) * g[1, 2]
+ r * np.cos(alpha) * g[2, 2])
return X, Y, Z
def bar_length(ecrd):
Lxyz = ecrd[1, :] - ecrd[0, :]
L = np.sqrt(Lxyz @ Lxyz)
return L
def rot_transf_2d(ecrd):
Lxy = ecrd[1, :] - ecrd[0, :]
L = np.sqrt(Lxy @ Lxy)
cosa, cosb = Lxy / L
G = np.array([[cosa, cosb, 0., 0., 0., 0.],
[-cosb, cosa, 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0.],
[0., 0., 0., cosa, cosb, 0.],
[0., 0., 0., -cosb, cosa, 0.],
[0., 0., 0., 0., 0., 1.]])
return G, L, cosa, cosb
def rot_transf_3d(ecrd, g):
Lxyz = ecrd[1, :] - ecrd[0, :]
L = np.sqrt(Lxyz @ Lxyz)
z = np.zeros((3, 3))
G = np.block([[g, z, z, z],
[z, g, z, z],
[z, z, g, z],
[z, z, z, g]])
return G, L