import numpy as nm
from sfepy.base.base import assert_
from sfepy.linalg import norm_l2_along_axis as norm
from sfepy.linalg import dot_sequences, insert_strided_axis
from sfepy.fem.poly_spaces import PolySpace
from sfepy.fem.mappings import VolumeMapping
from sfepy.mechanics.tensors import dim2sym
[docs]def create_mapping(coors, gel, order):
"""
Create mapping from transformed (in `x-y` plane) element faces to
reference element faces.
Parameters
----------
coors : array
The transformed coordinates of element nodes, shape `(n_el,
n_ep, dim)`. The function verifies that the all `z` components
are zero.
gel : GeometryElement instance
The geometry element corresponding to the faces.
order : int
The polynomial order of the mapping.
Returns
-------
mapping : VolumeMapping instance
The reference element face mapping.
"""
# Strip 'z' component (should be 0 now...).
assert_(nm.allclose(coors[:, :, -1], 0.0, rtol=1e-12, atol=1e-12))
coors = coors[:, :, :-1].copy()
# Mapping from transformed element to reference element.
sh = coors.shape
seq_coors = coors.reshape((sh[0] * sh[1], sh[2]))
seq_conn = nm.arange(seq_coors.shape[0], dtype=nm.int32)
seq_conn.shape = sh[:2]
mapping = VolumeMapping(seq_coors, seq_conn, gel=gel, order=1)
return mapping
[docs]def describe_geometry(ig, field, region, integral):
"""
Describe membrane geometry in a given region.
Parameters
----------
ig : int
The element group index.
field : Field instance
The field defining the FE approximation.
region : Region instance
The surface region to describe.
integral : Integral instance
The integral defining the quadrature points.
Returns
-------
mtx_t : array
The transposed transformation matrix :math:`T`, see
:func:`create_transformation_matrix`.
membrane_geo : CMapping instance
The mapping from transformed elements to a reference elements.
"""
# Coordinates of element vertices.
sg, _ = field.get_mapping(ig, region, integral, 'surface')
sd = field.aps[ig].surface_data[region.name]
coors = field.coors[sd.econn[:, :sg.n_ep]]
# Coordinate transformation matrix (transposed!).
mtx_t = create_transformation_matrix(coors)
# Transform coordinates to the local coordinate system.
coors_loc = dot_sequences((coors - coors[:, 0:1, :]), mtx_t)
# Mapping from transformed elements to reference elements.
gel = field.gel.surface_facet
vm = create_mapping(coors_loc, gel, 1)
qp = integral.get_qp(gel.name)
ps = PolySpace.any_from_args(None, gel, field.approx_order)
membrane_geo = vm.get_mapping(qp[0], qp[1], poly_space=ps)
return mtx_t, membrane_geo
[docs]def get_tangent_stress_matrix(stress, bfg):
"""
Get the tangent stress matrix of a thin incompressible 2D membrane
in 3D space, given a stress.
Parameters
----------
stress : array
The components `11, 22, 12` of the second Piola-Kirchhoff stress
tensor, shape `(n_el, n_qp, 3, 1)`.
bfg : array
The in-plane base function gradients, shape `(n_el, n_qp, dim-1,
n_ep)`.
Returns
-------
mtx : array
The tangent stress matrix, shape `(n_el, n_qp, dim*n_ep, dim*n_ep)`.
"""
n_el, n_qp, dim, n_ep = bfg.shape
dim += 1
mtx = nm.zeros((n_el, n_qp, dim * n_ep, dim * n_ep), dtype=nm.float64)
g1tg1 = dot_sequences(bfg[..., 0:1, :], bfg[..., 0:1, :], 'ATB')
g1tg2 = dot_sequences(bfg[..., 0:1, :], bfg[..., 1:2, :], 'ATB')
g2tg1 = dot_sequences(bfg[..., 1:2, :], bfg[..., 0:1, :], 'ATB')
g2tg2 = dot_sequences(bfg[..., 1:2, :], bfg[..., 1:2, :], 'ATB')
aux = stress[..., 0:1, :] * g1tg1 + stress[..., 2:3, :] * g1tg2 \
+ stress[..., 2:3, :] * g2tg1 + stress[..., 1:2, :] * g2tg2
mtx[..., 0 * n_ep : 1 * n_ep, 0 * n_ep : 1 * n_ep] = aux
mtx[..., 1 * n_ep : 2 * n_ep, 1 * n_ep : 2 * n_ep] = aux
mtx[..., 2 * n_ep : 3 * n_ep, 2 * n_ep : 3 * n_ep] = aux
return mtx
[docs]def get_invariants(mtx_c, c33):
"""
Get the first and second invariants of the right Cauchy-Green
deformation tensor describing deformation of an incompressible
membrane.
Parameters
----------
mtx_c ; array
The in-plane right Cauchy-Green deformation tensor
:math:`C_{ij}`, :math:`i, j = 1, 2`, shape `(n_el, n_qp, dim-1,
dim-1)`.
c33 : array
The component :math:`C_{33}` computed from the incompressibility
condition, shape `(n_el, n_qp)`.
Returns
-------
i1 : array
The first invariant of :math:`C_{ij}`.
i2 : array
The second invariant of :math:`C_{ij}`.
"""
i1 = mtx_c[..., 0, 0] + mtx_c[..., 1, 1] + c33
i2 = mtx_c[..., 0, 0] * mtx_c[..., 1,1] \
+ mtx_c[..., 1, 1] * c33 \
+ mtx_c[..., 0, 0] * c33 \
- mtx_c[..., 0, 1]**2
return i1, i2
[docs]def get_green_strain_sym3d(mtx_c, c33):
r"""
Get the 3D Green strain tensor in symmetric storage.
Parameters
----------
mtx_c ; array
The in-plane right Cauchy-Green deformation tensor
:math:`C_{ij}`, :math:`i, j = 1, 2`, shape `(n_el, n_qp, dim-1,
dim-1)`.
c33 : array
The component :math:`C_{33}` computed from the incompressibility
condition, shape `(n_el, n_qp)`.
Returns
-------
mtx_e : array
The membrane Green strain :math:`E_{ij} = \frac{1}{2} (C_{ij}) -
\delta_{ij}`, symmetric storage: items (11, 22, 33, 12, 13, 23),
shape `(n_el, n_qp, sym, 1)`.
"""
n_el, n_qp, dm, _ = mtx_c.shape
dim = dm + 1
sym = dim2sym(dim)
mtx_e = nm.empty((n_el, n_qp, sym, 1), dtype=mtx_c.dtype)
mtx_e[..., 0, 0] = 0.5 * (mtx_c[..., 0, 0] - 1.0)
mtx_e[..., 1, 0] = 0.5 * (mtx_c[..., 1, 1] - 1.0)
mtx_e[..., 2, 0] = 0.5 * (c33 - 1.0)
mtx_e[..., 3, 0] = 0.5 * mtx_c[..., 0, 1]
mtx_e[..., 4:, 0] = 0.0
return mtx_e
