Source code for sfepy.mechanics.shell10x

"""
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_transformation_matrix(coors): r""" Create a transposed coordinate transformation matrix, that transforms 3D coordinates of quadrilateral cell vertices so that the transformed vertices of a plane cell are in the :math:`x-y` plane. The rotation is performed w.r.t. the centres of quadrilaterals. Parameters ---------- coors : array The coordinates of cell vertices, shape `(n_el, 4, 3)`. Returns ------- mtx_t : array The transposed transformation matrix :math:`T`, i.e. :math:`X_{inplane} = T^T X_{3D}`. Notes ----- :math:`T = [t_1, t_2, n]`, where :math:`t_1`, :math:`t_2`, are unit in-plane (column) vectors and :math:`n` is the unit normal vector, all mutually orthonormal. """ # Local coordinate system. t1 = coors[:, 1, :] + coors[:, 2, :] - coors[:, 0, :] - coors[:, 3, :] t2 = coors[:, 2, :] + coors[:, 3, :] - coors[:, 0, :] - coors[:, 1, :] n = nm.cross(t1, t2) t2 = nm.cross(n, t1) t1 = t1 / norm(t1)[:, None] t2 = t2 / norm(t2)[:, None] n = n / norm(n)[:, None] # Coordinate transformation matrix (transposed!). mtx_t = nm.concatenate((t1[:, :, None], t2[:, :, None], n[:, :, None]), axis=2) return mtx_t
[docs] def transform_asm_matrices(out, mtx_t, blocks): """ Transform matrix assembling contributions to global coordinate system, one node at a time. Parameters ---------- out : array The array of matrices, transformed in-place. mtx_t : array The array of transposed transformation matrices :math:`T`, see :func:`create_transformation_matrix`. blocks : array The DOF blocks that are """ for ir in blocks: ir = ir[:, None] for ic in blocks: fn = out[:, 0, ir, ic] out[:, 0, ir, ic] = ddot(ddot(mtx_t, fn, 'AB'), mtx_t, 'ABT')
[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 create_strain_transform(mtx_ts): r""" Create strain tensor transformation matrices, given coordinate transformation matrices. Notes ----- Expresses :math:`T E T^T` in terms of symmetrix storage as :math:`Q e`, with the ordering of components: :math:`e = [e_{11}, e_{22}, e_{33}, 2 e_{12}, 2 e_{13}, 2 e_{23}]`. """ mtx_qs = nm.empty((mtx_ts.shape[0], mtx_ts.shape[1], 6, 6), dtype=nm.float64) for ie in range(mtx_ts.shape[0]): for iq in range(mtx_ts.shape[1]): l = mtx_ts[ie, iq, :, 0] m = mtx_ts[ie, iq, :, 1] n = mtx_ts[ie, iq, :, 2] mtx_q = [ [l[0]**2, m[0]**2, n[0]**2, l[0]*m[0], n[0]*l[0], m[0]*n[0]], [l[1]**2, m[1]**2, n[1]**2, l[1]*m[1], n[1]*l[1], m[1]*n[1]], [l[2]**2, m[2]**2, n[2]**2, l[2]*m[2], n[2]*l[2], m[2]*n[2]], [2*l[0]*l[1], 2*m[0]*m[1], 2*n[0]*n[1], l[0]*m[1]+l[1]*m[0], n[0]*l[1]+n[1]*l[0], m[0]*n[1]+m[1]*n[0]], [2*l[2]*l[0], 2*m[2]*m[0], 2*n[2]*n[0], l[2]*m[0]+l[0]*m[2], n[2]*l[0]+n[0]*l[2], m[2]*n[0]+m[0]*n[2]], [2*l[1]*l[2], 2*m[1]*m[2], 2*n[1]*n[2], l[1]*m[2]+l[2]*m[1], n[1]*l[2]+n[2]*l[1], m[1]*n[2]+m[2]*n[1]], ] mtx_qs[ie, iq] = mtx_q return mtx_qs
[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 create_drl_transform(ebs): """ Create the transformation matrix for locking of the drilling rotations. """ n_el = ebs.shape[0] mtx_drl = nm.tile(nm.eye(24, dtype=nm.float64), (n_el, 1, 1)) for ii in range(4): nh = ebs[:, ii, -1, :] mtx = nm.c_[nm.c_[1.0 - nh[:, 0]**2, -nh[:, 0]*nh[:, 1], -nh[:, 0]*nh[:, 2]], nm.c_[-nh[:, 0]*nh[:, 1], 1.0 - nh[:, 1]**2, -nh[:, 1]*nh[:, 2]], nh].reshape(n_el, 3, 3) mtx_drl[:, 12+ii:12+ii+9:4, 12+ii:12+ii+9:4] = mtx return mtx_drl
[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