Shape Function Generation#

You can get all the necessary functions for interpolation by calling the generate_class_functions method of any class. It returns the shape functions as symbolic matrices and also functions for fast numerical evaluation (red the corresponding docs for further information).

Note These functions are automatically generated runtime, but they can be useful if you want to build something on your own, or for educational/publication purposes.

Note Generation of shape functions and their derivatives is illustrated using 8-noded hexahedron cells, but the usage applies for all cells of the library.

[13]:
from sigmaepsilon.mesh.cells import H8

shp, dshp, shpf, shpmf, dshpf = H8.Geometry.generate_class_functions()

The first two of the returned items are symbolic matrices for the shape functions and their derivatives:

[14]:
shp
[14]:
$\displaystyle \left[\begin{matrix}- 0.125 r s t + 0.125 r s + 0.125 r t - 0.125 r + 0.125 s t - 0.125 s - 0.125 t + 0.125\\0.125 r s t - 0.125 r s - 0.125 r t + 0.125 r + 0.125 s t - 0.125 s - 0.125 t + 0.125\\- 0.125 r s t + 0.125 r s - 0.125 r t + 0.125 r - 0.125 s t + 0.125 s - 0.125 t + 0.125\\0.125 r s t - 0.125 r s + 0.125 r t - 0.125 r - 0.125 s t + 0.125 s - 0.125 t + 0.125\\0.125 r s t + 0.125 r s - 0.125 r t - 0.125 r - 0.125 s t - 0.125 s + 0.125 t + 0.125\\- 0.125 r s t - 0.125 r s + 0.125 r t + 0.125 r - 0.125 s t - 0.125 s + 0.125 t + 0.125\\0.125 r s t + 0.125 r s + 0.125 r t + 0.125 r + 0.125 s t + 0.125 s + 0.125 t + 0.125\\- 0.125 r s t - 0.125 r s - 0.125 r t - 0.125 r + 0.125 s t + 0.125 s + 0.125 t + 0.125\end{matrix}\right]$
[15]:
dshp
[15]:
$\displaystyle \left[\begin{matrix}- 0.125 s t + 0.125 s + 0.125 t - 0.125 & - 0.125 r t + 0.125 r + 0.125 t - 0.125 & - 0.125 r s + 0.125 r + 0.125 s - 0.125\\0.125 s t - 0.125 s - 0.125 t + 0.125 & 0.125 r t - 0.125 r + 0.125 t - 0.125 & 0.125 r s - 0.125 r + 0.125 s - 0.125\\- 0.125 s t + 0.125 s - 0.125 t + 0.125 & - 0.125 r t + 0.125 r - 0.125 t + 0.125 & - 0.125 r s - 0.125 r - 0.125 s - 0.125\\0.125 s t - 0.125 s + 0.125 t - 0.125 & 0.125 r t - 0.125 r - 0.125 t + 0.125 & 0.125 r s + 0.125 r - 0.125 s - 0.125\\0.125 s t + 0.125 s - 0.125 t - 0.125 & 0.125 r t + 0.125 r - 0.125 t - 0.125 & 0.125 r s - 0.125 r - 0.125 s + 0.125\\- 0.125 s t - 0.125 s + 0.125 t + 0.125 & - 0.125 r t - 0.125 r - 0.125 t - 0.125 & - 0.125 r s + 0.125 r - 0.125 s + 0.125\\0.125 s t + 0.125 s + 0.125 t + 0.125 & 0.125 r t + 0.125 r + 0.125 t + 0.125 & 0.125 r s + 0.125 r + 0.125 s + 0.125\\- 0.125 s t - 0.125 s - 0.125 t - 0.125 & - 0.125 r t - 0.125 r + 0.125 t + 0.125 & - 0.125 r s - 0.125 r + 0.125 s + 0.125\end{matrix}\right]$

The last three of the returned items are functions for the numerical evaluation of shape functions, the shape function matrix and the shape function derivatives.

[16]:
coords = H8.Geometry.master_coordinates()
coords
[16]:
array([[-1., -1., -1.],
       [ 1., -1., -1.],
       [ 1.,  1., -1.],
       [-1.,  1., -1.],
       [-1., -1.,  1.],
       [ 1., -1.,  1.],
       [ 1.,  1.,  1.],
       [-1.,  1.,  1.]])

The coordinates of the master element is described by a 8x3 matrix, since the cell has 8 nodes and three spatial dimensions.

[17]:
coords.shape
[17]:
(8, 3)

Evaluate these functions using the first two coordinates of the master element to get an idea about the shapes of the resulting arrays (see the docs for the details).

[18]:
shpf(coords[:2]).shape
[18]:
(2, 8)
[19]:
shpmf(coords[:2]).shape
[19]:
(2, 1, 8)
[20]:
dshpf(coords[:2]).shape
[20]:
(2, 8, 3)

If you want to evaluate shape function derivatives wrt. to the local coordinate frames of cells, you need to generate a mesh first:

[21]:
import numpy as np
from sigmaepsilon.mesh import PolyData, PointData, CartesianFrame
from sigmaepsilon.mesh.grid import grid
from sigmaepsilon.mesh.cells import H8

size = Lx, Ly, Lz = 800, 600, 100
shape = nx, ny, nz = 8, 6, 2
xbins = np.linspace(0, Lx, nx + 1)
ybins = np.linspace(0, Ly, ny + 1)
zbins = np.linspace(0, Lz, nz + 1)
bins = xbins, ybins, zbins
coords, topo = grid(bins=bins, eshape="H8")
frame = CartesianFrame(dim=3)

pd = PointData(coords=coords)
cd = H8(topo=topo, frames=frame)
mesh = PolyData(pd, cd)

Observing the shape of the topology array helps to understand the shapes of the arrays of the following blocks.

[22]:
topo.shape
[22]:
(96, 8)

To calculate the derivatives, you need the jacobian matrices of the cells evaluated at certain points. The following block evaluates the jacobian matrices of the cells at the nodes of the cells.

[23]:
jac = cd.jacobian_matrix()
jac.shape
[23]:
(96, 8, 3, 3)

Then, you can calculate the derivatives wrt. the local frames of the cells like this:

[24]:
pcoords = H8.Geometry.master_coordinates()
gdshp = H8.Geometry.shape_function_derivatives(pcoords[:2], jac=jac)
gdshp.shape
[24]:
(96, 2, 8, 3)