{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Display Eigenmodes of Vibration {#eigenmodes_example}\n\nThis example applies the `warp_by_vector\n<pyvista.DataSetFilters.warp_by_vector>`{.interpreted-text role=\"func\"}\nfilter to a cube whose eigenmodes have been computed using the Ritz\nmethod, as outlined in Visscher, William M., Albert Migliori, Thomas M.\nBell, et Robert A. Reinert. \\\"On the normal modes of free vibration of\ninhomogeneous and anisotropic elastic objects\\\". The Journal of the\nAcoustical Society of America 90, n.4 (October 1991): 2154-62.\n<https://asa.scitation.org/doi/10.1121/1.401643>\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, let\\'s solve the eigenvalue problem for a vibrating cube. We use\na crude approximation (by choosing a low max polynomial order) to get a\nfast computation.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from itertools import product\n\nimport numpy as np\nfrom scipy.linalg import eigh\n\nimport pyvista as pv\n\n\ndef analytical_integral_rppd(p, q, r, a, b, c):\n \"\"\"Returns the analytical value of the RPPD integral, i.e. the integral\n of x**p * y**q * z**r for (x, -a, a), (y, -b, b), (z, -c, c).\"\"\"\n if p < 0 or q < 0 or r < 0.0:\n return 0.0\n else:\n return (\n a ** (p + 1)\n * b ** (q + 1)\n * c ** (r + 1)\n * ((-1) ** p + 1)\n * ((-1) ** q + 1)\n * ((-1) ** r + 1)\n / ((p + 1) * (q + 1) * (r + 1))\n )\n\n\ndef make_cijkl_E_nu(E=200, nu=0.3):\n \"\"\"Makes cijkl from E and nu.\n Default values for steel are: E=200 GPa, nu=0.3.\"\"\"\n lambd = E * nu / (1 + nu) / (1 - 2 * nu)\n mu = E / 2 / (1 + nu)\n cij = np.zeros((6, 6))\n cij[(0, 1, 2), (0, 1, 2)] = lambd + 2 * mu\n cij[(0, 0, 1, 1, 2, 2), (1, 2, 0, 2, 0, 1)] = lambd\n cij[(3, 4, 5), (3, 4, 5)] = mu\n # check symmetry\n assert np.allclose(cij, cij.T)\n # convert to order 4 tensor\n coord_mapping = {\n (1, 1): 1,\n (2, 2): 2,\n (3, 3): 3,\n (2, 3): 4,\n (1, 3): 5,\n (1, 2): 6,\n (2, 1): 6,\n (3, 1): 5,\n (3, 2): 4,\n }\n\n cijkl = np.zeros((3, 3, 3, 3))\n for i, j, k, l in product(range(3), repeat=4):\n u = coord_mapping[i + 1, j + 1]\n v = coord_mapping[k + 1, l + 1]\n cijkl[i, j, k, l] = cij[u - 1, v - 1]\n return cijkl, cij\n\n\ndef get_first_N_above_thresh(N, freqs, thresh, decimals=3):\n \"\"\"Returns first N unique frequencies with amplitude above threshold based\n on first decimals.\"\"\"\n unique_freqs, unique_indices = np.unique(np.round(freqs, decimals=decimals), return_index=True)\n nonzero = unique_freqs > thresh\n unique_freqs, unique_indices = unique_freqs[nonzero], unique_indices[nonzero]\n return unique_freqs[:N], unique_indices[:N]\n\n\ndef assemble_mass_and_stiffness(N, F, geom_params, cijkl):\n \"\"\"This routine assembles the mass and stiffness matrix.\n It first builds an index of basis functions as a quadruplet of\n component and polynomial order for (x^p, y^q, z^r) of maximum order N.\n\n This routine only builds the symmetric part of the matrix to speed\n things up.\n \"\"\"\n # building coordinates\n triplets = [\n (p, q, r) for p in range(N + 1) for q in range(N - p + 1) for r in range(N - p - q + 1)\n ]\n assert len(triplets) == (N + 1) * (N + 2) * (N + 3) // 6\n\n quadruplets = []\n for i, triplet in product(range(3), triplets):\n quadruplets.append((i, *triplet))\n assert len(quadruplets) == 3 * (N + 1) * (N + 2) * (N + 3) // 6\n\n # assembling the mass and stiffness matrix in a single loop\n R = len(triplets)\n E = np.zeros((3 * R, 3 * R)) # the mass matrix\n G = np.zeros((3 * R, 3 * R)) # the stiffness matrix\n for index1, quad1 in enumerate(quadruplets):\n I, p1, q1, r1 = quad1\n for index2, quad2 in enumerate(quadruplets[index1:]):\n index2 = index2 + index1\n J, p2, q2, r2 = quad2\n G[index1, index2] = (\n cijkl[I, 1 - 1, J, 1 - 1]\n * p1\n * p2\n * F(p1 + p2 - 2, q1 + q2, r1 + r2, **geom_params)\n + cijkl[I, 1 - 1, J, 2 - 1]\n * p1\n * q2\n * F(p1 + p2 - 1, q1 + q2 - 1, r1 + r2, **geom_params)\n + cijkl[I, 1 - 1, J, 3 - 1]\n * p1\n * r2\n * F(p1 + p2 - 1, q1 + q2, r1 + r2 - 1, **geom_params)\n + cijkl[I, 2 - 1, J, 1 - 1]\n * q1\n * p2\n * F(p1 + p2 - 1, q1 + q2 - 1, r1 + r2, **geom_params)\n + cijkl[I, 2 - 1, J, 2 - 1]\n * q1\n * q2\n * F(p1 + p2, q1 + q2 - 2, r1 + r2, **geom_params)\n + cijkl[I, 2 - 1, J, 3 - 1]\n * q1\n * r2\n * F(p1 + p2, q1 + q2 - 1, r1 + r2 - 1, **geom_params)\n + cijkl[I, 3 - 1, J, 1 - 1]\n * r1\n * p2\n * F(p1 + p2 - 1, q1 + q2, r1 + r2 - 1, **geom_params)\n + cijkl[I, 3 - 1, J, 2 - 1]\n * r1\n * q2\n * F(p1 + p2, q1 + q2 - 1, r1 + r2 - 1, **geom_params)\n + cijkl[I, 3 - 1, J, 3 - 1]\n * r1\n * r2\n * F(p1 + p2, q1 + q2, r1 + r2 - 2, **geom_params)\n )\n G[index2, index1] = G[index1, index2] # since stiffness matrix is symmetric\n if I == J:\n E[index1, index2] = F(p1 + p2, q1 + q2, r1 + r2, **geom_params)\n E[index2, index1] = E[index1, index2] # since mass matrix is symmetric\n return E, G, quadruplets\n\n\nN = 8 # maximum order of x^p y^q z^r polynomials\nrho = 8.0 # g/cm^3\nl1, l2, l3 = 0.2, 0.2, 0.2 # all in cm\ngeometry_parameters = {'a': l1 / 2.0, 'b': l2 / 2.0, 'c': l3 / 2.0}\ncijkl, cij = make_cijkl_E_nu(200, 0.3) # Gpa, without unit\nE, G, quadruplets = assemble_mass_and_stiffness(\n N,\n analytical_integral_rppd,\n geometry_parameters,\n cijkl,\n)\n\n# solving the eigenvalue problem using symmetric solver\nw, vr = eigh(a=G, b=E)\nomegas = np.sqrt(np.abs(w) / rho) * 1e5 # convert back to Hz\nfreqs = omegas / (2 * np.pi)\n# expected values from (Bernard 2014, p.14),\n# error depends on polynomial order ``N``\nexpected_freqs_kHz = np.array([704.8, 949.0, 965.2, 1096.3, 1128.4, 1182.8, 1338.9, 1360.9])\ncomputed_freqs_kHz, mode_indices = get_first_N_above_thresh(8, freqs / 1e3, thresh=1, decimals=1)\nprint('found the following first unique eigenfrequencies:')\nfor ind, (freq1, freq2) in enumerate(zip(computed_freqs_kHz, expected_freqs_kHz)):\n error = np.abs(freq2 - freq1) / freq1 * 100.0\n print(f\"freq. {ind + 1:1}: {freq1:8.1f} kHz, expected: {freq2:8.1f} kHz, error: {error:.2f} %\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let\\'s display a mode on a mesh of the cube.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create the 3D NumPy array of spatially referenced data\n# (nx by ny by nz)\nnx, ny, nz = 30, 31, 32\n\nx = np.linspace(-l1 / 2.0, l1 / 2.0, nx)\ny = np.linspace(-l2 / 2.0, l2 / 2.0, ny)\nx, y = np.meshgrid(x, y)\nz = np.zeros_like(x) + l3 / 2.0\ngrid = pv.StructuredGrid(x, y, z)\n\nslices = []\nfor zz in np.linspace(-l3 / 2.0, l3 / 2.0, nz)[::-1]:\n slice_ = grid.points.copy()\n slice_[:, -1] = zz\n slices.append(slice_)\n\nvol = pv.StructuredGrid()\nvol.points = np.vstack(slices)\nvol.dimensions = [*grid.dimensions[0:2], nz]\n\nfor i, mode_index in enumerate(mode_indices):\n eigenvector = vr[:, mode_index]\n displacement_points = np.zeros_like(vol.points)\n for weight, (component, p, q, r) in zip(eigenvector, quadruplets):\n displacement_points[:, component] += (\n weight * vol.points[:, 0] ** p * vol.points[:, 1] ** q * vol.points[:, 2] ** r\n )\n if displacement_points.max() > 0.0:\n displacement_points /= displacement_points.max()\n vol[f'eigenmode_{i:02}'] = displacement_points\n\nwarpby = 'eigenmode_00'\nwarped = vol.warp_by_vector(warpby, factor=0.04)\nwarped.translate([-1.5 * l1, 0.0, 0.0], inplace=True)\npl = pv.Plotter()\npl.add_mesh(vol, style='wireframe', scalars=warpby, show_scalar_bar=False)\npl.add_mesh(warped, scalars=warpby)\npl.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, let\\'s make a gallery of the first 8 unique eigenmodes.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"pl = pv.Plotter(shape=(2, 4))\nfor i, j in product(range(2), range(4)):\n pl.subplot(i, j)\n current_index = 4 * i + j\n vector = f\"eigenmode_{current_index:02}\"\n pl.add_text(\n f\"mode {current_index}, freq. {computed_freqs_kHz[current_index]:.1f} kHz\",\n font_size=10,\n )\n pl.add_mesh(vol.warp_by_vector(vector, factor=0.03), scalars=vector, show_scalar_bar=False)\npl.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}