Harmonic (frequency-response) analysis
======================================

:meth:`~femorph_solver.Model.harmonic_solve` and the underlying
:func:`femorph_solver.solvers.harmonic.solve_harmonic_modal` evaluate the
steady-state frequency-response function

.. math::

    \bigl(K \;-\; \omega^{2}\,M \;+\; i\,\omega\,C\bigr)\, u(\omega) \;=\; F

at every excitation frequency in a user-supplied list.  The solve is
**modal-superposition**: the load and the response are projected onto
the lowest ``n_modes`` mass-orthonormal eigenmodes of the undamped
system, and a per-mode complex receptance is back-substituted.

Algorithm
---------

The undamped eigenproblem :math:`K \varphi_{n} = \omega_{n}^{2} M \varphi_{n}`
is solved once (via :doc:`modal`) and the eigenvectors are normalised
so :math:`\Phi^{T} M \Phi = I`.  The modal response coordinate at each
excitation frequency :math:`\omega` is

.. math::

    \eta_{n}(\omega) \;=\;
        \frac{\varphi_{n}^{T}\,F}
             {\omega_{n}^{2} - \omega^{2}
              + 2\,j\,\zeta_{n}\,\omega_{n}\,\omega}

and the physical response is the modal sum
:math:`u(\omega) = \Phi\,\eta(\omega)`.  Every excitation frequency
re-uses the same eigenpairs — the cost of a 200-point sweep is
dominated by the one-time modal solve, not the per-frequency
back-substitution.

API
---

The shortest path lives on :class:`~femorph_solver.Model`.  Build the
model, fix the BCs, stage the load, then call
:meth:`~femorph_solver.Model.harmonic_solve`:

.. code-block:: python

    import numpy as np
    import femorph_solver as fs
    from femorph_solver import ELEMENTS

    m = fs.Model.from_grid(grid)
    m.assign(ELEMENTS.HEX8(integration="enhanced_strain"),
             {"EX": 2.0e11, "PRXY": 0.30, "DENS": 7850.0})
    m.fix(where=clamped_face, dof=["UX", "UY", "UZ"])
    m.apply_force(tip_node, fz=-1.0)

    res = m.harmonic_solve(
        frequencies=np.linspace(0.0, 1000.0, 200),  # Hz
        n_modes=20,
        modal_damping_ratio=0.02,
    )

    res.frequency        # (n_freq,) Hz
    res.displacement     # (n_freq, N) complex
    res.frf_at(dof_idx)  # (n_freq,) complex receptance at one DOF
    res.modal            # the underlying ModalResult that spanned the basis

For low-level work directly on the K / M matrices, the same solver is
exposed as :func:`femorph_solver.solvers.harmonic.solve_harmonic_modal` and
takes the prescribed-DOF dictionary explicitly.

Damping
-------

Two damping models are accepted; both can be combined and their
contributions add.

Modal damping ratio
~~~~~~~~~~~~~~~~~~~

A scalar (uniform across modes) or an array of length ``n_modes``::

    res = m.harmonic_solve(frequencies=freq, n_modes=20,
                           modal_damping_ratio=0.02)

    res = m.harmonic_solve(frequencies=freq, n_modes=4,
                           modal_damping_ratio=[0.01, 0.02, 0.03, 0.05])

Rayleigh damping
~~~~~~~~~~~~~~~~

Pass an ``(α, β)`` pair so the per-mode ratio becomes
:math:`\zeta_{n} = \alpha/(2\,\omega_{n}) + \beta\,\omega_{n}/2`::

    res = m.harmonic_solve(frequencies=freq, n_modes=20,
                           rayleigh=(0.05, 1e-4))

Pure rigid-body modes (:math:`\omega_{n} = 0`) drop the
:math:`\alpha/(2\omega_{n})` term — singular by construction — and
keep whatever ``modal_damping_ratio`` was supplied at that index.

Convergence and basis size
--------------------------

The truncation error at frequency :math:`\omega` is governed by every
mode whose :math:`\omega_{n}` exceeds :math:`\omega`.  In practice:

- Retain at least the first mode whose :math:`f_{n}` exceeds the
  largest excitation frequency in your sweep.
- Add another factor-of-two cushion when the response near a higher
  pole matters (truncation manifests as a soft underestimate of the
  off-resonance amplitude, not a hard miss).
- Static-correction terms are not yet implemented; below the first
  retained pole the truncation residual is bounded but not zero.

Reusing a precomputed modal basis
---------------------------------

Parametric forcing sweeps are cheap when the eigenpairs don't change.
Run the modal solve once, then thread the same :class:`ModalResult`
into successive harmonic calls::

    modal = m.modal_solve(n_modes=20)

    for case in load_cases:
        m.set_load(case.F)
        res = m.harmonic_solve(
            frequencies=case.frequencies,
            modal_result=modal,
            modal_damping_ratio=case.zeta,
        )
        case.dump(res)

Persistence
-----------

:meth:`HarmonicResult.save` writes the complex per-frequency response
to a single ``.pv`` file alongside the underlying frequency grid and
applied damping ratios.  The on-disk format splits each per-frequency
complex array into separate ``displacement_re_NNNN`` and
``displacement_im_NNNN`` point-data entries since pyvista's field
arrays don't carry complex dtype natively; reading the file back via
:func:`femorph_solver.result.read` returns a
:class:`~femorph_solver.result.HarmonicResultDisk` with a
:meth:`~femorph_solver.result.HarmonicResultDisk.displacement` accessor
that recombines the real and imaginary halves.

See also
--------

- :doc:`modal` — the underlying eigenproblem.
- :doc:`transient` — Newmark-β time-domain integration.
- :doc:`choosing-a-solver` — which sparse backend factors the inner
  shift-invert when the modal solve runs at large size.