Cyclic-symmetry modal analysis
==============================

A rotor built from ``N`` identical sectors spanning 360° admits a
*cyclic-symmetry* modal reduction: every full-rotor eigenmode satisfies

.. math::

   u_{\text{sector } j+1}(\text{global}) = e^{i k \alpha}\, R(\alpha)\,
       u_{\text{sector } j}(\text{global})

with :math:`\alpha = 2\pi / N` and :math:`R(\alpha)` the spatial
rotation about the symmetry axis.  The integer
:math:`k \in \{0, 1, \dots, N/2\}` is the *harmonic index* (or nodal
diameter).  Solving one sector per ``k`` reproduces the full-rotor
spectrum at a cost of :math:`1 / N` per harmonic — for a 15-sector
bladed rotor the eight needed sector solves already beat a single
full-rotor eigensolve.

.. contents::
   :local:
   :depth: 2


High-level API
--------------

The cyclic-aware front door is :class:`femorph_solver.CyclicModel`:
it wraps a base-sector :class:`femorph_solver.Model` plus the cyclic
parameters (``n_sectors`` and the rotation axis) and detects the
low / high cyclic-face node pairing internally.  Downstream code does
not need to roll its own KD-tree or hand-build DOF index arrays.

.. code-block:: python

   import femorph_solver as fs

   cm = fs.CyclicModel.from_pv(
       fs.examples.cyclic_bladed_rotor_sector_path(),
       n_sectors=15,
       axis="z",       # default; pass "x"/"y"/"z" or a 3-vector
   )
   results = cm.modal_solve(n_modes=4)
   for r in results:
       print(f"k={r.harmonic_index}  f={r.frequency}")

Each entry of ``results`` is a
:class:`~femorph_solver.solvers.cyclic.CyclicModalResult` for one harmonic
index, holding ``omega_sq``, ``frequency``, and the complex
base-sector ``mode_shapes``.  When you want the full-rotor frequency
multiset (each intermediate-``k`` mode counted twice), call
``cm.aggregated_frequencies(n_modes=4)``.


Manual face override
--------------------

When auto-detection misclassifies the cyclic faces — for example, on a
mesh whose interface nodes don't sit within ``pair_tol`` of their
rotated partners because of upstream meshing tolerance — pass the
low / high faces explicitly as grid-point-index arrays:

.. code-block:: python

   cm = fs.CyclicModel.from_pv(
       fs.examples.cyclic_bladed_rotor_sector_path(),
       n_sectors=15,
       low_face=low_indices,
       high_face=high_indices,
   )

The arrays must have the same length and pass an ``R(α)``-image-within-
``pair_tol`` consistency check; misuses raise ``ValueError`` with a
clear message naming the offending pair.


Loading a Model from a ``.pv`` file
-----------------------------------

The bundled rotor sector ships as a single native ``.pv`` file:
geometry plus element-type assignments, material properties, and
unit-system stamp baked in.  This is the same format
:meth:`~femorph_solver.Model.save` writes, so any ``Model`` you
serialise can come back through :meth:`~femorph_solver.Model.from_pv`
ready-to-solve, with no APDL setup:

.. code-block:: python

   import femorph_solver as fs

   m = fs.Model.from_pv("rotor_sector.pv")
   r = m.modal_solve(n_modes=4)         # full-sector modal (no cyclic BC)

   cm = fs.CyclicModel.from_pv("rotor_sector.pv", n_sectors=15)
   r = cm.modal_solve(n_modes=4)        # cyclic sweep across harmonics


Underneath: the building blocks
-------------------------------

:func:`femorph_solver.solvers.cyclic.solve_cyclic_modal` is the lower-
level entry point if you have a model assembled outside this package
or want to drive the cyclic reduction directly:

.. code-block:: python

   from femorph_solver.solvers.cyclic import solve_cyclic_modal

   results = solve_cyclic_modal(
       K,                    # (N, N) sector stiffness, real SPD
       M,                    # (N, N) sector mass, real SPD
       low_node_dofs,        # (P, d) global DOFs at P paired low-face nodes
       high_node_dofs,       # (P, d) paired high-face DOFs (same labels, same order)
       n_sectors=N,
       n_modes=10,
       pair_rotation=R_alpha,
       harmonic_indices=None,
   )

The constraint matrix is assembled directly in COO form: low-face and
high-face DOFs are identified, the high-face DOFs are eliminated via
the phase-rotation relation, and the reduced Hermitian generalised
problem is solved with :func:`scipy.linalg.eigh` (dense, subset by
index).  For ``k = 0`` and ``k = N/2`` (even ``N``) the phase is
:math:`\pm 1` and the reduced system is real; intermediate ``k`` give
a complex Hermitian reduction.

:class:`CyclicModel` does this assembly + DOF lookup for you: the
inner :class:`~femorph_solver.Model` exposes ``stiffness_matrix()``
and ``mass_matrix()``, the wrapper translates its cached point-pair
arrays into global DOF indices, and the spatial rotation is built
from the user-supplied axis / angle.  Call ``cm.modal_solve(...)``
and you get the same ``CyclicModalResult`` list back.


Harmonic-index counting
-----------------------

The cyclic sweep is equivalent to the full-rotor eigenproblem — every
full-rotor frequency falls in exactly one harmonic index — but modes
are *degenerate* in pairs for intermediate ``k``:

- ``k = 0``: single real eigenmodes (the "in-phase" family)
- ``k = N/2`` (even ``N`` only): single real eigenmodes (anti-phase)
- ``0 < k < N/2``: each eigenvalue represents a travelling-wave pair in
  the full rotor; count it **twice** when aggregating to the full-rotor
  multiset.

Rigid-body modes are partitioned across ``k = 0`` (axial translation
+ axial rotation) and ``k = 1`` (the in-plane translation pair).


Why the spatial rotation matters
--------------------------------

Stiffness / mass are usually assembled in a global XYZ frame, so the
phase relation between neighbouring sectors picks up not just the
scalar ``exp(i k α)`` but also the rotation :math:`R(\alpha)` that
takes sector ``j`` to sector ``j+1``.  :class:`CyclicModel` builds it
from the ``axis`` / ``angle`` kwargs (default: ``α = 2π / n_sectors``
about ``+z``) and threads it through the lower-level solver
automatically; only when DOFs are already in a per-sector cylindrical
frame would you bypass :class:`CyclicModel` and pass
``pair_rotation=None`` to :func:`solve_cyclic_modal` directly.

Elements with both translational and rotational DOFs (shells, beams)
use a block-diagonal ``R`` of two copies of the spatial rotation — one
for the translations, one for the rotations.  The lower-level solver
itself is element-agnostic: it reads only the ``(P, d)`` DOF-index
arrays and the ``(d, d)`` rotation.


Bundled example — bladed-rotor sector
-------------------------------------

The femorph-solver wheel ships a 15-sector × 24° bladed-rotor base
sector (101 ``HEX8`` cells) at
``femorph_solver.examples.cyclic_bladed_rotor_sector_path()``.  The
file is a complete :meth:`~femorph_solver.Model.save` snapshot: the
HEX8 element-type and isotropic material assignments are baked in,
so the demo above runs without any external download or APDL setup.

Cyclic faces on this mesh are not ``θ=const`` planes — they're curved
surfaces that map to each other under a 24° rotation about z.
:class:`CyclicModel` resolves the pairing uniformly on the rotated
point-cloud match.

The deeper-history version of this example (driven directly off
``mapdl_archive.examples.sector_archive_file`` with hand-rolled
``cKDTree`` pair detection) is exercised in
``tests/elements/solid185/test_solid185_cyclic_sector_archive.py`` as
a regression against the legacy code path.