Note
Click here to download the full example code
Displaying eigenmodes of vibration using warp_by_vector
¶
This example applies the warp_by_vector
filter to a cube whose eigenmodes
have been computed using the Ritz method, as outlined in Visscher, William M.,
Albert Migliori, Thomas M. Bell, et Robert A. Reinert. « On the normal modes
of free vibration of inhomogeneous and anisotropic elastic objects ». The
Journal of the Acoustical Society of America 90, nᵒ 4 (october 1991): 2154‑62.
https://doi.org/10.1121/1.401643.
First, let’s solve the eigenvalue problem for a vibrating cube. We use a crude approximation (by choosing a low max polynomial order) to get a fast computation.
import numpy as np
from scipy.linalg import eigh
import pyvista as pv
def analytical_integral_rppd(p, q, r, a, b, c):
"""Returns the analytical value of the RPPD integral, i.e. the integral
of x**p * y**q * z**r for (x, -a, a), (y, -b, b), (z, -c, c)."""
if p < 0:
return 0.
elif q < 0:
return 0.
elif r < 0.:
return 0.
else:
return a ** (p + 1) * b ** (q + 1) * c ** (r + 1) * \
((-1) ** p + 1) * ((-1) ** q + 1) * ((-1) ** r + 1) \
/ ((p + 1) * (q + 1) * (r + 1))
def make_cijkl_E_nu(E=200, nu=0.3):
"""Makes cijkl from E and nu.
Default values for steel are: E=200 GPa, nu=0.3."""
lambd = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)
cij = np.zeros((6, 6))
cij[(0, 1, 2), (0, 1, 2)] = lambd + 2 * mu
cij[(0, 0, 1, 1, 2, 2), (1, 2, 0, 2, 0, 1)] = lambd
cij[(3, 4, 5), (3, 4, 5)] = mu
# check symmetry
assert np.allclose(cij, cij.T)
# convert to order 4 tensor
coord_mapping = {(1, 1): 1,
(2, 2): 2,
(3, 3): 3,
(2, 3): 4,
(1, 3): 5,
(1, 2): 6,
(2, 1): 6,
(3, 1): 5,
(3, 2): 4}
cijkl = np.zeros((3, 3, 3, 3))
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
u = coord_mapping[(i + 1, j + 1)]
v = coord_mapping[(k + 1, l + 1)]
cijkl[i, j, k, l] = cij[u - 1, v - 1]
return cijkl, cij
def get_first_N_above_thresh(N, freqs, thresh, decimals=3):
"""Returns first N unique frequencies with amplitude above threshold based
on first decimals."""
unique_freqs, unique_indices = np.unique(
np.round(freqs, decimals=decimals), return_index=True)
nonzero = unique_freqs > thresh
unique_freqs, unique_indices = unique_freqs[nonzero], unique_indices[
nonzero]
return unique_freqs[:N], unique_indices[:N]
def assemble_mass_and_stiffness(N, F, geom_params, cijkl):
"""This routine assembles the mass and stiffness matrix.
It first builds an index of basis functions as a quadruplet of
component and polynomial order for (x^p, y^q, z^r) of maximum order N.
This routine only builds the symmetric part of the matrix to speed
things up.
"""
# building coordinates
triplets = []
for p in range(N + 1):
for q in range(N - p + 1):
for r in range(N - p - q + 1):
triplets.append((p, q, r))
assert len(triplets) == (N + 1) * (N + 2) * (N + 3) // 6
quadruplets = []
for i in range(3):
for triplet in triplets:
quadruplets.append((i, *triplet))
assert len(quadruplets) == 3 * (N + 1) * (N + 2) * (N + 3) // 6
# assembling the mass and stiffness matrix in a single loop
R = len(triplets)
E = np.zeros((3 * R, 3 * R)) # the mass matrix
G = np.zeros((3 * R, 3 * R)) # the stiffness matrix
for index1, quad1 in enumerate(quadruplets):
I, p1, q1, r1 = quad1
for index2, quad2 in enumerate(quadruplets[index1:]):
index2 = index2 + index1
J, p2, q2, r2 = quad2
G[index1, index2] = cijkl[I, 1 - 1, J, 1 - 1] * p1 * p2 * F(
p1 + p2 - 2, q1 + q2, r1 + r2, **geom_params) + \
cijkl[I, 1 - 1, J, 2 - 1] * p1 * q2 * F(
p1 + p2 - 1, q1 + q2 - 1, r1 + r2,
**geom_params) + \
cijkl[I, 1 - 1, J, 3 - 1] * p1 * r2 * F(
p1 + p2 - 1, q1 + q2, r1 + r2 - 1,
**geom_params) + \
cijkl[I, 2 - 1, J, 1 - 1] * q1 * p2 * F(
p1 + p2 - 1, q1 + q2 - 1, r1 + r2,
**geom_params) + \
cijkl[I, 2 - 1, J, 2 - 1] * q1 * q2 * F(
p1 + p2, q1 + q2 - 2, r1 + r2, **geom_params) + \
cijkl[I, 2 - 1, J, 3 - 1] * q1 * r2 * F(
p1 + p2, q1 + q2 - 1, r1 + r2 - 1,
**geom_params) + \
cijkl[I, 3 - 1, J, 1 - 1] * r1 * p2 * F(
p1 + p2 - 1, q1 + q2, r1 + r2 - 1,
**geom_params) + \
cijkl[I, 3 - 1, J, 2 - 1] * r1 * q2 * F(
p1 + p2, q1 + q2 - 1, r1 + r2 - 1,
**geom_params) + \
cijkl[I, 3 - 1, J, 3 - 1] * r1 * r2 * F(
p1 + p2, q1 + q2, r1 + r2 - 2, **geom_params)
G[index2, index1] = G[
index1, index2] # since stiffness matrix is symmetric
if I == J:
E[index1, index2] = F(p1 + p2, q1 + q2, r1 + r2, **geom_params)
E[index2, index1] = E[
index1, index2] # since mass matrix is symmetric
return E, G, quadruplets
N = 8 # maximum order of x^p y^q z^r polynomials
rho = 8.0 # g/cm^3
l1, l2, l3 = .2, .2, .2 # all in cm
geometry_parameters = {'a': l1 / 2., 'b': l2 / 2., 'c': l3 / 2.}
cijkl, cij = make_cijkl_E_nu(200, 0.3) # Gpa, without unit
E, G, quadruplets = assemble_mass_and_stiffness(N, analytical_integral_rppd,
geometry_parameters, cijkl)
# solving the eigenvalue problem using symmetric solver
w, vr = eigh(a=G, b=E)
omegas = np.sqrt(np.abs(w) / rho) * 1e5 # convert back to Hz
freqs = omegas / (2 * np.pi)
# expected values from (Bernard 2014, p.14),
# error depends on polynomial order ``N``
expected_freqs_kHz = np.array(
[704.8, 949., 965.2, 1096.3, 1128.4, 1182.8, 1338.9, 1360.9])
computed_freqs_kHz, mode_indices = get_first_N_above_thresh(8, freqs / 1e3,
thresh=1,
decimals=1)
print('found the following first unique eigenfrequencies:')
for ind, (freq1, freq2) in enumerate(
zip(computed_freqs_kHz, expected_freqs_kHz)):
error = np.abs(freq2 - freq1) / freq1 * 100.
print(
f"freq. {ind + 1:1}: {freq1:8.1f} kHz," + \
f" expected: {freq2:8.1f} kHz, error: {error:.2f} %")
Out:
found the following first unique eigenfrequencies:
freq. 1: 705.1 kHz, expected: 704.8 kHz, error: 0.04 %
freq. 2: 949.1 kHz, expected: 949.0 kHz, error: 0.01 %
freq. 3: 965.7 kHz, expected: 965.2 kHz, error: 0.05 %
freq. 4: 1096.3 kHz, expected: 1096.3 kHz, error: 0.00 %
freq. 5: 1128.6 kHz, expected: 1128.4 kHz, error: 0.02 %
freq. 6: 1183.9 kHz, expected: 1182.8 kHz, error: 0.09 %
freq. 7: 1339.0 kHz, expected: 1338.9 kHz, error: 0.01 %
freq. 8: 1361.8 kHz, expected: 1360.9 kHz, error: 0.07 %
Now, let’s display a mode on a mesh of the cube.
# Create the 3D NumPy array of spatially referenced data
# (nx by ny by nz)
nx, ny, nz = 30, 31, 32
x = np.linspace(-l1 / 2., l1 / 2., nx)
y = np.linspace(-l2 / 2., l2 / 2., ny)
x, y = np.meshgrid(x, y)
z = np.zeros_like(x) + l3 / 2.
grid = pv.StructuredGrid(x, y, z)
slices = []
for zz in np.linspace(-l3 / 2., l3 / 2., nz)[::-1]:
slice = grid.points.copy()
slice[:, -1] = zz
slices.append(slice)
vol = pv.StructuredGrid()
vol.points = np.vstack(slices)
vol.dimensions = [*grid.dimensions[0:2], nz]
for i, mode_index in enumerate(mode_indices):
eigenvector = vr[:, mode_index]
displacement_points = np.zeros_like(vol.points)
for weight, (component, p, q, r) in zip(eigenvector, quadruplets):
displacement_points[:, component] += weight * vol.points[:, 0] ** p * \
vol.points[:, 1] ** q * \
vol.points[:, 2] ** r
if displacement_points.max() > 0.:
displacement_points /= displacement_points.max()
vol[f'eigenmode_{i:02}'] = displacement_points
warpby = 'eigenmode_00'
warped = vol.warp_by_vector(warpby, factor=0.04)
warped.translate([-1.5 * l1, 0., 0.])
p = pv.Plotter()
p.add_mesh(vol, style='wireframe', scalars=warpby)
p.add_mesh(warped, scalars=warpby)
p.show()

Out:
[(0.5624591927763861, 0.7150308248544035, 0.7150308248544028),
(-0.15257163207801672, 7.216449660063518e-16, 0.0),
(0.0, 0.0, 1.0)]
Finally, let’s make a gallery of the first 8 unique eigenmodes.
p = pv.Plotter(shape=(2, 4))
for i in range(2):
for j in range(4):
p.subplot(i, j)
current_index = 4 * i + j
vector = f"eigenmode_{current_index:02}"
p.add_text(
f"mode {current_index}," + \
f" freq. {computed_freqs_kHz[current_index]:.1f} kHz",
font_size=10)
p.add_mesh(vol.warp_by_vector(vector, factor=0.03), scalars=vector)
p.show()

Out:
[(0.7392698613030427, 0.7392698613030582, 0.7392698613030666),
(-1.8957058145474548e-14, -3.594347042223944e-15, 4.884981308350689e-15),
(0.0, 0.0, 1.0)]
Total running time of the script: ( 0 minutes 9.347 seconds)