Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
Lb = normrow(Qb)
Lc = normrow(Qc)
LQn = normrow(Qn)
Lmu = normrow(mu)
a = arccos((La**2 + Lb**2 - Lc**2) / (2 * La * Lb))
k = 2 * sin(a) / Lc
ex = Qn / tile(LQn, (1, 3)) # temporary simplification
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
cma = cross(Mc, Qa)
cmb = cross(Mc, Qb)
ua = cma / tile(normrow(cma), (1, 3))
ub = cmb / tile(normrow(cmb), (1, 3))
c1 = cross(Qa, ua)
c2 = cross(Qb, ub)
Lc1 = normrow(c1)
Lc2 = normrow(c2)
Ms = sum(Mc**2, 1)[:, newaxis]
Sa = ua * tile(Ms * Lc1 / (La * sum(Mc * c1, 1)[:, newaxis]), (1, 3))
Sb = ub * tile(Ms * Lc2 / (Lb * sum(Mc * c2, 1)[:, newaxis]), (1, 3))
Sa[isnan(Sa)] = 0
Sb[isnan(Sb)] = 0
S[inds, :] += Sa
S[indi, :] -= Sa + Sb
S[indf, :] += Sb
# Add node junction duplication for when elements cross each other
# mu[0, :] = -1.25*x[0, :] + 1.5*x[1, :] - 0.25*x[2, :]
# mu[-1, :] = 0.25*x[-3, :] - 1.5*x[-2, :] + 1.25*x[-1, :]
ex = Qn / tile(LQn, (1, 3))
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
cma = cross(Mc, Qa)
cmb = cross(Mc, Qb)
ua = cma / tile(normrow(cma), (1, 3))
ub = cmb / tile(normrow(cmb), (1, 3))
c1 = cross(Qa, ua)
c2 = cross(Qb, ub)
Lc1 = normrow(c1)
Lc2 = normrow(c2)
Ms = sum(Mc**2, 1)[:, newaxis]
Sa = ua * tile(Ms * Lc1 / (La * sum(Mc * c1, 1)[:, newaxis]), (1, 3))
Sb = ub * tile(Ms * Lc2 / (Lb * sum(Mc * c2, 1)[:, newaxis]), (1, 3))
Sa[isnan(Sa)] = 0
Sb[isnan(Sb)] = 0
S[inds, :] += Sa
S[indi, :] -= Sa + Sb
S[indf, :] += Sb
return S
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
cma = cross(Mc, Qa)
cmb = cross(Mc, Qb)
ua = cma / tile(normrow(cma), (1, 3))
ub = cmb / tile(normrow(cmb), (1, 3))
c1 = cross(Qa, ua)
c2 = cross(Qb, ub)
Lc1 = normrow(c1)
Lc2 = normrow(c2)
Ms = sum(Mc**2, 1)[:, newaxis]
Sa = ua * tile(Ms * Lc1 / (La * sum(Mc * c1, 1)[:, newaxis]), (1, 3))
Sb = ub * tile(Ms * Lc2 / (Lb * sum(Mc * c2, 1)[:, newaxis]), (1, 3))
Sa[isnan(Sa)] = 0
Sb[isnan(Sb)] = 0
S[inds, :] += Sa
S[indi, :] -= Sa + Sb
S[indf, :] += Sb
return S
Ct = C.transpose()
Ci = C[:, free]
Cit = Ci.transpose()
Ct2 = Ct.copy()
Ct2.data **= 2
# --------------------------------------------------------------------------
# if none of the initial lengths are set,
# set the initial lengths to the current lengths
# --------------------------------------------------------------------------
if all(linit == 0):
linit = normrow(C.dot(x))
# --------------------------------------------------------------------------
# initial values
# --------------------------------------------------------------------------
q = ones((num_e, 1), dtype=float)
l = normrow(C.dot(x))
f = q * l
v = zeros((num_v, 3), dtype=float)
r = zeros((num_v, 3), dtype=float)
# --------------------------------------------------------------------------
# helpers
# --------------------------------------------------------------------------
def rk(x0, v0, steps=2):
def a(t, v):
dx = v * t
x[free] = x0[free] + dx[free]
# update residual forces
r[free] = p[free] - D.dot(x)
return cb * r / mass
if steps == 1:
Returns:
array: Updated beam nodal shears.
"""
S *= 0
Xs = X[inds, :]
Xi = X[indi, :]
Xf = X[indf, :]
Qa = Xi - Xs
Qb = Xf - Xi
Qc = Xf - Xs
Qn = cross(Qa, Qb)
mu = 0.5 * (Xf - Xs)
La = normrow(Qa)
Lb = normrow(Qb)
Lc = normrow(Qc)
LQn = normrow(Qn)
Lmu = normrow(mu)
a = arccos((La**2 + Lb**2 - Lc**2) / (2 * La * Lb))
k = 2 * sin(a) / Lc
ex = Qn / tile(LQn, (1, 3)) # temporary simplification
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
cma = cross(Mc, Qa)
cmb = cross(Mc, Qb)
ua = cma / tile(normrow(cma), (1, 3))
ub = cmb / tile(normrow(cmb), (1, 3))
c1 = cross(Qa, ua)
c2 = cross(Qb, ub)
a = arccos((La**2 + Lb**2 - Lc**2) / (2 * La * Lb))
k = 2 * sin(a) / Lc
ex = Qn / tile(LQn, (1, 3)) # temporary simplification
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
cma = cross(Mc, Qa)
cmb = cross(Mc, Qb)
ua = cma / tile(normrow(cma), (1, 3))
ub = cmb / tile(normrow(cmb), (1, 3))
c1 = cross(Qa, ua)
c2 = cross(Qb, ub)
Lc1 = normrow(c1)
Lc2 = normrow(c2)
Ms = sum(Mc**2, 1)[:, newaxis]
Sa = ua * tile(Ms * Lc1 / (La * sum(Mc * c1, 1)[:, newaxis]), (1, 3))
Sb = ub * tile(Ms * Lc2 / (Lb * sum(Mc * c2, 1)[:, newaxis]), (1, 3))
Sa[isnan(Sa)] = 0
Sb[isnan(Sb)] = 0
S[inds, :] += Sa
S[indi, :] -= Sa + Sb
S[indf, :] += Sb
# Add node junction duplication for when elements cross each other
# mu[0, :] = -1.25*x[0, :] + 1.5*x[1, :] - 0.25*x[2, :]
# mu[-1, :] = 0.25*x[-3, :] - 1.5*x[-2, :] + 1.25*x[-1, :]
return S
def _beam_shear(S, X, inds, indi, indf, EIx, EIy):
S *= 0
Xs = X[inds, :]
Xi = X[indi, :]
Xf = X[indf, :]
Qa = Xi - Xs
Qb = Xf - Xi
Qc = Xf - Xs
Qn = cross(Qa, Qb)
mu = 0.5 * (Xf - Xs)
La = normrow(Qa)
Lb = normrow(Qb)
Lc = normrow(Qc)
LQn = normrow(Qn)
Lmu = normrow(mu)
a = arccos((La**2 + Lb**2 - Lc**2) / (2 * La * Lb))
k = 2 * sin(a) / Lc
ex = Qn / tile(LQn, (1, 3))
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
EIy (array): Nodal EIy flexural stiffnesses.
Returns:
array: Updated beam nodal shears.
"""
S *= 0
Xs = X[inds, :]
Xi = X[indi, :]
Xf = X[indf, :]
Qa = Xi - Xs
Qb = Xf - Xi
Qc = Xf - Xs
Qn = cross(Qa, Qb)
mu = 0.5 * (Xf - Xs)
La = normrow(Qa)
Lb = normrow(Qb)
Lc = normrow(Qc)
LQn = normrow(Qn)
Lmu = normrow(mu)
a = arccos((La**2 + Lb**2 - Lc**2) / (2 * La * Lb))
k = 2 * sin(a) / Lc
ex = Qn / tile(LQn, (1, 3)) # temporary simplification
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
cma = cross(Mc, Qa)
cmb = cross(Mc, Qb)
ua = cma / tile(normrow(cma), (1, 3))
ub = cmb / tile(normrow(cmb), (1, 3))
Returns:
array: Updated beam nodal shears.
"""
S *= 0
Xs = X[inds, :]
Xi = X[indi, :]
Xf = X[indf, :]
Qa = Xi - Xs
Qb = Xf - Xi
Qc = Xf - Xs
Qn = cross(Qa, Qb)
mu = 0.5 * (Xf - Xs)
La = normrow(Qa)
Lb = normrow(Qb)
Lc = normrow(Qc)
LQn = normrow(Qn)
Lmu = normrow(mu)
a = arccos((La**2 + Lb**2 - Lc**2) / (2 * La * Lb))
k = 2 * sin(a) / Lc
ex = Qn / tile(LQn, (1, 3)) # temporary simplification
ez = mu / tile(Lmu, (1, 3))
ey = cross(ez, ex)
K = tile(k / LQn, (1, 3)) * Qn
Kx = tile(sum(K * ex, 1)[:, newaxis], (1, 3)) * ex
Ky = tile(sum(K * ey, 1)[:, newaxis], (1, 3)) * ey
Mc = EIx * Kx + EIy * Ky
cma = cross(Mc, Qa)
cmb = cross(Mc, Qb)
ua = cma / tile(normrow(cma), (1, 3))
ub = cmb / tile(normrow(cmb), (1, 3))
c1 = cross(Qa, ua)