import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from .settings import *
def stress_2d_ele_tags_only(ele_tags):
Stress2dEleClasstags = set([EleClassTag.tri3n,
EleClassTag.tri6n,
EleClassTag.quad4n,
EleClassTag.quad8n,
EleClassTag.quad9n])
ele_classtags = ops.getEleClassTags()
idx = [i for i, e in enumerate(ele_classtags) if e in Stress2dEleClasstags]
ele_tags_2d_tris_quads_only = [ele_tags[i] for i in idx]
return ele_tags_2d_tris_quads_only
[docs]def sig_out_per_node(how_many='all'):
"""Return a 2d numpy array of stress components per OpenSees node.
Three first stress components (sxx, syy, sxy) are calculated and
extracted from OpenSees, while the rest svm (Huber-Mieses-Hencky),
two principal stresses (s1, s2) and directional angle are calculated
as postprocessed quantities.
Args:
how_many (str): supported options are: 'all' - all components,
'sxx', 'syy', 'sxy', 'svm' (or 'vmis'), 's1', 's2', 'angle'.
Returns:
sig_out (ndarray): a 2d array of stress components per node with
the following components: sxx, syy, sxy, svm, s1, s2, angle.
Size (n_nodes x 7).
Examples:
sig_out = opsv.sig_out_per_node()
Notes:
s1, s2: principal stresses
angle: angle of the principal stress s1
"""
ele_tags = ops.getEleTags()
node_tags = ops.getNodeTags()
n_nodes = len(node_tags)
ele_classtag = ops.getEleClassTags(ele_tags[0])[0]
# initialize helper arrays
sig_out = np.zeros((n_nodes, 7))
nodes_tag_count = np.zeros((n_nodes, 2), dtype=int)
nodes_tag_count[:, 0] = node_tags
nen = np.shape(ops.eleNodes(ele_tags[0]))[0]
for i, ele_tag in enumerate(ele_tags):
ele_node_tags = ops.eleNodes(ele_tag)
tmp_list = [0]*nen
for j, ele_node_tag in enumerate(ele_node_tags):
tmp_list[j] = node_tags.index(ele_node_tag)
nodes_tag_count[tmp_list, 1] += 1
if ele_classtag == EleClassTag.SSPquad:
sig_ip_el = ops.eleResponse(ele_tag, 'stress')
# sigM_nd = sigM_ip
sigM_np = np.tile(sig_ip_el, (nen, 1))
else:
sig_nd_el = ops.eleResponse(ele_tag, 'stressAtNodes')
sigM_nd = np.reshape(sig_nd_el, (-1, 3))
# sxx yy xy components
sig_out[tmp_list, 0] += sigM_nd[:nen, 0]
sig_out[tmp_list, 1] += sigM_nd[:nen, 1]
sig_out[tmp_list, 2] += sigM_nd[:nen, 2]
indxs, = np.where(nodes_tag_count[:, 1] > 1)
# n_indxs < n_nodes: e.g. 21<25 (bous), 2<6 (2el) etc.
n_indxs = np.shape(indxs)[0]
# divide summed stresses by the number of common nodes
sig_out[indxs, :] = \
sig_out[indxs, :]/nodes_tag_count[indxs, 1].reshape(n_indxs, 1)
if how_many == 'all' or how_many == 'svm' or how_many == 'vmis':
# warning reshape from (pts,ncomp) to (ncomp,pts)
vm_out = vm_stress(np.transpose(sig_out[:, :3]))
sig_out[:, 3] = vm_out
if (how_many == 'all' or how_many == 's1' or how_many == 's2' or how_many == 'angle'): # noqa: E501
princ_sig_out = princ_stress(np.transpose(sig_out[:, :3]))
sig_out[:, 4:7] = np.transpose(princ_sig_out)
print('-- WARNING!!! full sig_out matrix calculated here --')
return sig_out
def sig_component_per_node(stress_str):
"""Return a 2d numpy array of stress components per OpenSees node.
Three first stress components (sxx, syy, sxy) are calculated and
extracted from OpenSees, while the rest svm (Huber-Mieses-Hencky),
two principal stresses (s1, s2) and directional angle are calculated
as postprocessed quantities.
Args:
how_many (str): supported options are: 'all' - all components,
'sxx', 'syy', 'sxy', 'svm' (or 'vmis'), 's1', 's2', 'angle'.
Returns:
sig_out (ndarray): a 2d array of stress components per node with
the following components: sxx, syy, sxy, svm, s1, s2, angle.
Size (n_nodes x 7).
Examples:
sig_out = opsv.sig_out_per_node()
Notes:
s1, s2: principal stresses
angle: angle of the principal stress s1
"""
ele_tags_all = ops.getEleTags()
ele_tags = stress_2d_ele_tags_only(ele_tags_all)
ele_classtag = ops.getEleClassTags(ele_tags[0])[0]
node_tags = ops.getNodeTags()
n_nodes = len(node_tags)
# initialize helper arrays
sig_out = np.zeros((n_nodes, 4))
nodes_tag_count = np.zeros((n_nodes, 2), dtype=int)
nodes_tag_count[:, 0] = node_tags
nen = np.shape(ops.eleNodes(ele_tags[0]))[0]
for i, ele_tag in enumerate(ele_tags):
ele_node_tags = ops.eleNodes(ele_tag)
tmp_list = [0]*nen
for j, ele_node_tag in enumerate(ele_node_tags):
tmp_list[j] = node_tags.index(ele_node_tag)
nodes_tag_count[tmp_list, 1] += 1
if ele_classtag == EleClassTag.SSPquad:
sig_ip_el = ops.eleResponse(ele_tag, 'stress')
# sigM_nd = sigM_ip
sigM_np = np.tile(sig_ip_el, (nen, 1))
else:
sig_nd_el = ops.eleResponse(ele_tag, 'stressAtNodes')
sigM_nd = np.reshape(sig_nd_el, (-1, 3))
# sxx yy xy components
sig_out[tmp_list, 0] += sigM_nd[:nen, 0]
sig_out[tmp_list, 1] += sigM_nd[:nen, 1]
sig_out[tmp_list, 2] += sigM_nd[:nen, 2]
indxs, = np.where(nodes_tag_count[:, 1] > 1)
# n_indxs < n_nodes: e.g. 21<25 (bous), 2<6 (2el) etc.
n_indxs = np.shape(indxs)[0]
# divide summed stresses by the number of common nodes
sig_out[indxs, :] = \
sig_out[indxs, :]/nodes_tag_count[indxs, 1].reshape(n_indxs, 1)
if stress_str == 'sxx':
sig_out_vec = sig_out[:, 0]
elif stress_str == 'syy':
sig_out_vec = sig_out[:, 1]
elif stress_str == 'sxy':
sig_out_vec = sig_out[:, 2]
elif stress_str == 'svm' or stress_str == 'vmis':
# warning reshape from (pts,ncomp) to (ncomp,pts)
sig_out_vec = vm_stress(np.transpose(sig_out[:, :3]))
elif (stress_str == 's1' or stress_str == 's2' or stress_str == 'angle'):
princ_sig_out = princ_stress(np.transpose(sig_out[:, :3]))
if stress_str == 's1':
# sig_out_vec = np.transpose(princ_sig_out)[:, 0]
sig_out_vec = princ_sig_out[0, :]
elif stress_str == 's2':
sig_out_vec = princ_sig_out[1, :]
elif stress_str == 'angle':
sig_out_vec = princ_sig_out[2, :]
return sig_out_vec
def princ_stress(sig):
"""Return a tuple (s1, s2, angle): principal stresses (plane stress) and angle
Args:
sig (ndarray): input array of stresses at nodes: sxx, syy, sxy (tau)
Returns:
out (ndarray): 1st row is first principal stress s1, 2nd row is second
principal stress s2, 3rd row is the angle of s1
"""
sx, sy, tau = sig[0], sig[1], sig[2]
ds = (sx-sy)/2
R = np.sqrt(ds**2 + tau**2)
s1 = (sx+sy)/2. + R
s2 = (sx+sy)/2. - R
angle = np.arctan2(tau, ds)/2
out = np.vstack((s1, s2, angle))
return out
def vm_stress(sig):
n_sig_comp, n_pts = np.shape(sig)
if n_sig_comp > 3:
x, y, z, xy, xz, yz = sig
else:
x, y, xy = sig
z, xz, yz = 0., 0., 0.
_a = 0.5*((x-y)**2 + (y-z)**2 + (z-x)**2 + 6.*(xy**2 + xz**2 + yz**2))
return np.sqrt(_a)
def quads_to_4tris(quads_conn, nds_crd, nds_val):
"""
Get triangles connectivity, coordinates and new values at quad centroids.
Args:
quads_conn (ndarray):
nds_crd (ndarray):
nds_val (ndarray):
Returns:
tris_conn, nds_c_crd, nds_c_val (tuple):
Notes:
Triangles connectivity array is based on
quadrilaterals connectivity.
Each quad is split into four triangles.
New nodes are created at the quad centroid.
See also:
function: quads_to_8tris_9n, quads_to_8tris_8n
"""
n_quads, _ = quads_conn.shape
n_nds, _ = nds_crd.shape
# coordinates and values at quad centroids _c_
nds_c_crd = np.zeros((n_quads, 2))
nds_c_val = np.zeros(n_quads)
tris_conn = np.zeros((4*n_quads, 3), dtype=int)
for i, quad_conn in enumerate(quads_conn):
j = 4*i
n0, n1, n2, n3 = quad_conn
# quad centroids
nds_c_crd[i] = np.array([np.sum(nds_crd[[n0, n1, n2, n3], 0])/4.,
np.sum(nds_crd[[n0, n1, n2, n3], 1])/4.])
nds_c_val[i] = np.sum(nds_val[[n0, n1, n2, n3]])/4.
# triangles connectivity
tris_conn[j] = np.array([n0, n1, n_nds+i])
tris_conn[j+1] = np.array([n1, n2, n_nds+i])
tris_conn[j+2] = np.array([n2, n3, n_nds+i])
tris_conn[j+3] = np.array([n3, n0, n_nds+i])
return tris_conn, nds_c_crd, nds_c_val
def bricks_to_24tris(bricks_conn, nds_crd, nds_val, disps=None):
"""
Get triangles connectivity, coordinates and new vals at brick face centroids.
Args:
bricks_conn (ndarray):
nds_crd (ndarray):
nds_val (ndarray):
Returns:
tris_conn, nds_c_crd, nds_c_val (tuple):
Notes:
Triangles connectivity array is based on
stdBricks connectivity.
Each of 6 brick faces is split into four triangles.
New nodes are created at the face centroid.
See also:
function: bricks_to_8tris_9n, bricks_to_8tris_8n
"""
n_bricks, _ = bricks_conn.shape
n_nds, _ = nds_crd.shape
# coordinates and values at brick centroids _c_
nds_c_crd = np.zeros((n_bricks*6, 3))
nds_c_val = np.zeros(n_bricks*6)
disps_c = None
if disps is not None:
disps_c = np.zeros((n_bricks*6, 3))
tris_conn = np.zeros((24*n_bricks, 3), dtype=int)
for i, brick_conn in enumerate(bricks_conn):
j = 24*i
n0, n1, n2, n3, n4, n5, n6, n7 = brick_conn
# brick centroids
nds_c_crd[i*6] = np.array([np.sum(nds_crd[[n0, n1, n5, n4], 0])/4.,
np.sum(nds_crd[[n0, n1, n5, n4], 1])/4.,
np.sum(nds_crd[[n0, n1, n5, n5], 2])/4.])
nds_c_crd[i*6+1] = np.array([np.sum(nds_crd[[n1, n2, n6, n5], 0])/4.,
np.sum(nds_crd[[n1, n2, n6, n5], 1])/4.,
np.sum(nds_crd[[n1, n2, n6, n5], 2])/4.])
nds_c_crd[i*6+2] = np.array([np.sum(nds_crd[[n2, n3, n7, n6], 0])/4.,
np.sum(nds_crd[[n2, n3, n7, n6], 1])/4.,
np.sum(nds_crd[[n2, n3, n7, n6], 2])/4.])
nds_c_crd[i*6+3] = np.array([np.sum(nds_crd[[n3, n0, n4, n7], 0])/4.,
np.sum(nds_crd[[n3, n0, n4, n7], 1])/4.,
np.sum(nds_crd[[n3, n0, n4, n7], 2])/4.])
nds_c_crd[i*6+4] = np.array([np.sum(nds_crd[[n4, n5, n6, n7], 0])/4.,
np.sum(nds_crd[[n4, n5, n6, n7], 1])/4.,
np.sum(nds_crd[[n4, n5, n6, n7], 2])/4.])
nds_c_crd[i*6+5] = np.array([np.sum(nds_crd[[n0, n1, n2, n3], 0])/4.,
np.sum(nds_crd[[n0, n1, n2, n3], 1])/4.,
np.sum(nds_crd[[n0, n1, n2, n3], 2])/4.])
nds_c_val[6*i] = np.sum(nds_val[[n0, n1, n5, n4]])/4.
nds_c_val[6*i+1] = np.sum(nds_val[[n1, n2, n6, n5]])/4.
nds_c_val[6*i+2] = np.sum(nds_val[[n2, n3, n7, n6]])/4.
nds_c_val[6*i+3] = np.sum(nds_val[[n3, n0, n4, n7]])/4.
nds_c_val[6*i+4] = np.sum(nds_val[[n4, n5, n6, n7]])/4.
nds_c_val[6*i+5] = np.sum(nds_val[[n0, n1, n2, n3]])/4.
# triangles connectivity
tris_conn[j] = np.array([n0, n1, n_nds+i*6])
tris_conn[j+1] = np.array([n1, n5, n_nds+i*6])
tris_conn[j+2] = np.array([n5, n4, n_nds+i*6])
tris_conn[j+3] = np.array([n4, n0, n_nds+i*6])
tris_conn[j+4] = np.array([n1, n2, n_nds+i*6+1])
tris_conn[j+5] = np.array([n2, n6, n_nds+i*6+1])
tris_conn[j+6] = np.array([n6, n5, n_nds+i*6+1])
tris_conn[j+7] = np.array([n5, n1, n_nds+i*6+1])
tris_conn[j+8] = np.array([n2, n3, n_nds+i*6+2])
tris_conn[j+9] = np.array([n3, n7, n_nds+i*6+2])
tris_conn[j+10] = np.array([n7, n6, n_nds+i*6+2])
tris_conn[j+11] = np.array([n6, n2, n_nds+i*6+2])
tris_conn[j+12] = np.array([n3, n0, n_nds+i*6+3])
tris_conn[j+13] = np.array([n0, n4, n_nds+i*6+3])
tris_conn[j+14] = np.array([n4, n7, n_nds+i*6+3])
tris_conn[j+15] = np.array([n7, n3, n_nds+i*6+3])
tris_conn[j+16] = np.array([n4, n5, n_nds+i*6+4])
tris_conn[j+17] = np.array([n5, n6, n_nds+i*6+4])
tris_conn[j+18] = np.array([n6, n7, n_nds+i*6+4])
tris_conn[j+19] = np.array([n7, n4, n_nds+i*6+4])
tris_conn[j+20] = np.array([n0, n1, n_nds+i*6+5])
tris_conn[j+21] = np.array([n1, n2, n_nds+i*6+5])
tris_conn[j+22] = np.array([n2, n3, n_nds+i*6+5])
tris_conn[j+23] = np.array([n3, n0, n_nds+i*6+5])
if disps is not None:
disps_c[6*i] = np.sum(disps[[n0, n1, n5, n4]], axis=0)/4.
disps_c[6*i+1] = np.sum(disps[[n1, n2, n6, n5]], axis=0)/4.
disps_c[6*i+2] = np.sum(disps[[n2, n3, n7, n6]], axis=0)/4.
disps_c[6*i+3] = np.sum(disps[[n3, n0, n4, n7]], axis=0)/4.
disps_c[6*i+4] = np.sum(disps[[n4, n5, n6, n7]], axis=0)/4.
disps_c[6*i+5] = np.sum(disps[[n0, n1, n2, n3]], axis=0)/4.
return tris_conn, nds_c_crd, nds_c_val, disps_c
# brick20n bricks to tris
def bricks_to_48tris(bricks_conn, nds_crd, nds_val, disps=None):
"""
Get triangles connectivity, coordinates and new vals at brick face centroids,
for brick20n
Args:
bricks_conn (ndarray):
nds_crd (ndarray):
nds_val (ndarray):
Returns:
tris_conn, nds_c_crd, nds_c_val (tuple):
Notes:
Triangles connectivity array is based on
stdBricks connectivity.
Each of 6 brick faces is split into four triangles.
New nodes are created at the face centroid.
See also:
function: bricks_to_24tris
"""
n_bricks, _ = bricks_conn.shape
n_nds, _ = nds_crd.shape
# coordinates and values at brick centroids _c_
nds_c_crd = np.zeros((n_bricks*6, 3))
nds_c_val = np.zeros(n_bricks*6)
disps_c = None
if disps is not None:
disps_c = np.zeros((n_bricks*6, 3))
tris_conn = np.zeros((48*n_bricks, 3), dtype=int)
for i, brick_conn in enumerate(bricks_conn):
j = 48*i
n0, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, \
n14, n15, n16, n17, n18, n19 = brick_conn
# brick centroids
nds_c_crd[i*6] = np.array([np.sum(nds_crd[[n0, n1, n5, n4], 0])/4.,
np.sum(nds_crd[[n0, n1, n5, n4], 1])/4.,
np.sum(nds_crd[[n0, n1, n5, n5], 2])/4.])
nds_c_crd[i*6+1] = np.array([np.sum(nds_crd[[n1, n2, n6, n5], 0])/4.,
np.sum(nds_crd[[n1, n2, n6, n5], 1])/4.,
np.sum(nds_crd[[n1, n2, n6, n5], 2])/4.])
nds_c_crd[i*6+2] = np.array([np.sum(nds_crd[[n2, n3, n7, n6], 0])/4.,
np.sum(nds_crd[[n2, n3, n7, n6], 1])/4.,
np.sum(nds_crd[[n2, n3, n7, n6], 2])/4.])
nds_c_crd[i*6+3] = np.array([np.sum(nds_crd[[n3, n0, n4, n7], 0])/4.,
np.sum(nds_crd[[n3, n0, n4, n7], 1])/4.,
np.sum(nds_crd[[n3, n0, n4, n7], 2])/4.])
nds_c_crd[i*6+4] = np.array([np.sum(nds_crd[[n4, n5, n6, n7], 0])/4.,
np.sum(nds_crd[[n4, n5, n6, n7], 1])/4.,
np.sum(nds_crd[[n4, n5, n6, n7], 2])/4.])
nds_c_crd[i*6+5] = np.array([np.sum(nds_crd[[n0, n1, n2, n3], 0])/4.,
np.sum(nds_crd[[n0, n1, n2, n3], 1])/4.,
np.sum(nds_crd[[n0, n1, n2, n3], 2])/4.])
nds_c_val[6*i] = np.sum(nds_val[[n0, n1, n5, n4]])/4.
nds_c_val[6*i+1] = np.sum(nds_val[[n1, n2, n6, n5]])/4.
nds_c_val[6*i+2] = np.sum(nds_val[[n2, n3, n7, n6]])/4.
nds_c_val[6*i+3] = np.sum(nds_val[[n3, n0, n4, n7]])/4.
nds_c_val[6*i+4] = np.sum(nds_val[[n4, n5, n6, n7]])/4.
nds_c_val[6*i+5] = np.sum(nds_val[[n0, n1, n2, n3]])/4.
# triangles connectivity
tris_conn[j] = np.array([n0, n8, n_nds+i*6])
tris_conn[j+1] = np.array([n8, n1, n_nds+i*6])
tris_conn[j+2] = np.array([n1, n17, n_nds+i*6])
tris_conn[j+3] = np.array([n17, n5, n_nds+i*6])
tris_conn[j+4] = np.array([n5, n12, n_nds+i*6])
tris_conn[j+5] = np.array([n12, n4, n_nds+i*6])
tris_conn[j+6] = np.array([n4, n16, n_nds+i*6])
tris_conn[j+7] = np.array([n16, n0, n_nds+i*6])
tris_conn[j+8] = np.array([n1, n9, n_nds+i*6+1])
tris_conn[j+9] = np.array([n9, n2, n_nds+i*6+1])
tris_conn[j+10] = np.array([n2, n18, n_nds+i*6+1])
tris_conn[j+11] = np.array([n18, n6, n_nds+i*6+1])
tris_conn[j+12] = np.array([n6, n13, n_nds+i*6+1])
tris_conn[j+13] = np.array([n13, n5, n_nds+i*6+1])
tris_conn[j+14] = np.array([n5, n17, n_nds+i*6+1])
tris_conn[j+15] = np.array([n17, n1, n_nds+i*6+1])
tris_conn[j+16] = np.array([n2, n10, n_nds+i*6+2])
tris_conn[j+17] = np.array([n10, n3, n_nds+i*6+2])
tris_conn[j+18] = np.array([n3, n19, n_nds+i*6+2])
tris_conn[j+19] = np.array([n19, n7, n_nds+i*6+2])
tris_conn[j+20] = np.array([n7, n14, n_nds+i*6+2])
tris_conn[j+21] = np.array([n14, n6, n_nds+i*6+2])
tris_conn[j+22] = np.array([n6, n18, n_nds+i*6+2])
tris_conn[j+23] = np.array([n18, n2, n_nds+i*6+2])
tris_conn[j+24] = np.array([n3, n11, n_nds+i*6+3])
tris_conn[j+25] = np.array([n11, n0, n_nds+i*6+3])
tris_conn[j+26] = np.array([n0, n16, n_nds+i*6+3])
tris_conn[j+27] = np.array([n16, n4, n_nds+i*6+3])
tris_conn[j+28] = np.array([n4, n15, n_nds+i*6+3])
tris_conn[j+29] = np.array([n15, n7, n_nds+i*6+3])
tris_conn[j+30] = np.array([n7, n19, n_nds+i*6+3])
tris_conn[j+31] = np.array([n19, n3, n_nds+i*6+3])
tris_conn[j+32] = np.array([n4, n12, n_nds+i*6+4])
tris_conn[j+33] = np.array([n12, n5, n_nds+i*6+4])
tris_conn[j+34] = np.array([n5, n13, n_nds+i*6+4])
tris_conn[j+35] = np.array([n13, n6, n_nds+i*6+4])
tris_conn[j+36] = np.array([n6, n14, n_nds+i*6+4])
tris_conn[j+37] = np.array([n14, n7, n_nds+i*6+4])
tris_conn[j+38] = np.array([n7, n15, n_nds+i*6+4])
tris_conn[j+39] = np.array([n15, n4, n_nds+i*6+4])
tris_conn[j+40] = np.array([n0, n8, n_nds+i*6+5])
tris_conn[j+41] = np.array([n8, n1, n_nds+i*6+5])
tris_conn[j+42] = np.array([n1, n9, n_nds+i*6+5])
tris_conn[j+43] = np.array([n9, n2, n_nds+i*6+5])
tris_conn[j+44] = np.array([n2, n10, n_nds+i*6+5])
tris_conn[j+45] = np.array([n10, n3, n_nds+i*6+5])
tris_conn[j+46] = np.array([n3, n11, n_nds+i*6+5])
tris_conn[j+47] = np.array([n11, n0, n_nds+i*6+5])
if disps is not None:
disps_c[6*i] = np.sum(disps[[n0, n1, n5, n4,
n8, n17, n12, n16]], axis=0)/4.
disps_c[6*i+1] = np.sum(disps[[n1, n2, n6, n5,
n9, n18, n13, n17]], axis=0)/4.
disps_c[6*i+2] = np.sum(disps[[n2, n3, n7, n6,
n10, n19, n14, n18]], axis=0)/4.
disps_c[6*i+3] = np.sum(disps[[n3, n0, n4, n7,
n11, n16, n15, n19]], axis=0)/4.
disps_c[6*i+4] = np.sum(disps[[n4, n5, n6, n7,
n12, n13, n14, n15]], axis=0)/4.
disps_c[6*i+5] = np.sum(disps[[n0, n1, n2, n3,
n8, n9, n10, n11]], axis=0)/4.
return tris_conn, nds_c_crd, nds_c_val, disps_c
def tetra4n_to_4tris(tetras4n_conn):
"""Get triangles connectivity.
Four-node tetrahedron is subdivided into four triangles
Args:
tetra4n_conn (ndarray):
Returns:
tris_conn_subdiv (ndarray):
"""
n_tetras, _ = tetras4n_conn.shape
# n_nds, _ = nds_crd.shape
tris_conn_subdiv = np.zeros((4*n_tetras, 3), dtype=int)
for i, tetra4n_conn in enumerate(tetras4n_conn):
j = 4*i
n0, n1, n2, n3 = tetra4n_conn
# triangles connectivity
tris_conn_subdiv[j] = np.array([n0, n1, n2])
tris_conn_subdiv[j+1] = np.array([n0, n3, n1])
tris_conn_subdiv[j+2] = np.array([n0, n2, n3])
tris_conn_subdiv[j+3] = np.array([n1, n2, n3])
return tris_conn_subdiv
def tetra10n_to_16tris(tetras10n_conn):
"""Get triangles connectivity.
Six-node triangle is subdivided into four triangles
Args:
tetra10n_conn (ndarray):
Returns:
tris_conn_subdiv (ndarray):
"""
n_tetras, _ = tetras10n_conn.shape
# n_nds, _ = nds_crd.shape
tris_conn_subdiv = np.zeros((16*n_tetras, 3), dtype=int)
for i, tetra10n_conn in enumerate(tetras10n_conn):
j = 16*i
n0, n1, n2, n3, n4, n5, n6, n7, n8, n9 = tetra10n_conn
# triangles connectivity
tris_conn_subdiv[j] = np.array([n0, n4, n6])
tris_conn_subdiv[j+1] = np.array([n1, n5, n4])
tris_conn_subdiv[j+2] = np.array([n4, n5, n6])
tris_conn_subdiv[j+3] = np.array([n2, n6, n5])
tris_conn_subdiv[j+4] = np.array([n0, n7, n4])
tris_conn_subdiv[j+5] = np.array([n4, n7, n8])
tris_conn_subdiv[j+6] = np.array([n1, n4, n8])
tris_conn_subdiv[j+7] = np.array([n3, n8, n7])
tris_conn_subdiv[j+8] = np.array([n0, n6, n7])
tris_conn_subdiv[j+9] = np.array([n6, n9, n7])
tris_conn_subdiv[j+10] = np.array([n2, n9, n6])
tris_conn_subdiv[j+11] = np.array([n3, n7, n9])
tris_conn_subdiv[j+12] = np.array([n1, n5, n8])
tris_conn_subdiv[j+13] = np.array([n5, n9, n8])
tris_conn_subdiv[j+14] = np.array([n2, n9, n5])
tris_conn_subdiv[j+15] = np.array([n3, n8, n9])
return tris_conn_subdiv
def tris6n_to_4tris(tris_conn):
"""Get triangles connectivity.
Six-node triangle is subdivided into four triangles
Args:
tris_conn (ndarray):
Returns:
tris_conn_subdiv (ndarray):
"""
n_tris, _ = tris_conn.shape
# n_nds, _ = nds_crd.shape
tris_conn_subdiv = np.zeros((4*n_tris, 3), dtype=int)
for i, tri_conn in enumerate(tris_conn):
j = 4*i
n0, n1, n2, n3, n4, n5 = tri_conn
# triangles connectivity
tris_conn_subdiv[j] = np.array([n0, n3, n5])
tris_conn_subdiv[j+1] = np.array([n3, n1, n4])
tris_conn_subdiv[j+2] = np.array([n3, n4, n5])
tris_conn_subdiv[j+3] = np.array([n5, n4, n2])
return tris_conn_subdiv
# def plot_mesh_2d(nds_crd, eles_conn, lw=0.4, ec='k', ax=False):
def plot_mesh_2d(nds_crd, eles_conn, lw=0.4, ec='k'):
"""
Plot 2d mesh (quads or triangles) outline.
"""
nen = np.shape(eles_conn)[1]
if nen == 3 or nen == 4:
for ele_conn in eles_conn:
x = nds_crd[ele_conn, 0]
y = nds_crd[ele_conn, 1]
plt.fill(x, y, edgecolor=ec, lw=lw, fill=False)
elif nen == 6:
for ele_conn in eles_conn:
x = nds_crd[[ele_conn[0], ele_conn[3], ele_conn[1], ele_conn[4],
ele_conn[2], ele_conn[5]], 0]
y = nds_crd[[ele_conn[0], ele_conn[3], ele_conn[1], ele_conn[4],
ele_conn[2], ele_conn[5]], 1]
plt.fill(x, y, edgecolor=ec, lw=lw, fill=False)
elif nen == 9:
for ele_conn in eles_conn:
x = nds_crd[[ele_conn[0], ele_conn[4], ele_conn[1], ele_conn[5],
ele_conn[2], ele_conn[6], ele_conn[3], ele_conn[7]],
0]
y = nds_crd[[ele_conn[0], ele_conn[4], ele_conn[1], ele_conn[5],
ele_conn[2], ele_conn[6], ele_conn[3], ele_conn[7]],
1]
plt.fill(x, y, edgecolor=ec, lw=lw, fill=False)
elif nen == 8:
for ele_conn in eles_conn:
x = nds_crd[[ele_conn[0], ele_conn[4], ele_conn[1], ele_conn[5],
ele_conn[2], ele_conn[6], ele_conn[3], ele_conn[7]],
0]
y = nds_crd[[ele_conn[0], ele_conn[4], ele_conn[1], ele_conn[5],
ele_conn[2], ele_conn[6], ele_conn[3], ele_conn[7]],
1]
plt.fill(x, y, edgecolor=ec, lw=lw, fill=False)
# def plot_stress_2d(nds_val, mesh_outline=1, cmap='turbo', levels=50,
# fig_wi_he=False, fig_lbrt=False, ax=False):
[docs]def plot_stress_2d(nds_val, mesh_outline=1, cmap='turbo', levels=50):
"""
Plot stress distribution of a 2d elements of a 2d model.
Args:
nds_val (ndarray): the values of a stress component, which can
be extracted from sig_out array (see sig_out_per_node
function)
mesh_outline (int): 1 - mesh is plotted, 0 - no mesh plotted.
cmap (str): Matplotlib color map (default is 'turbo')
Usage:
See demo_quads_4x4.py example.
"""
node_tags, ele_tags_all = ops.getNodeTags(), ops.getEleTags()
ele_tags = stress_2d_ele_tags_only(ele_tags_all)
n_nodes, n_eles = len(node_tags), len(ele_tags)
ele_classtag = ops.getEleClassTags(ele_tags[0])[0]
if (ele_classtag == EleClassTag.tri3n):
nen = 3
elif (ele_classtag == EleClassTag.tri6n):
nen = 6
elif (ele_classtag == EleClassTag.quad4n):
nen = 4
elif (ele_classtag == EleClassTag.quad8n):
nen = 8
elif (ele_classtag == EleClassTag.quad9n):
nen = 9
# nen = np.shape(ops.eleNodes(ele_tags[0]))[0]
# idiom coordinates as ordered in node_tags
# use node_tags.index(tag) for correspondence
nds_crd = np.zeros((n_nodes, 2))
for i, node_tag in enumerate(node_tags):
nds_crd[i] = ops.nodeCoord(node_tag)
# from utils / sig_out_per_node
# fixme: if this can be simplified
# index (starts from 0) to node_tag correspondence
# (a) data in np.array of integers
# nodes_tag_count = np.zeros((n_nodes, 2), dtype=int)
# nodes_tag_count[:, 0] = node_tags
#
# correspondence indx and node_tag is in node_tags.index
# after testing remove the above
eles_conn = np.zeros((n_eles, nen), dtype=int)
for i, ele_tag in enumerate(ele_tags):
ele_node_tags = ops.eleNodes(ele_tag)
for j, ele_node_tag in enumerate(ele_node_tags):
eles_conn[i, j] = node_tags.index(ele_node_tag)
if (ele_classtag == EleClassTag.tri3n):
tris_conn = eles_conn
nds_crd_all = nds_crd
nds_val_all = nds_val
elif (ele_classtag == EleClassTag.quad4n):
tris_conn, nds_c_crd, nds_c_val = \
quads_to_4tris(eles_conn, nds_crd, nds_val)
nds_crd_all = np.vstack((nds_crd, nds_c_crd))
nds_val_all = np.hstack((nds_val, nds_c_val))
elif (ele_classtag == EleClassTag.tri6n):
tris_conn = tris6n_to_4tris(eles_conn)
nds_crd_all = nds_crd
nds_val_all = nds_val
elif (ele_classtag == EleClassTag.quad8n):
tris_conn, nds_c_crd, nds_c_val = quads_to_8tris_8n(eles_conn,
nds_crd, nds_val)
nds_crd_all = np.vstack((nds_crd, nds_c_crd))
nds_val_all = np.hstack((nds_val, nds_c_val))
elif (ele_classtag == EleClassTag.quad9n):
tris_conn = quads_to_8tris_9n(eles_conn)
nds_crd_all = nds_crd
nds_val_all = nds_val
# 1. plot contour maps
triangulation = tri.Triangulation(nds_crd_all[:, 0],
nds_crd_all[:, 1],
tris_conn)
plt.tricontourf(triangulation, nds_val_all, levels=levels, cmap=cmap)
# 2. plot original mesh (quad) without subdivision into triangles
if mesh_outline:
# plot_mesh_2d(nds_crd, eles_conn, ax=ax)
plot_mesh_2d(nds_crd, eles_conn)
plt.colorbar()
plt.axis('equal')
plt.grid(False)
# see also quads_to_8tris_9n
def quads_to_8tris_8n(quads_conn, nds_crd, nds_val):
"""
Get triangles connectivity, coordinates and new values at quad centroids.
Args:
quads_conn (ndarray):
nds_crd (ndarray):
nds_val (ndarray):
Returns:
tris_conn, nds_c_crd, nds_c_val (tuple):
Notes:
Triangles connectivity array is based on
quadrilaterals connectivity.
Each quad is split into eight triangles.
New nodes are created at the quad centroid.
See also:
function: quads_to_8tris_9n, quads_to_4tris
"""
n_quads, _ = quads_conn.shape
n_nds, _ = nds_crd.shape
# coordinates and values at quad centroids _c_
nds_c_crd = np.zeros((n_quads, 2))
nds_c_val = np.zeros(n_quads)
tris_conn = np.zeros((8*n_quads, 3), dtype=int)
for i, quad_conn in enumerate(quads_conn):
j = 8*i
n0, n1, n2, n3, n4, n5, n6, n7 = quad_conn
# quad centroids
# nds_c_crd[i] = np.array([np.sum(nds_crd[[n0, n1, n2, n3], 0])/4.,
# np.sum(nds_crd[[n0, n1, n2, n3], 1])/4.])
# nds_c_val[i] = np.sum(nds_val[[n0, n1, n2, n3]])/4.
nds_c_crd[i] = quad8n_val_at_center(nds_crd[[n0, n1, n2, n3,
n4, n5, n6, n7]])
nds_c_val[i] = quad8n_val_at_center(nds_val[[n0, n1, n2, n3,
n4, n5, n6, n7]])
# triangles connectivity
tris_conn[j] = np.array([n0, n4, n_nds+i])
tris_conn[j+1] = np.array([n4, n1, n_nds+i])
tris_conn[j+2] = np.array([n1, n5, n_nds+i])
tris_conn[j+3] = np.array([n5, n2, n_nds+i])
tris_conn[j+4] = np.array([n2, n6, n_nds+i])
tris_conn[j+5] = np.array([n6, n3, n_nds+i])
tris_conn[j+6] = np.array([n3, n7, n_nds+i])
tris_conn[j+7] = np.array([n7, n0, n_nds+i])
return tris_conn, nds_c_crd, nds_c_val
# see also quads_to_8tris_8n
def quads_to_8tris_9n(quads_conn):
"""
Get triangles connectivity, coordinates and new values at quad centroids.
Args:
quads_conn (ndarray):
nds_crd (ndarray):
nds_val (ndarray):
Returns:
tris_conn, nds_c_crd, nds_c_val (tuple):
Notes:
Triangles connectivity array is based on
quadrilaterals connectivity.
Each quad is split into eight triangles.
New nodes are created at the quad centroid.
See also:
function: quads_to_8tris_8n, quads_to_4tris
"""
n_quads, _ = quads_conn.shape
# n_nds, _ = nds_crd.shape
tris_conn = np.zeros((8*n_quads, 3), dtype=int)
for i, quad_conn in enumerate(quads_conn):
j = 8*i
n0, n1, n2, n3, n4, n5, n6, n7, n8 = quad_conn
# center nodeds already present in quad 9n
# triangles connectivity
tris_conn[j] = np.array([n0, n4, n8])
tris_conn[j+1] = np.array([n4, n1, n8])
tris_conn[j+2] = np.array([n1, n5, n8])
tris_conn[j+3] = np.array([n5, n2, n8])
tris_conn[j+4] = np.array([n2, n6, n8])
tris_conn[j+5] = np.array([n6, n3, n8])
tris_conn[j+6] = np.array([n3, n7, n8])
tris_conn[j+7] = np.array([n7, n0, n8])
# return tris_conn, nds_c_crd, nds_c_val
return tris_conn
def quad8n_val_at_center(vals):
"""
Calculate values at the center of 8-node quad element.
"""
val_c1 = -np.mean(vals[[0, 1, 2, 3]], axis=0) + 2*np.mean(vals[[4, 5, 6, 7]], axis=0) # noqa: E501
# val_c1 = -np.sum(vals[[0, 1, 2, 3]], axis=0)/4 + np.sum(vals[[4, 5, 6, 7]], axis=0)/2 # noqa: E501
# val_c2 = -np.sum(vals[[0, 1, 2, 3]], axis=0)/4 + np.sum(vals[[4, 5, 6, 7]], axis=0)/2 # noqa: E501
# val_c3 = np.mean(vals[[4, 5, 6, 7]], axis=0)
# val_c4 = np.mean(vals[[0, 1, 2, 3]], axis=0)
# val_c5 = np.mean(vals, axis=0)
# return val_c1, val_c2, val_c3, val_c4, val_c5
return val_c1
def plot_stress(stress_str, mesh_outline=1, cmap='turbo', levels=50):
"""Plot stress distribution of the model.
Args:
stress_str (string): stress component string. Available options are:
'sxx', 'syy', 'sxy', 'vmis', 's1', 's2', 'alpha'
mesh_outline (int): 1 - mesh is plotted, 0 - no mesh plotted.
cmap (str): Matplotlib color map (default is 'turbo')
levels (int): number and positions of the contour lines / regions.
Usage:
::
opsv.plot_stress('vmis')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()
See also:
:ref:`opsvis_sig_out_per_node`
"""
ndim = ops.getNDM()[0]
if ndim == 2:
_plot_stress_2d(stress_str, mesh_outline, cmap, levels)
# if axis_off:
# plt.axis('off')
# not implemented yet
# elif ndim == 3:
# _plot_stress_3d(stress_str, mesh_outline, cmap, levels)
else:
print(f'\nWarning! ndim: {ndim} not implemented yet.')
# plt.show() # call this from main py file for more control
def _plot_stress_2d(stress_str, mesh_outline, cmap, levels):
"""See documentation for plot_stress command"""
# node_tags = ops.getNodeTags()
# ele_tags = ops.getEleTags()
# n_nodes = len(node_tags)
# second version - better - possible different types
# of elements (mix of quad and tri)
# for ele_tag in ele_tags:
# nen = np.shape(ops.eleNodes(ele_tag))[0]
# avoid calculating and storing all stress components
# sig_out = sig_out_per_node(stress_str)
# switcher = {'sxx': 0,
# 'syy': 1,
# 'sxy': 2,
# 'svm': 3,
# 'vmis': 3,
# 's1': 4,
# 's2': 5,
# 'angle': 6}
# nds_val = sig_out[:, switcher[stress_str]]
nds_val = sig_component_per_node(stress_str)
plot_stress_2d(nds_val, mesh_outline, cmap, levels)