SMI-Lab-Inha · Sep 30, 2026 · Sep 29, 2026 · Sep 29, 2026 · Sep 29, 2026 · Sep 30, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -8,7 +8,109 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/).
 
 ## [Unreleased]
 
-(nothing yet)
+### Fixed
+
+- **Dense symmetric solves depended on how many modes were requested.**
+  The dense path reduced `K x = λ M x` through a Cholesky factor of the
+  mass matrix, and a backward-stable symmetric eigensolver gets each
+  eigenvalue of the reduced matrix only to an *absolute* accuracy of
+  about `eps · λ_max`. The ElastoDyn adapters floor rotary inertia at
+  1e-6, which puts `λ_max / λ_min` near 1e15, so the lowest modes carried
+  errors of several percent that moved with the `subset_by_index`
+  window. On the bundled NREL 5MW blade at 12.1 rpm the first flap mode
+  read **0.7154 / 0.7096 / 0.7661 / 0.7282 Hz** for `n_modes` = 4 / 6 /
+  10 / 20, and first edge 1.1221 / 1.1184 / 1.0884 / 1.1303 Hz; the tower
+  fore-aft / side-side pair moved and swapped. Every model below the
+  500-DOF sparse threshold took this path — every ElastoDyn blade and
+  tower.
+
+  The dense path now solves the inverted pencil `M x = μ (K + s M) x`,
+  `μ = 1 / (λ + s)`, which puts the wanted modes at the top of the
+  spectrum where the same absolute accuracy is a relative one — the
+  shift-invert transformation the sparse path already used. The shift is
+  zero when `K` is positive definite and is otherwise grown from an
+  estimate of the lowest eigenvalues until `K + s M` factorises, so
+  rigid-body modes and the negative eigenvalues of a column loaded past
+  buckling are still returned. A request reaching into the upper half of
+  the spectrum takes those modes from the mass-reduced form, which is the
+  accurate one there. Results are now independent of `n_modes` and agree
+  with the sparse path to roundoff: blade flap1 / edge1 / flap2 =
+  **0.72865 / 1.09754 / 2.00842 Hz** at every window, tower 1st pair
+  0.33483 / 0.33561 Hz, backward error ~1e-10 (was ~5e-2). Frequencies of
+  every dense-path model change by up to the old error — on the NREL 5MW
+  blade up to ~5 % on the lowest flap / edge modes, on the tower ~0.3 %.
+
+  When one strongly negative eigenvalue forces a shift much larger than
+  the soft end of the spectrum, those soft modes all map to the same `μ`
+  to working precision. The window is then grown until its edge falls in
+  a resolvable gap, and each such cluster is re-solved by Rayleigh-Ritz
+  on its own subspace, so the soft modes come back in order and
+  independent of `n_modes` (`K = diag(−1e16, 1, 2, 3, 4)` returns
+  `[−1e16, 1]` for two modes, not `[−1e16, 4]`).
+
+  The mass-reduced form is kept as the fallback for a pencil no shift
+  makes definite, and the residual retry of 1.18.0 still stands behind
+  it; the near-singular-mass case it was built for now solves directly
+  on the dense path without a retry or a warning.
+  `tests/test_campbell.py::test_first_flap_lifts_endpoint_to_endpoint`,
+  which sat at 2.95 % against its 3 % gate on some BLAS builds, now
+  measures 7.6 %.
+- **Floating models built from OpenFAST decks lost their roll and pitch
+  modes.** `Tower.from_elastodyn_with_mooring` took the restoring from
+  the WAMIT `.hst`, which holds only the hydrostatic part
+  `ρ g I_wp + ρ g V z_B`. OpenFAST adds the body weight through
+  ElastoDyn's gravity; the model did not, so a ballast-stabilised spar
+  read as unstable — OC3 Hywind `C44 = C55 = −5.0e9 N·m/rad` — and the
+  general eigensolver silently dropped the two negative-eigenvalue modes,
+  returning four of the six rigid-body modes. The constructor now adds the
+  weight restoring `−g Σ mᵢ (zᵢ − z_ref)` of the platform, the tower and
+  the RNA to `hydro_K` (+6.17e9 N·m/rad on OC3: +6.58e9 of ballast,
+  −0.41e9 of tower and RNA). It is added only when a HydroDyn deck is
+  given, since weight without the buoyancy that balances it does not
+  describe a floating body. `Tower.from_windio_floating` adds the same
+  term in both of its tiers, so the deck-backed tier stays equivalent to
+  `from_elastodyn_with_mooring`. OC3 Hywind from the r-test decks now
+  gives surge / sway / heave / roll / pitch / yaw = 0.00806 / 0.00806 /
+  0.03241 / **0.03383 / 0.03383** / 0.04121 Hz against Jonkman (2010)
+  0.0080 / 0.0080 / 0.0324 / 0.0342 / 0.0343 Hz (roll and pitch were
+  absent before). Yaw is the catenary-only value; the published 0.1210 Hz
+  includes an additional 9.834e7 N·m/rad crowfoot spring the r-test
+  MoorDyn deck does not carry. The coupled first tower pair moves from
+  0.4887 / 0.4903 to 0.4922 / 0.4938 Hz.
+- **A negative-stiffness rigid-body mode is no longer dropped silently.**
+  `Tower.run` on a free-free floating model now warns (`UserWarning`)
+  when `hydro_K + mooring_K` has a negative direction, naming the
+  dominant DOFs: the platform is statically unstable there and the
+  general solver leaves those modes out of the result. The bundled
+  `08_nrel5mw_oc4semi` sample triggers it — the OC4 DeepCwind deck keeps
+  its ballast in HydroDyn filled members rather than `PtfmMass`, which
+  the platform model does not include, so its roll / pitch restoring is
+  negative and the modes labelled roll / pitch there are tower bending.
+- **`Tower.add_point_mass` on a deck-built floater now updates the
+  weight restoring.** The weight term is formed when the model is
+  constructed, so a lump added afterwards reached the FEM mass matrix but
+  not the roll / pitch restoring. The lump's `−m g (z − z_ref)` is now
+  added to `hydro_K` as well, which matches a model rebuilt with the lump
+  in place.
+- **NumPy 1.26 support.** Several code paths called `np.trapezoid`,
+  which only exists from NumPy 2.0, so on the supported `numpy>=1.26`
+  every `Tower.run()` (through the pre-solve checks), every ElastoDyn
+  blade read and every floating constructor failed with
+  `AttributeError`. They now share one version-independent trapezoid
+  helper, and a test keeps the NumPy-2-only spelling out of the package.
+- **The first tower pair of the IEA-3.4 and IFE UPSCALE 25MW decks is no
+  longer degenerate.** Both towers are symmetric in stiffness, but the
+  lumped RNA is not (pitch inertia 1.9× and 3.0× roll inertia, plus an
+  `ixz` product), which splits the first FA / SS pair by 0.09 % and 1.0 %.
+  The dense solver above, the sparse path and an `eigh(M, K)` solve all
+  agree on the split. The old dense path returned the pair degenerate
+  only because its error on the lowest modes was larger than the split,
+  and two integration tests had pinned that; they now pin the split.
+- **WindIO discovery ignored every yaml under a path containing
+  "OpenFAST".** The exclusion meant for the `OpenFAST/` deck directory of
+  an RWT layout matched the absolute path, so a checkout living under
+  e.g. `~/OpenFAST-GUI/` found no ontology at all. Only the part of the
+  path below the search root is examined now.
 
 ## [1.18.0] — 2026-08-13
 

diff --git a/VALIDATION.md b/VALIDATION.md
@@ -102,7 +102,9 @@ metrics:
 | Discrete mid-span point mass (issue #35) | cantilever carrying one lump at station $a$: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1) | 1st frequency; mesh-position independence | < 0.5 %; coarse-vs-fine < 0.2 % | (within tol) | [`tests/fem/test_gravity_and_point_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_gravity_and_point_mass.py) | no |
 | Distributed Winkler soil bed vs the lumped mudline condensation (issue #118) | Psaroudakis et al. (2021) / Yu & Amdahl (2023) Eq. 25 is the exact static condensation of a constant-EI pile on a bed of rate $k = D_P E_{SO}$ | coupled 1st frequency, distributed bed vs lumped springs | < 1 % | 0.4 % (the embedded pile inertia the condensed form drops) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no |
 | Distributed Winkler bed converges on the rigid mudline clamp (issue #118) | limit $E_{SO} \to \infty$; residual compliance scales as the elastic length $(4EI/k)^{1/4}$ | 1st frequency at $E_{SO} \times 10^8$ vs the clamped model | < 1 % | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no |
-| Near-singular mass matrix does not silently corrupt the low modes | cantilever carrying one lump: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1), at a 4e5:1 lump-to-beam mass ratio | 1st frequency, dense and sparse dispatch sizes | < 0.5 % | (within tol; the unguarded symmetric solve is 137 % out) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no |
+| Dense symmetric solve is independent of the requested mode count | NREL 5MW ElastoDyn blade (12.1 rpm) and tower, bundled deck; shift-invert spectral transformation (Ericsson & Ruhe 1980) | lowest modes for 4 vs 20 requested; dense vs sparse shift-invert | `rtol = 1e-6` | identical to roundoff (blade flap1 0.72865 Hz at every window; was 0.715 / 0.710 / 0.766 / 0.728 Hz for 4 / 6 / 10 / 20) | [`tests/fem/test_dense_inverted_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_dense_inverted_solver.py) | no |
+| Near-singular mass solved directly on the dense path | cantilever carrying one lump: $f = \sqrt{3EI/(m L^3)}/2\pi$ (Blevins 1979, Table 8-1), 4e5:1 lump-to-beam mass ratio | 1st frequency, no residual retry, backward error | < 0.5 %; residual < 1e-6 | (within tol) | [`tests/fem/test_dense_inverted_solver.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_dense_inverted_solver.py) | no |
+| Near-singular mass matrix does not silently corrupt the low modes | cantilever carrying one lump: $f = \sqrt{3EI/(m a^3)}/2\pi$ (Blevins 1979, Table 8-1), at a 4e5:1 lump-to-beam mass ratio | 1st frequency, dense and sparse dispatch sizes, through the residual retry behind the mass-reduced fallback | < 0.5 % | (within tol; the unguarded mass-reduced solve is 137 % out) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no |
 | Rigid-body modes are not lost to the solver guard | construction (rank-deficient $K$, well-conditioned $M$; requested subsets of 1, 3, 6 and 10 modes) | zero modes retained, requested mode count returned | exact | (within tol) | [`tests/fem/test_ill_conditioned_mass.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/fem/test_ill_conditioned_mass.py) | no |
 | Deck-reader `n_nodes` refinement (issue #58) | Euler-Bernoulli closed form, modes 1-3, plus self-convergence $n{=}100$ vs $200$ | bending frequencies | < 1 %; self-convergence < 0.2 % | (within tol) | [`tests/test_refine_mesh.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_refine_mesh.py) | no |
 
@@ -114,6 +116,8 @@ by title + publisher:
 - Blevins (1979). *Formulas for Natural Frequency and Mode Shape*. Krieger Publishing. (Textbook; no DOI.)
 - Karnovsky & Lebed (2001). *Formulas for Structural Dynamics*. McGraw-Hill. (Textbook; no DOI.)
 - Timoshenko & Gere (1961). *Theory of Elastic Stability*, 2nd ed. McGraw-Hill. (Textbook; no DOI.)
+- Ericsson & Ruhe (1980). *The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems*. *Mathematics of Computation* 35(152), 1251–1268. DOI [10.1090/S0025-5718-1980-0583502-2](https://doi.org/10.1090/S0025-5718-1980-0583502-2).
+- Faltinsen (1990). *Sea Loads on Ships and Offshore Structures*. Cambridge University Press. (Textbook; no DOI.)
 - Greenhill (1881). *Determination of the greatest height consistent with stability that a vertical pole or mast can be made*. *Proceedings of the Cambridge Philosophical Society* 4, 65–73. (No DOI; tabulated in Timoshenko & Gere 1961, Art. 2.16.)
 - Shadlou & Bhattacharya (2016). *Dynamic stiffness of monopiles supporting offshore wind turbine generators*. *Soil Dynamics and Earthquake Engineering* 88, 15–32. DOI [10.1016/j.soildyn.2016.04.002](https://doi.org/10.1016/j.soildyn.2016.04.002).
 - Yu & Amdahl (2023). *A Rayleigh-Ritz solution for high order natural frequencies and eigenmodes of monopile supported offshore wind turbines considering tapered towers and soil-pile interactions*. *Marine Structures* 92, 103482. DOI [10.1016/j.marstruc.2023.103482](https://doi.org/10.1016/j.marstruc.2023.103482).
@@ -223,6 +227,9 @@ is behavioural / contract-style.
 | `MooringSystem.stiffness_matrix` — OC3 yaw stiffness (catenary-only) | OC3 r-test MoorDyn `.dat` | `K[5,5] ≈ 1.156e7 N·m/rad` (catenary lines only; OC3's delta-line crowfoot adds the bulk of the published 9.83e7) | `rel = 5%` | (within tol) | [`tests/test_mooring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_mooring.py) | yes |
 | `MooringSystem.stiffness_matrix` — 3-fold-symmetry signatures | OC3 r-test MoorDyn `.dat` | K[0,0] = K[1,1], K[3,3] = K[4,4], K[0,4]·K[1,3] < 0 | `rel = 1e-4` | (within tol) | [`tests/test_mooring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_mooring.py) | yes |
 | `Tower.from_elastodyn_with_mooring` — end-to-end OC3 coupled solve | OC3 ElastoDyn + r-test MoorDyn `.dat` | hub\_conn=2, PlatformSupport populated, 1st tower-bending in 0.4–0.6 Hz | `rel = 5%` | within 1.2 % of Jonkman 2010's 0.482 Hz | [`tests/test_mooring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_mooring.py) | yes |
+| Floating weight restoring $-g \sum m_i (z_i - z_{ref})$ in roll / pitch (platform + tower + RNA) | closed form; classical $C_{44} = \rho g V (z_B - z_G) + \rho g I_{wp}$ (Faltinsen 1990, ch. 2) | 6×6 weight matrix, reference-point shift, mass scaling, point masses | `rel = 1e-9` | exact | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | no |
+| Negative platform restoring is reported, not silently dropped | construction (OC3 Hywind sample with the weight term removed) | `UserWarning` naming roll / pitch | always | (within tol) | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | no |
+| **`Tower.from_elastodyn_with_mooring` — OC3 Hywind rigid-body frequencies** (WAMIT `.hst` + weight of platform, tower and RNA) | Jonkman (2010) NREL/TP-500-47535, FAST full-system natural frequencies: surge/sway 0.0080, heave 0.0324, roll 0.0342, pitch 0.0343 Hz | six labelled rigid-body modes | surge/sway/roll/pitch 3 %, heave 2 % | surge/sway 0.8 %, heave 0.0 %, roll 1.1 %, pitch 1.4 % (roll/pitch were missing before the weight term) | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | yes |
 | `MudlineFoundation.from_soil_properties` — Yu & Amdahl (2023) Table 9 DTU 10 MW flexible pile, E\_SO = 30 MPa | Yu and Amdahl (2023) Marine Structures 92, 103482, Table 9 row 2 | K\_hh, K\_hr, K\_rr via Shadlou and Bhattacharya (2016) formulas | `rel = 10%` | within 3.0 % on all three springs | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no |
 | `MudlineFoundation.from_soil_properties` — Yu & Amdahl (2023) Table 9 DTU 10 MW rigid pile, E\_SO = 30 MPa | Yu and Amdahl (2023) Table 9 row 1 | K\_hh, K\_hr, K\_rr via Shadlou and Bhattacharya (2016) rigid-pile formulas | `rel = 5%` | within 2.9 % on all three springs | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no |
 | `MudlineFoundation.as_mooring_K` — cross-coupling sign matches Jonkman 2010 OC3 convention | construction | K\[0,4\] < 0, K\[1,3\] > 0, K\[0,4\]·K\[1,3\] < 0 | always | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no |

diff --git a/docs/conventions.rst b/docs/conventions.rst
@@ -289,6 +289,17 @@ Common pitfalls
 * **Leaving ``mooring_K`` in the platform/hull frame.** Like hydro, it
   must be referenced to the tower axis.
 
+* **Passing a WAMIT ``.hst`` as the whole restoring.** The ``.hst`` is
+  hydrostatic only (``ρ g I_wp + ρ g V z_B``); the body weight's
+  ``−m g z_G`` in roll / pitch is applied by OpenFAST through ElastoDyn
+  gravity. ``hydro_K + mooring_K`` must carry both. The deck
+  constructors add the weight of platform, tower and RNA — each mass
+  at its height above the ``ref_msl`` reference — for you; a
+  hand-built ``PlatformSupport`` has to include it itself. Without it
+  a ballast-stabilised spar is negative in roll / pitch (OC3 Hywind:
+  ``−5.0e9 N·m/rad``), and ``Tower.run`` warns that those rigid-body
+  modes are unstable rather than returning them.
+
 * **Confusing OpenFAST's two frames.** ``PtfmCMzt`` (CM) and
   ``PtfmRefzt`` (reference) are both in the tower-base *t*-frame;
   pyBmodes' ``ref_msl`` is ``PtfmRefzt``. Don't add one to the other.

diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/build.py b/src/pybmodes/_examples/sample_inputs/reference_turbines/build.py
@@ -89,6 +89,7 @@
 if SRC_DIR.is_dir() and str(SRC_DIR) not in sys.path:
     sys.path.insert(0, str(SRC_DIR))
 
+from pybmodes._numeric import trapezoid  # noqa: E402
 from pybmodes.io.bmi import PlatformSupport  # noqa: E402
 from pybmodes.io.elastodyn_reader import (  # noqa: E402
     read_elastodyn_blade,
@@ -172,7 +173,7 @@ def _blade_total_mass(blade) -> float:
     span = np.asarray(blade.bl_fract, dtype=float)
     rho = np.asarray(blade.b_mass_den, dtype=float)
     # bl_fract is normalized [0, 1]; multiply by blade length elsewhere.
-    return float(np.trapezoid(rho, span))
+    return trapezoid(rho, span)
 
 
 def _emit_section_properties_table(

diff --git a/src/pybmodes/_numeric.py b/src/pybmodes/_numeric.py
@@ -0,0 +1,35 @@
+# Copyright 2024-2026 Jae Hoon Seo
+# Marine Structural Mechanics and Integrity Lab (SMI Lab), Inha University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Small numerical helpers shared across the package."""
+
+from __future__ import annotations
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+__all__ = ["trapezoid"]
+
+
+def trapezoid(y: ArrayLike, x: ArrayLike) -> float:
+    """Trapezoidal integral of ``y`` over ``x``.
+
+    Implemented directly rather than through ``np.trapezoid``, which only
+    exists from NumPy 2.0, or ``np.trapz``, which NumPy 2 deprecates, so
+    the package works across the whole advertised ``numpy>=1.26`` range.
+    """
+    y = np.asarray(y, dtype=float)
+    x = np.asarray(x, dtype=float)
+    return float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(x)))
diff --git a/src/pybmodes/checks.py b/src/pybmodes/checks.py
@@ -102,6 +102,7 @@
 
 import numpy as np
 
+from pybmodes._numeric import trapezoid
 from pybmodes.options import DEFAULT_CHECK_OPTIONS as _CHECK_OPTIONS
 
 if TYPE_CHECKING:
@@ -269,7 +270,7 @@ def _check_section_properties_finite(
 ) -> None:
     """Flag any non-finite (NaN / ±Inf) entry in the numeric section-
     property fields. ERROR-severity because every downstream
-    consumer (``np.trapezoid``, ``np.linalg.eigh``, the FE assembly)
+    consumer (the trapezoid mass integral, ``np.linalg.eigh``, the FE assembly)
     silently produces NaN-filled outputs on NaN inputs.
 
     The per-field checks below this one (``_check_span_monotonic`` /
@@ -395,7 +396,7 @@ def _check_rna_vs_tower_mass(
     mass_den = np.asarray(sp.mass_den, dtype=float)
     if span_phys.size < 2 or not np.all(np.diff(span_phys) > 0):
         return  # other checks will flag the bad span; skip cleanly here
-    tower_mass = float(np.trapezoid(mass_den, span_phys))
+    tower_mass = trapezoid(mass_den, span_phys)
     if tower_mass <= 0.0:
         return  # mass-density check will catch this
     rna = float(bmi.tip_mass.mass)