2323mass / moment-of-inertia metric that makes a translation amplitude
2424comparable to a rotation amplitude).
2525
26- The classifier is deliberately *conservative*: it names a mode only
27- when one platform DOF clearly dominates the base-node modal kinetic
28- energy. A flexible tower mode, or a strongly coupled / eigensolver-
29- rotated near-degenerate pair (surge≈sway, roll≈pitch), is left
30- ``None`` rather than mislabelled — consistent with the project's
31- "only name what's unambiguous" stance. Empirically every genuine
32- rigid-body mode on the validated floating decks (OC3 Hywind, the
33- IEA-15 / IEA-22 / OC4 / UPSCALE samples) is overwhelmingly single-DOF
34- and is named; the first flexible bending pair and above stay ``None``.
26+ The DOF labels are assigned *globally*: the six rigid-body candidates
27+ and the six platform DOFs are matched by a Hungarian optimal
28+ assignment over the per-DOF modal-kinetic-energy fractions, so the
29+ mode that best expresses each DOF wins that label. The earlier greedy
30+ "argmax per mode, drop later duplicates" rule mislabelled asymmetric
31+ platforms — an early mode whose mass-weighted energy was inflated by a
32+ parasitic high-inertia DOF (a small coupled yaw rotation reads as huge
33+ energy) would steal that DOF's label, starve the true owner into
34+ ``None``, and cascade onto the remaining modes (issue #93). The global
35+ assignment removes that order/greed sensitivity; on symmetric decks,
36+ where every rigid-body mode is ~98 % single-DOF, it reduces to the
37+ same per-DOF matching as before.
38+
39+ Before the assignment, frequency-degenerate rigid pairs are rotated
40+ onto axis-aligned directions. On a (bi)symmetric platform surge≈sway
41+ and roll≈pitch share an eigenvalue, so the eigensolver may return any
42+ rotation of that 2-D eigenspace — a non-deterministic basis (BLAS
43+ thread-order dependent) that would otherwise leave a 45°-mixed pair
44+ below the dominance threshold and unnamed. Rotating each degenerate
45+ pair onto its platform axes (the rigid-body analog of the FA/SS
46+ resolver in ``pybmodes.elastodyn.params``) makes the labelling
47+ deterministic. Asymmetric platforms break the degeneracy, so this step
48+ is a no-op there and the global assignment carries the work.
49+
50+ The classifier stays deliberately *conservative*: after assignment a
51+ mode is named only when its assigned platform DOF clearly dominates
52+ the base-node modal kinetic energy (``_DOMINANCE_THRESHOLD``). A
53+ flexible tower mode, or a genuinely strongly-coupled pair whose energy
54+ splits across DOFs, is left ``None`` rather than mislabelled —
55+ consistent with the project's "only name what's unambiguous" stance.
56+ Empirically every genuine rigid-body mode on the validated floating
57+ decks (OC3 Hywind, the IEA-15 / IEA-22 / OC4 / UPSCALE samples) is
58+ overwhelmingly single-DOF and is named; the first flexible bending
59+ pair and above stay ``None``.
3560"""
3661
3762from __future__ import annotations
6994# the base and would otherwise be mis-named.
7095_N_RIGID = 6
7196
97+ # Two rigid-body modes whose relative frequency gap is below this are
98+ # treated as a degenerate pair: on a (bi)symmetric platform surge≈sway
99+ # and roll≈pitch share an eigenvalue, so the eigensolver returns an
100+ # arbitrary rotation of that 2-D eigenspace (non-deterministic across
101+ # BLAS thread orderings — the same hazard commit-fixed for the FA/SS
102+ # tower pair). Matches ``_DEGENERATE_FREQ_RTOL`` in
103+ # ``pybmodes.elastodyn.params``.
104+ _DEGENERATE_FREQ_RTOL = 1e-4
105+
106+
107+ def _platform_dof_energy (b : np .ndarray , Mp : np .ndarray ) -> np .ndarray :
108+ """Per-platform-DOF modal kinetic energy of a base 6-vector ``b``
109+ (FEM order), reordered to platform DOF order. Not normalised.
110+
111+ The platform mass matrix is the metric that puts a translation
112+ amplitude and a rotation amplitude on a comparable footing
113+ (``b_i·(M_p b)_i``).
114+ """
115+ e = np .abs (b * (Mp @ b ))
116+ return np .asarray (e [_FEM_TO_PLATFORM ], dtype = float )
117+
118+
119+ def _align_degenerate_rigid_pairs (
120+ base : np .ndarray , freqs : np .ndarray , Mp : np .ndarray
121+ ) -> np .ndarray :
122+ """Rotate near-degenerate rigid-body base-motion pairs onto axis-
123+ aligned directions.
124+
125+ ``base`` is ``(6, n_rigid)`` of tower-base 6-vectors (FEM order),
126+ ``freqs`` the matching frequencies. Walks consecutive pairs; when a
127+ pair is frequency-degenerate (``_DEGENERATE_FREQ_RTOL``) it rotates
128+ the 2-D eigenspace — any in-plane rotation is an equally valid
129+ eigenpair — to put the pair's dominant platform DOF entirely in the
130+ first slot and the orthogonal partner in the second. The closed-form
131+ angle ``θ = ½·arctan2(2ab, a²−b²)`` (``a``, ``b`` the two modes'
132+ components in the target FEM DOF) maximises that DOF's content in
133+ the first slot, identical in form to the FA/SS resolver in
134+ ``pybmodes.elastodyn.params``.
135+
136+ The rotation is accepted only when it *cleanly* separates the pair —
137+ both rotated modes dominated (``_DOMINANCE_THRESHOLD``) by *different*
138+ platform DOFs. Otherwise the pair is genuinely coupled (an asymmetric
139+ floater whose degeneracy is broken anyway, so this branch rarely
140+ triggers there) and the originals are kept for the global assignment
141+ to handle. The input is not mutated; a fresh copy is returned.
142+ """
143+ out = base .copy ()
144+ n = out .shape [1 ]
145+ i = 0
146+ while i < n - 1 :
147+ denom = max (abs (float (freqs [i ])), abs (float (freqs [i + 1 ])), 1e-12 )
148+ if abs (float (freqs [i ]) - float (freqs [i + 1 ])) / denom <= _DEGENERATE_FREQ_RTOL :
149+ bi = out [:, i ].copy ()
150+ bj = out [:, i + 1 ].copy ()
151+ # Target the pair's combined-dominant platform DOF; f1 is the
152+ # FEM DOF that platform DOF maps from.
153+ combined = _platform_dof_energy (bi , Mp ) + _platform_dof_energy (bj , Mp )
154+ f1 = int (_FEM_TO_PLATFORM [int (np .argmax (combined ))])
155+ a , b = float (bi [f1 ]), float (bj [f1 ])
156+ theta = 0.5 * float (np .arctan2 (2.0 * a * b , a * a - b * b ))
157+ c , s = np .cos (theta ), np .sin (theta )
158+ ri = c * bi + s * bj
159+ rj = - s * bi + c * bj
160+
161+ ei = _platform_dof_energy (ri , Mp )
162+ ej = _platform_dof_energy (rj , Mp )
163+ ti , tj = float (ei .sum ()), float (ej .sum ())
164+ if ti > 0.0 and tj > 0.0 :
165+ fi , fj = ei / ti , ej / tj
166+ ki , kj = int (np .argmax (fi )), int (np .argmax (fj ))
167+ if (
168+ ki != kj
169+ and fi [ki ] >= _DOMINANCE_THRESHOLD
170+ and fj [kj ] >= _DOMINANCE_THRESHOLD
171+ ):
172+ out [:, i ], out [:, i + 1 ] = ri , rj
173+ i += 2
174+ continue
175+ i += 1
176+ return out
177+
72178
73179def classify_platform_modes (
74180 eigvecs : np .ndarray ,
75181 active_dofs : np .ndarray ,
76182 nselt : int ,
77183 platform_mass : np .ndarray ,
184+ frequencies : np .ndarray | None = None ,
78185) -> list [str | None ]:
79186 """Return a per-mode label list naming the platform rigid-body
80187 modes (``surge`` / … / ``yaw``) or ``None`` where no single
@@ -94,11 +201,21 @@ def classify_platform_modes(
94201 DOF order (``PlatformND.mass`` from
95202 :func:`pybmodes.fem.nondim.nondim_platform`). Supplies the
96203 mass / inertia metric for the energy weighting.
204+ frequencies : (n_modes,) modal frequencies (any unit), ascending.
205+ When given, frequency-degenerate rigid-body pairs (surge≈sway,
206+ roll≈pitch on a symmetric platform) are rotated onto axis-
207+ aligned directions before labelling so the result is
208+ deterministic regardless of the arbitrary basis the eigensolver
209+ returns within a degenerate eigenspace. When ``None`` the
210+ rotation is skipped (the labels then rely on the global
211+ assignment alone).
97212
98213 Caller must invoke this only for a floating model
99214 (``hub_conn == 2`` with a ``PlatformSupport``); for any other
100215 model there are no rigid-body modes to name.
101216 """
217+ from scipy .optimize import linear_sum_assignment
218+
102219 ndt = NESH * nselt + 6
103220 n_modes = eigvecs .shape [1 ]
104221
@@ -114,40 +231,43 @@ def classify_platform_modes(
114231 # candidates are the first _N_RIGID columns.
115232 n_rigid = min (_N_RIGID , n_modes )
116233
117- labels : list [str | None ] = []
118- used : set [str ] = set ()
119- for m in range (n_modes ):
120- if m >= n_rigid :
121- labels .append (None ) # flexible tower mode
122- continue
123- b = base [:, m ]
124- # Per-DOF modal kinetic energy contribution b_i·(M_p b)_i. The
125- # platform mass matrix is the metric that puts translation and
126- # rotation amplitudes on a comparable footing.
127- e = np .abs (b * (Mp @ b ))
234+ labels : list [str | None ] = [None ] * n_modes
235+ if n_rigid == 0 :
236+ return labels
237+
238+ # Resolve symmetric-platform degeneracies first: rotate surge≈sway
239+ # and roll≈pitch pairs onto axis-aligned directions so the labelling
240+ # doesn't depend on which (equally valid) rotation the eigensolver
241+ # happened to return inside the degenerate eigenspace.
242+ work = base [:, :n_rigid ].astype (float , copy = True )
243+ if frequencies is not None :
244+ work = _align_degenerate_rigid_pairs (
245+ work , np .asarray (frequencies , dtype = float )[:n_rigid ], Mp
246+ )
247+
248+ # Per-DOF modal-kinetic-energy fractions for each rigid-body
249+ # candidate, in platform-DOF order: score[m, k] is the fraction of
250+ # mode m's base-node modal kinetic energy carried by platform DOF k.
251+ score = np .zeros ((n_rigid , 6 ))
252+ for m in range (n_rigid ):
253+ e = _platform_dof_energy (work [:, m ], Mp )
128254 total = float (e .sum ())
129255 if total <= 0.0 or not np .isfinite (total ):
130- labels .append (None )
131- continue
132- frac = e [_FEM_TO_PLATFORM ] / total # platform-DOF order
133- k = int (np .argmax (frac ))
134- if frac [k ] < _DOMINANCE_THRESHOLD :
135- labels .append (None ) # coupled / rotated pair
136- continue
137- name = _PLATFORM_DOF_NAMES [k ]
138- if name in used :
139- # A 6-DOF platform has exactly one rigid-body mode per DOF.
140- # Seeing the same dominant DOF twice within the lowest six
141- # means the one-mode-per-DOF assumption has already failed
142- # (a degenerate / rotated pair, or a flexible mode leaking
143- # in) — emitting a suffixed physical label like "surge (2)"
144- # would mislead a downstream plot/report into treating the
145- # duplicate as meaningful. Stay conservative: leave it
146- # ``None``, consistent with this classifier's "name only
147- # the unambiguous" contract.
148- labels .append (None )
149- continue
150- used .add (name )
151- labels .append (name )
256+ continue # inert row → stays zero
257+ score [m ] = e / total
258+
259+ # Global one-mode-per-DOF matching (Hungarian, maximising total
260+ # energy fraction), so the mode that best expresses each platform
261+ # DOF wins that label instead of an earlier mode greedily stealing
262+ # it and starving the true owner (issue #93). The square 6×6 case
263+ # gives a perfect matching, so no DOF is ever named twice; a short
264+ # rigid block (fewer than six modes requested) matches the
265+ # available rows. Each assignment is then gated by the dominance
266+ # threshold: a genuinely coupled / rotated pair whose energy splits
267+ # across DOFs falls below it and is left ``None``.
268+ rows , cols = linear_sum_assignment (score , maximize = True )
269+ for m , k in zip (rows , cols ):
270+ if score [m , k ] >= _DOMINANCE_THRESHOLD :
271+ labels [m ] = _PLATFORM_DOF_NAMES [k ]
152272
153273 return labels
0 commit comments