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


API
---

:func:`femorph_solver.solvers.cyclic.solve_cyclic_modal` takes the base-sector
stiffness / mass plus a one-to-one DOF pairing between the low-angle
and high-angle cyclic faces:

.. 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,  # (d, d) spatial rotation, or None for local frames
       harmonic_indices=None,  # defaults to [0, 1, …, N//2]
   )

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``.

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.


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``.  Pass it as the ``(d, d)``
matrix acting on each paired node's ``d`` translational DOFs (3 for a
3D solid sweeping about z).  Omit it (``None``) only when face DOFs
are already in a per-sector cylindrical frame, where :math:`R(\alpha)`
reduces to the identity.

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 solver itself is
element-agnostic: it reads only the ``(P, d)`` DOF-index arrays and
the ``(d, d)`` rotation.


Example — bladed-rotor sector
-----------------------------

The femorph-solver test suite ships the 15-sector bladed rotor sector from
``mapdl_archive.examples.sector_archive_file`` (101 SOLID185 hex
elements per sector, ``α = 24°``) along with a port of femorph's
rotation-matching pair detection.  The end-to-end flow:

.. code-block:: python

   import mapdl_archive
   import numpy as np
   from scipy.spatial import cKDTree

   import femorph_solver
   from femorph_solver.solvers.cyclic import solve_cyclic_modal

   a = mapdl_archive.Archive(mapdl_archive.examples.sector_archive_file)
   grid = a.grid.extract_cells(np.where(a.grid.celltypes == 12)[0]).cast_to_unstructured_grid()

   # Auto-detect N and low/high node pairing.
   pts = np.asarray(grid.points)
   tree = cKDTree(pts)
   N = 15                                  # or scan N via rotation+match
   alpha = 2 * np.pi / N
   R = np.array([[np.cos(alpha), -np.sin(alpha), 0.0],
                 [np.sin(alpha),  np.cos(alpha), 0.0],
                 [0.0,            0.0,           1.0]])
   d, idx = tree.query(pts @ R.T, k=1)
   low_pt = np.where(d < 1e-4)[0]
   high_pt = idx[low_pt]

   # Build the sector model and look up UX/UY/UZ DOF indices per pair.
   m = femorph_solver.Model.from_grid(grid)
   m.et(185, "SOLID185")
   m.mp("EX", 1, 16.9e6);  m.mp("PRXY", 1, 0.31);  m.mp("DENS", 1, 4.1408e-4)

   dof_map = m.dof_map()
   def node_dofs(nn):
       nn_col = dof_map[:, 0]
       return np.array([int(np.where((nn_col == n) & (dof_map[:, 1] == j))[0][0])
                        for n in nn for j in (0, 1, 2)]).reshape(-1, 3)

   nns = np.asarray(grid.point_data["ansys_node_num"])
   low_nd = node_dofs(nns[low_pt])
   high_nd = node_dofs(nns[high_pt])

   results = solve_cyclic_modal(
       m.stiffness_matrix(), m.mass_matrix(),
       low_nd, high_nd,
       n_sectors=N, n_modes=4,
       pair_rotation=R,
   )
   for r in results:
       print(f"k={r.harmonic_index}  f={r.frequency}")

That full example is exercised in
``tests/elements/solid185/test_solid185_cyclic_sector_archive.py``,
which verifies the sector sweep reproduces the full 15-sector rotor's
spectrum to within 0.5 % on the first 10 elastic modes.