{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Cyclic-symmetry mode family across every harmonic\n\nA 12-sector annular disk's modal spectrum decomposes per\nharmonic index $k \\in \\{0, 1, \\ldots, N/2\\}$.  Each\nharmonic's lowest non-rigid mode picks out a distinct nodal\ndiameter pattern: $k = 0$ is the breathing\n(axisymmetric) mode, $k = 1$ is the rigid-tilt-like\nmode, $k = 2$ is the diametrical-2 (two stationary\nradial nodes), $k = 3$ adds a third nodal diameter,\nand so on up to $k = N / 2 = 6$, the highest harmonic\nsustainable on a 12-sector rotor.\n\nThis example sweeps every harmonic, picks the lowest non-\nrigid mode at each $k$, and renders all seven side-by-\nside in a single 4\u00d72 viewport grid \u2014 the **cyclic-symmetry\nmode family** of the disk.\n\n## References\n\n* Thomas, D. L. (1979) \"Dynamics of rotationally periodic\n  structures,\" *J. Sound Vib.* 66 (4), 585\u2013597.\n* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,\n  \u00a710.3.4 (cyclic-symmetry modal).\n* Wildheim, S. J. (1979) \"Excitation of rotationally periodic\n  structures,\" *J. Appl. Mech.* 46, 891\u2013893.\n* Singh, M. P. and Vakakis, A. F. (1993) \"On the dynamics of\n  cyclic-symmetric structures,\" *Mechanism and Machine\n  Theory* 28 (5), 627\u2013637 (mode classification by harmonic\n  index).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport numpy as np\nimport pyvista as pv\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Geometry \u2014 12-sector annular disk\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_SECTORS = 12\nALPHA = 2.0 * np.pi / N_SECTORS\n\nR_INNER = 0.30\nR_OUTER = 0.55\nT_DISK = 0.02\nN_T, N_R, N_Z = 6, 4, 1\n\nE = 2.0e11\nNU = 0.3\nRHO = 7850.0\n\ntheta = np.linspace(0.0, ALPHA, N_T + 1)\nr = np.linspace(R_INNER, R_OUTER, N_R + 1)\nz = np.linspace(0.0, T_DISK, N_Z + 1)\nR, T, Z = np.meshgrid(r, theta, z, indexing=\"ij\")\nX = R * np.cos(T)\nY = R * np.sin(T)\ngrid = pv.StructuredGrid(X, Y, Z).cast_to_unstructured_grid()\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.HEX8, material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Cyclic-face node pairing\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pts = np.asarray(grid.points)\nnode_nums = np.asarray(grid.point_data[\"ansys_node_num\"], dtype=np.int64)\nangle = np.arctan2(pts[:, 1], pts[:, 0])\ntol = 1e-6\nlow_mask = np.abs(angle) < tol\nhigh_mask = np.abs(angle - ALPHA) < tol\nlow_pts, high_pts = pts[low_mask], pts[high_mask]\norder_low = np.lexsort(\n    np.column_stack([np.linalg.norm(low_pts[:, :2], axis=1), low_pts[:, 2]]).T[::-1]\n)\norder_high = np.lexsort(\n    np.column_stack([np.linalg.norm(high_pts[:, :2], axis=1), high_pts[:, 2]]).T[::-1]\n)\nlow_paired = np.where(low_mask)[0][order_low]\nhigh_paired = np.where(high_mask)[0][order_high]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Run the harmonic sweep\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cm = femorph_solver.CyclicModel.wrap(\n    m,\n    n_sectors=N_SECTORS,\n    low_face=low_paired,\n    high_face=high_paired,\n)\nresults = cm.solve_modal(n_modes=4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Pick the lowest non-rigid mode per harmonic\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "family: list[tuple[int, float, np.ndarray]] = []\nfor r_k in results:\n    freqs = np.asarray(r_k.frequency, dtype=np.float64)\n    shapes = np.asarray(r_k.mode_shapes)\n    keep = np.where(freqs > 1.0)[0]\n    if keep.size:\n        idx = int(keep[0])\n        phi = shapes[:, idx].reshape(-1, 3)\n        family.append((r_k.harmonic_index, float(freqs[idx]), phi))\n\nprint(f\"recovered {len(family)} non-rigid modes (one per harmonic):\")\nfor k, f, _ in family:\n    print(f\"  k = {k:>2}   f = {f:8.2f} Hz\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render all harmonics in a single 4 \u00d7 2 viewport grid\n\nFor visual consistency, normalise each mode shape to unit\npeak amplitude so the warp factor reads the same regardless\nof the eigsolver's mass-orthonormalisation scaling.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_panels = len(family)\nrows = (n_panels + 1) // 2\nplotter = pv.Plotter(shape=(rows, 2), off_screen=True, window_size=(960, 240 * rows))\n\nfor idx, (k, f, phi) in enumerate(family):\n    row, col = divmod(idx, 2)\n    plotter.subplot(row, col)\n    phi_re = np.real(phi)\n    peak = float(np.max(np.abs(phi_re))) or 1.0\n    warp = grid.copy()\n    warp.points = grid.points + (T_DISK * 4.0 / peak) * phi_re\n    warp[\"UZ\"] = phi_re[:, 2]\n    plotter.add_mesh(grid, style=\"wireframe\", color=\"grey\", opacity=0.4)\n    plotter.add_mesh(warp, scalars=\"UZ\", cmap=\"coolwarm\", show_edges=False, show_scalar_bar=False)\n    plotter.add_text(f\"k = {k}: f = {f:.0f} Hz\", position=\"upper_edge\", font_size=10)\n\n# Hide any unused viewports.\nfor empty_idx in range(n_panels, rows * 2):\n    row, col = divmod(empty_idx, 2)\n    plotter.subplot(row, col)\n    plotter.add_text(\"(unused)\", position=\"upper_edge\", font_size=8, color=\"white\")\n\nplotter.link_views()\nplotter.view_isometric()\nplotter.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.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}