{
  "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 modal on a rotor sector\n\nA rotor with $N$ identical sectors spanning 360\u00b0 has a spectrum\nthat factors by *harmonic index* $k \\in \\{0, 1, \\dots, N/2\\}$.\nAnalysing one base sector plus the constraint\n\n\\begin{align}u_\\text{high} = e^{i k \\alpha}\\, R(\\alpha)\\, u_\\text{low},\n    \\qquad \\alpha = 2\\pi / N,\\end{align}\n\nrecovers the full-rotor spectrum at a fraction of the DOF count.\n\nThis example builds a simple 4-sector rotor from a pyvista annular\nsector, wraps it in a :class:`~femorph_solver.CyclicModel`, and runs\n:meth:`~femorph_solver.CyclicModel.solve_modal` across all harmonic\nindices.  ``CyclicModel`` auto-detects the low/high cyclic faces from\nthe rotation axis and angle, so the caller never has to assemble a\nmanual DOF pair list.\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 as fs\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 90\u00b0 annular sector\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_SECTORS = 4\nALPHA = 2.0 * np.pi / N_SECTORS\n\nN_T, N_R, N_Z = 12, 6, 3\ntheta = np.linspace(0.0, ALPHA, N_T + 1)\nr = np.linspace(0.8, 1.0, N_R + 1)\nz = np.linspace(0.0, 0.1, N_Z + 1)\n\n# (r, \u03b8, z) order gives a CCW-from-above hex winding in the (r, \u03b8)\n# plane \u2014 swapping r and \u03b8 would flip the winding inside-out.\nR, T, Z = np.meshgrid(r, theta, z, indexing=\"ij\")\nX = R * np.cos(T)\nY = R * np.sin(T)\n\ngrid = pv.StructuredGrid(X, Y, Z).cast_to_unstructured_grid()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Wrap as a CyclicModel\n``CyclicModel.from_grid`` builds the underlying base-sector\n:class:`~femorph_solver.Model`, detects the cyclic face pairs from\nthe rotation axis (default ``\"z\"``) + ``n_sectors``, and validates\nthat every low-face point has a rotated counterpart on the high\nface within ``pair_tol``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cm = fs.CyclicModel.from_grid(grid, n_sectors=N_SECTORS)\n\n# Per-sector physics \u2014 element kernel + isotropic material \u2014 is set\n# on the inner Model exactly as for a non-cyclic analysis.\ncm.assign(ELEMENTS.HEX8, material={\"EX\": 2.1e11, \"PRXY\": 0.30, \"DENS\": 7850.0})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve every harmonic index\n``modal_solve`` defaults to all harmonics ``k = 0 \u2026 N/2``, returning\none :class:`~femorph_solver.solvers.cyclic.CyclicModalResult` per\nharmonic.  Each result carries ``frequency`` and ``mode_shapes``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "results = cm.solve_modal(n_modes=4)\n\nprint(f\"Cyclic modal on {N_SECTORS}-sector rotor:\")\nfor res in results:\n    fs_str = \", \".join(f\"{f:.2f}\" for f in res.frequency[:4])\n    print(f\"  k = {res.harmonic_index}: f [Hz] = [{fs_str}]\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualise mode 0 of k = 1\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "r1 = next(res for res in results if res.harmonic_index == 1)\nphi = r1.mode_shapes[:, 0].reshape(-1, 3)\ngrid.point_data[\"mode_k1\"] = np.abs(phi)\n\nplotter = pv.Plotter(off_screen=True)\nplotter.add_mesh(\n    grid.warp_by_vector(\"mode_k1\", factor=0.05 / (np.max(np.abs(phi)) + 1e-30)),\n    scalars=np.linalg.norm(phi, axis=1),\n    cmap=\"viridis\",\n    scalar_bar_args={\"title\": f\"|mode| (k=1, f={r1.frequency[0]:.1f} Hz)\"},\n)\nplotter.add_axes()\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
}