Source code for sfepy.mechanics.membranes

from __future__ import absolute_import
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.discrete import PolySpace
from sfepy.discrete.fem.mappings import VolumeMapping
from sfepy.mechanics.tensors import dim2sym
from six.moves import range

[docs]def create_transformation_matrix(coors):
"""
Create a transposed coordinate transformation matrix, that
transforms 3D coordinates of element face nodes so that the
transformed nodes are in the x-y plane. The rotation is performed
w.r.t. the first node of each face.

Parameters
----------
coors : array
The coordinates of element nodes, shape (n_el, n_ep, dim).

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[:, 0, :]
t2 = coors[:, -1, :] - coors[:, 0, :]
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_vectors(out, mtx_t):
"""
Transform vector assembling contributions to global coordinate system, one
node at a time.

Parameters
----------
out : array
The array of vectors, transformed in-place.
mtx_t : array
The transposed transformation matrix :math:T, see
:func:create_transformation_matrix.
"""
n_ep = out.shape[2] // mtx_t.shape[2]
for iep in range(n_ep):
ir = slice(iep, None, n_ep)
fn = out[:, 0, ir, 0]
fn[:] = dot_sequences(mtx_t, fn, 'AB')

[docs]def transform_asm_matrices(out, mtx_t):
"""
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 transposed transformation matrix :math:T, see
:func:create_transformation_matrix.
"""
n_ep = out.shape[-1] // mtx_t.shape[-1]
for iepr in range(n_ep):
ir = slice(iepr, None, n_ep)
for iepc in range(n_ep):
ic = slice(iepc, None, n_ep)
fn = out[:, 0, ir, ic]
fn[:] = dot_sequences(dot_sequences(mtx_t, fn, 'AB'), mtx_t, 'ABT')

[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(field, region, integral):
"""
Describe membrane geometry in a given region.

Parameters
----------
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(region, integral, 'surface')
sd = field.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)
membrane_geo.bf[:] = ps.eval_base(qp[0])

return mtx_t, membrane_geo

[docs]def describe_deformation(el_disps, bfg):
"""
Describe deformation of a thin incompressible 2D membrane in 3D
space, composed of flat finite element faces.

The coordinate system of each element (face), i.e. the membrane
mid-surface, should coincide with the x, y axes of the x-y
plane.

Parameters
----------
el_disps : array
The displacements of element nodes, shape (n_el, n_ep, dim).
bfg : array
The in-plane base function gradients, shape (n_el, n_qp, dim-1,
n_ep).

Returns
-------
mtx_c ; array
The in-plane right Cauchy-Green deformation tensor
:math:C_{ij}, :math:i, j = 1, 2.
c33 : array
The component :math:C_{33} computed from the incompressibility
condition.
mtx_b : array
The discrete Green strain variation operator.
"""
sh = bfg.shape
n_ep = sh[3]

dim = el_disps.shape[2]
sym2 = dim2sym(dim-1)

# Repeat el_disps by number of quadrature points.
el_disps_qp = insert_strided_axis(el_disps, 1, bfg.shape[1])

# Transformed (in-plane) displacement gradient with
# shape (n_el, n_qp, 2 (-> a), 3 (-> i)), du_i/dX_a.
du = dot_sequences(bfg, el_disps_qp)

# Deformation gradient F w.r.t. in plane coordinates.
# F_{ia} = dx_i / dX_a,
# a \in {1, 2} (rows), i \in {1, 2, 3} (columns).
mtx_f = du + nm.eye(dim - 1, dim, dtype=du.dtype)

# Right Cauchy-Green deformation tensor C.
# C_{ab} = F_{ka} F_{kb}, a, b \in {1, 2}.
mtx_c = dot_sequences(mtx_f, mtx_f, 'ABT')

# C_33 from incompressibility.
c33 = 1.0 / (mtx_c[..., 0, 0] * mtx_c[..., 1, 1]
- mtx_c[..., 0, 1]**2)

# Discrete Green strain variation operator.
mtx_b = nm.empty((sh[0], sh[1], sym2, dim * n_ep), dtype=nm.float64)
mtx_b[..., 0, 0*n_ep:1*n_ep] = bfg[..., 0, :] * mtx_f[..., 0, 0:1]
mtx_b[..., 0, 1*n_ep:2*n_ep] = bfg[..., 0, :] * mtx_f[..., 0, 1:2]
mtx_b[..., 0, 2*n_ep:3*n_ep] = bfg[..., 0, :] * mtx_f[..., 0, 2:3]
mtx_b[..., 1, 0*n_ep:1*n_ep] = bfg[..., 1, :] * mtx_f[..., 1, 0:1]
mtx_b[..., 1, 1*n_ep:2*n_ep] = bfg[..., 1, :] * mtx_f[..., 1, 1:2]
mtx_b[..., 1, 2*n_ep:3*n_ep] = bfg[..., 1, :] * mtx_f[..., 1, 2:3]
mtx_b[..., 2, 0*n_ep:1*n_ep] = bfg[..., 1, :] * mtx_f[..., 0, 0:1] \
+ bfg[..., 0, :] * mtx_f[..., 1, 0:1]
mtx_b[..., 2, 1*n_ep:2*n_ep] = bfg[..., 0, :] * mtx_f[..., 1, 1:2] \
+ bfg[..., 1, :] * mtx_f[..., 0, 1:2]
mtx_b[..., 2, 2*n_ep:3*n_ep] = bfg[..., 0, :] * mtx_f[..., 1, 2:3] \
+ bfg[..., 1, :] * mtx_f[..., 0, 2:3]

return mtx_c, c33, mtx_b

[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