"""
Functions implementing the shell10x element.
"""
from __future__ import absolute_import
import numpy as nm
from sfepy.linalg import norm_l2_along_axis as norm
from sfepy.linalg import dot_sequences as ddot
from six.moves import range
[docs]
def create_elastic_tensor(young, poisson, shear_correction=True):
"""
Create the elastic tensor with the applied shear correction (the default)
for the shell10x element.
"""
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
# Static condensation so that \sigma_{33} = 0.
mtx = stiffness_from_youngpoisson(3, young, poisson, plane='stress')
mtx[2, :3] = mtx[:2, 2] = 0.0
if shear_correction:
mtx[4, 4] *= 5.0 / 6.0
mtx[5, 5] *= 5.0 / 6.0
return mtx
[docs]
def create_local_bases(coors):
"""
Create local orthonormal bases in each vertex of quadrilateral cells.
Parameters
----------
coors : array
The coordinates of cell vertices, shape `(n_el, 4, 3)`.
Returns
-------
ebs : array
The local bases, shape `(n_el, 4, 3, 3)`. The basis vectors are rows of
the (..., 3, 3) blocks.
"""
ebs = nm.zeros((coors.shape[0], 4, 3, 3), dtype=nm.float64)
ebs[:, 0, 0, :] = coors[:, 1, :] - coors[:, 0, :]
ebs[:, 0, 1, :] = coors[:, 3, :] - coors[:, 0, :]
ebs[:, 0, 2, :] = nm.cross(ebs[:, 0, 0, :], ebs[:, 0, 1, :])
ebs[:, 1, 0, :] = coors[:, 1, :] - coors[:, 0, :]
ebs[:, 1, 1, :] = coors[:, 2, :] - coors[:, 1, :]
ebs[:, 1, 2, :] = nm.cross(ebs[:, 1, 0, :], ebs[:, 1, 1, :])
ebs[:, 2, 0, :] = coors[:, 2, :] - coors[:, 3, :]
ebs[:, 2, 1, :] = coors[:, 2, :] - coors[:, 1, :]
ebs[:, 2, 2, :] = nm.cross(ebs[:, 2, 0, :], ebs[:, 2, 1, :])
ebs[:, 3, 0, :] = coors[:, 2, :] - coors[:, 3, :]
ebs[:, 3, 1, :] = coors[:, 3, :] - coors[:, 0, :]
ebs[:, 3, 2, :] = nm.cross(ebs[:, 3, 0, :], ebs[:, 3, 1, :])
e2 = nm.array([0, 1, 0], dtype=nm.float64)
for ii in range(4):
ebs[:, ii, 0, :] = nm.cross(e2, ebs[:, ii, 2, :])
ebs[:, ii, 1, :] = nm.cross(ebs[:, ii, 2, :], ebs[:, ii, 0, :])
ebs[:, ii] = ebs[:, ii] / nm.linalg.norm(ebs[:, ii], axis=2)[..., None]
return ebs
[docs]
def create_rotation_ops(ebs):
"""
Create operators associated to rotation DOFs.
Parameters
----------
ebs : array
The local bases, shape `(n_el, 4, 3, 3)`.
Returns
-------
rops : array
The rotation operators, shape `(n_el, 4, 3, 3)`.
"""
rops = nm.empty((ebs.shape[0], 4, 3, 3), dtype=nm.float64)
for ii in range(4):
aux = nm.einsum('cj,ck->cjk', ebs[:, ii, 0], ebs[:, ii, 1])
rops[:, ii] = aux - aux.transpose((0, 2, 1))
return rops
[docs]
def get_mapping_data(ebs, rops, ps, coors_loc, qp_coors, qp_weights,
special_dx3=False):
r"""
Compute reference element mapping data for shell10x elements.
Notes
-----
The code assumes that the quadrature points are w.r.t. (:math:`t` =
thickness of the shell) :math:`[0, 1] \times [0, 1] \times [-t/2, t/2]`
reference cell and the quadrature weights are multiplied by :math:`t`.
"""
n_el = coors_loc.shape[0]
n_qp = qp_weights.shape[0]
bfu = ps.eval_basis(qp_coors[:, :2].copy())
bfgu = ps.eval_basis(qp_coors[:, :2].copy(), diff=True)
nh = ebs[..., -1, :]
# h = tilde dzeta = t * (dzeta - 0.5), d = dzeta in [0, 1]
# dh/dd = t
h = qp_coors[:, -1, None, None]
dxdxi = nm.empty((n_el, n_qp, 3, 3), dtype=nm.float64)
coors_loc_3d = coors_loc[:, None, ...] + (nh[:, None, ...] * h)
dxdxi[..., :2, :] = nm.einsum('qij,cqjk->cqik', bfgu, coors_loc_3d)
# The factor of t from dh/dd is not here because quadrature weights were
# multiplied by t. This seems more accurate than having t here and the
# usual weights.
# dx_col/dxi_row
dxdxi[..., 2:, :] = nm.einsum('qij,cjk->cqik', bfu, nh) # * 2
if special_dx3:
dxdxi[..., 2:, :] = [0.0, 0.0, 1.0]
det = nm.linalg.det(dxdxi)
# *2 to have [-1, 1] reference cell numbers...
# dxi_col/dx_row
dxidx = nm.linalg.inv(dxdxi)
if special_dx3:
return dxidx, det
bfg = nm.zeros((n_el, n_qp, 3, 3, 24), dtype=nm.float64)
bfg[..., 0, :2, :4] = bfgu
bfg[..., 1, :2, 4:8] = bfgu
bfg[..., 2, :2, 8:12] = bfgu
n = None
# bfg[:, 0, :2, 12:16] = (bfgu * rops[..., 0, 0][:, n, n, :]) * h = 0
bfg[..., 0, :2, 16:20] = (bfgu * rops[..., 0, 1][:, n, n, :]) * h
bfg[..., 0, :2, 20:24] = (bfgu * rops[..., 0, 2][:, n, n, :]) * h
# bfg[:, 0, 2:, 12:16] = bfu * rops[..., 0, 0][:, n, n, :] = 0
bfg[..., 0, 2:, 16:20] = bfu * rops[..., 0, 1][:, n, n, :]
bfg[..., 0, 2:, 20:24] = bfu * rops[..., 0, 2][:, n, n, :]
bfg[..., 1, :2, 12:16] = (bfgu * rops[..., 1, 0][:, n, n, :]) * h
# bfg[:, 1, :2, 16:20] = (bfgu * rops[..., 1, 1][:, n, n, :]) * h = 0
bfg[..., 1, :2, 20:24] = (bfgu * rops[..., 1, 2][:, n, n, :]) * h
bfg[..., 1, 2:, 12:16] = bfu * rops[..., 1, 0][:, n, n, :]
# bfg[:, 1, 2:, 16:20] = bfu * rops[..., 1, 1][:, n, n, :] = 0
bfg[..., 1, 2:, 20:24] = bfu * rops[..., 1, 2][:, n, n, :]
bfg[..., 2, :2, 12:16] = (bfgu * rops[..., 2, 0][:, n, n, :]) * h
bfg[..., 2, :2, 16:20] = (bfgu * rops[..., 2, 1][:, n, n, :]) * h
# bfg[:, 2, :2, 20:24] = (bfgu * rops[..., 2, 2][:, n, n, :]) * h = 0
bfg[..., 2, 2:, 12:16] = bfu * rops[..., 2, 0][:, n, n, :]
bfg[..., 2, 2:, 16:20] = bfu * rops[..., 2, 1][:, n, n, :]
# bfg[:, 2, 2:, 20:24] = bfu * rops[..., 2, 2][:, n, n, :] = 0
# ... related to [0, 1] vs. [-1, 1]:
# bfg[:, :, 2, 12:] /= 0.5 when *2 in dxidx above
bfgm = nm.einsum('cqij,cqkjl->cqikl', dxidx, bfg)
return coors_loc_3d, bfu, bfgm, dxidx, det
[docs]
def get_dsg_strain(coors_loc, qp_coors):
r"""
Compute DSG strain components.
Returns
-------
dsg : array
The strain matrix components corresponding to :math:`e_{13}`,
:math:`e_{23}`, shape `(n_el, n_qp, 2, 24)`.
Notes
-----
Involves :math:`w`, :math:`\alpha`, :math:`\beta` DOFs.
"""
dsg = nm.zeros((coors_loc.shape[0], qp_coors.shape[0], 2, 24),
dtype=nm.float64)
x1, x2, x3, x4 = coors_loc.transpose((1, 0, 2))
o = nm.ones((coors_loc.shape[0], 1), dtype=nm.float64)
v13_1 = nm.c_[-o, o, -0.5 * (x2[:, 1] - x1[:, 1]),
-0.5 * (x2[:, 1] - x1[:, 1]),
0.5 * (x2[:, 0] - x1[:, 0]),
0.5 * (x2[:, 0] - x1[:, 0])]
i13_1 = nm.r_[8, 9, 12, 13, 16, 17]
v13_2 = nm.c_[o, -o, -0.5 * (x3[:, 1] - x4[:, 1]),
-0.5 * (x3[:, 1] - x4[:, 1]),
0.5 * (x3[:, 0] - x4[:, 0]),
0.5 * (x3[:, 0] - x4[:, 0])]
i13_2 = nm.r_[10, 11, 14, 15, 18, 19]
v23_1 = nm.c_[-o, o, -0.5 * (x4[:, 1] - x1[:, 1]),
-0.5 * (x4[:, 1] - x1[:, 1]),
0.5 * (x4[:, 0] - x1[:, 0]),
0.5 * (x4[:, 0] - x1[:, 0])]
i23_1 = nm.r_[8, 11, 12, 15, 16, 19]
v23_2 = nm.c_[-o, o, -0.5 * (x3[:, 1] - x2[:, 1]),
-0.5 * (x3[:, 1] - x2[:, 1]),
0.5 * (x3[:, 0] - x2[:, 0]),
0.5 * (x3[:, 0] - x2[:, 0])]
i23_2 = nm.r_[9, 10, 13, 14, 17, 18]
c1, c2 = qp_coors[:, 0], qp_coors[:, 1]
# Assuming qp in [-1, 1]
# dsg[:, 0, i13_1] = 0.25 * (1.0 - c2)[:, None] * v13_1
# dsg[:, 0, i13_2] = 0.25 * (1.0 + c2)[:, None] * v13_2
# dsg[:, 1, i23_1] = 0.25 * (1.0 - c1)[:, None] * v23_1
# dsg[:, 1, i23_2] = 0.25 * (1.0 + c1)[:, None] * v23_2
# Assuming qp in [0, 1]
n = None
dsg[..., 0, i13_1] = (1.0 - c2)[n, :, n] * v13_1[:, n, :]
dsg[..., 0, i13_2] = (c2)[n, :, n] * v13_2[:, n, :]
dsg[..., 1, i23_1] = (1.0 - c1)[n, :, n] * v23_1[:, n, :]
dsg[..., 1, i23_2] = (c1)[n, :, n] * v23_2[:, n, :]
return dsg
[docs]
def create_strain_matrix(bfgm, dxidx, dsg):
"""
Create the strain operator matrix.
"""
tg = create_strain_transform(dxidx)
mtx_b = nm.concatenate((bfgm[..., 0:1, 0, :],
bfgm[..., 1:2, 1, :],
bfgm[..., 2:3, 2, :],
bfgm[..., 0:1, 1, :] + bfgm[..., 1:2, 0, :],
nm.einsum('ciaj,cijk->ciak', tg[..., 4:5, 4:], dsg),
nm.einsum('ciaj,cijk->ciak', tg[..., 5:6, 4:], dsg)),
2)
return mtx_b
[docs]
def add_eas_dofs(mtx_b, qp_coors, det, det0, dxidx0):
"""
Add additional strain components [Andelfinger and Ramm] (7 parameters to be
condensed out).
"""
tg0 = create_strain_transform(dxidx0)
# mtx_b is the same as Bk -> mtx_ee must be the same as Ba
# -> qp coors transformation below is needed.
qp_coors = 8 * (qp_coors - 0.5)
c1, c2 = qp_coors[:, 0], qp_coors[:, 1]
cc = c1 * c2
mtx_e = nm.zeros((qp_coors.shape[0], 6, 7), dtype=nm.float64)
mtx_e[:, 0, 0] = mtx_e[:, 3, 2] = c1
mtx_e[:, 1, 1] = mtx_e[:, 3, 3] = c2
mtx_e[:, 0, 4] = mtx_e[:, 1, 5] = mtx_e[:, 3, 6] = cc
dd = (det0 / det)[..., None, None]
mtx_ee = dd * nm.einsum('cij,qjk->cqik', tg0[:, 0], mtx_e)
mtx_be = nm.concatenate((mtx_b, mtx_ee), 3)
return mtx_be
[docs]
def rotate_elastic_tensor(mtx_d, bfu, ebs):
"""
Rotate the elastic tensor into the local coordinate system of each cell.
The local coordinate system results from interpolation of `ebs` with the
bilinear basis.
"""
# ! mtx_ts is not orthonormal for non-planar cells!
mtx_ts = nm.einsum('qi,cijk->cqjk', bfu[:, 0, :], ebs)
# print nm.einsum('cqji,cqjk->cqik', mtx_ts, mtx_ts)
# ? double entries here?
tt = create_strain_transform(mtx_ts)
mtx_dr = ddot(ddot(tt, mtx_d, 'ATB'), tt, 'AB')
return mtx_dr
[docs]
def lock_drilling_rotations(mtx, ebs, coefs):
"""
Lock the drilling rotations in the stiffness matrix.
"""
mtx_drl = create_drl_transform(ebs)
mtx_idrl = nm.linalg.inv(mtx_drl)
mtx_tr = ddot(ddot(mtx_drl, mtx), mtx_idrl)
idrl = nm.arange(20, 24)
mtx_tr[:, idrl, idrl] = coefs[:, None]
mtx2 = ddot(ddot(mtx_idrl, mtx_tr), mtx_drl)
return mtx2