API reference 

Notes

Arbitrary discrete mid-span point masses are not yet modelled (a separate FEM-assembly extension with its own validation track — see issue #35); outfitting_factor covers the dominant distributed non-structural mass and tip_mass the tower-top lump.

classmethod from_windio(yaml_path, *, component='tower', thickness_interp='linear', hub_conn=1, tip_mass=None, n_nodes=None)[source]

Build a tower (or monopile) model directly from a WindIO ontology .yaml (issue #35).

Parses the structural subset — components.<component>’s outer_shape.outer_diameter, structure.layers wall thickness, structure.outfitting_factor, reference_axis — plus the referenced entry in the top-level materials list, and feeds it to from_geometry().

Parameters:

yaml_path (path to a WindIO ontology file.)
component ("tower" (default) or "monopile".)
thickness_interp (how a layer thickness grid maps onto the) – FEM stations — "linear" (WindIO-native piecewise- linear, default) or "piecewise_constant" (WISDEM-style constant-per-segment). The choice measurably moves the 2nd tower-bending polynomial coefficients; see tests/test_windio.py.
hub_conn (root BC (default 1 cantilever; use 3 for a) – soil-flexible monopile).
tip_mass (optional tower-top RNA lump (issue #35) — a) – pybmodes.io.bmi.TipMassProps or a bare float (RNA mass in kg). Replaces the tower._bmi.tip_mass = … workaround and mirrors from_windio_floating()’s rna_tip. None -> zero tip mass.
n_nodes (optional FE-mesh refinement (issue #35) — re-grid) – the tower onto n_nodes evenly-spaced stations (geometry linearly interpolated, properties recomputed exactly), to resolve higher tower-bending mode shapes. None keeps the WindIO grid. The WindIO blade path has the analogous n_span.

Return type:

Notes

Requires the optional [windio] extra (PyYAML). This is the tubular tower / monopile path; for a WindIO blade composite layup use pybmodes.models.RotatingBlade.from_windio() (PreComp-class thin-wall cross-section reduction).

classmethod from_windio_floating(yaml_path, *, component_tower='tower', water_depth=None, hydrodyn_dat=None, moordyn_dat=None, elastodyn_dat=None, platform_support=None, rna_tip=None, n_nodes=None, rho=1025.0, g=9.80665)[source]

Coupled floating tower+platform model from a WindIO .yaml (issue #35).

The WindIO-native analogue of from_elastodyn_with_mooring(). The tower beam always comes from components.tower (the validated Phase-1 tubular path — machine-exact vs the ElastoDyn tower). The platform has two fidelity tiers:

Industry-grade (companion decks present). When a hydrodyn_dat / moordyn_dat / elastodyn_dat is supplied, that leg uses the complete deck model — WAMIT A_inf + C_hst, the full MoorDyn system (its own anchor/fairlead geometry and line properties), and the ElastoDyn PtfmMass/RIner (incl. trim ballast) + lumped RNA + draft convention. With all three present this is byte-identical to the BModes-JJ-validated from_elastodyn_with_mooring() except the tower is the WindIO one (≈ 0.0003 % reference grade).
Screening preview (yaml-only legs). Any leg without a deck falls to the WindIO model — member-waterplane C_hst (geometry-exact, ≈ 1.6 %), Morison + end-cap A_inf (heave is screening-only — Morison ≠ potential flow, as RAFT/WISDEM also find), WindIO catenary mooring geometry, and structural + fixed-ballast inertia (no trim ballast). A UserWarning names every screening leg; it is not industry-grade for the platform and is for fast pre-deck previewing only.
Injected platform (``platform_support=``). When a pybmodes.io.bmi.PlatformSupport is passed, the floater is taken verbatim from it — its own A_inf / C_hst / mooring_K / 6×6 inertia / draft / ref_msl. The tower geometry still comes from the WindIO .yaml; nothing about the floating substructure is read from the yaml or any deck. This is the “floater designed separately” workflow — a frequency-domain tool / WAMIT export / published 6×6 set feeding the same BModes-JJ-validated free-free FEM that reproduces OC3 Hywind to ≈ 0.0003 %. The tower beam length stays the WindIO flexible_length independent of the supplied draft (the radius + draft cancellation the deck path also relies on — see make_params). No screening warning (the caller owns the platform fidelity). Mutually exclusive with the companion decks; optionally pass rna_tip for the tower-top RNA lump (default: bare tower top, no RNA). The frequencies inherit the validated PlatformSupport assembly; this path adds no new numerics.

n_nodes re-grids the tower beam onto that many evenly-spaced stations (geometry linearly interpolated, closed-form tube properties recomputed exactly), mirroring from_windio() / from_geometry() (issue #58 — uniform mesh-refinement kwarg across the WindIO/geometry constructors). None keeps the WindIO grid. It refines only the tower discretisation; the platform assembly is unaffected.

Sets hub_conn = 2 / tow_support = 1 and reuses the existing BModes-JJ-validated free-free PlatformSupport FEM unchanged. Needs the optional [windio] extra. For ElastoDyn polynomial generation use the cantilever from_windio() regardless of platform (see cases/ECOSYSTEM_FINDING.md).

Parameters:

yaml_path (str | pathlib.Path)
component_tower (str)
water_depth (float | None)
hydrodyn_dat (str | pathlib.Path | None)
moordyn_dat (str | pathlib.Path | None)
elastodyn_dat (str | pathlib.Path | None)
platform_support (PlatformSupport | None)
rna_tip (TipMassProps | None)
n_nodes (int | None)
rho (float)
g (float)

Return type:

classmethod from_windio_with_monopile(yaml_path, *, component_tower='tower', component_monopile='monopile', thickness_interp='linear', tip_mass=None, n_nodes=None)[source]

Build a combined monopile + tower fixed-bottom cantilever from a WindIO ontology .yaml (issue #92).

from_windio() reduces a single tube; this constructor reduces the monopile and tower components separately (each keeps its own wall schedule and steel grade) and splices them bottom-to-top at the transition piece — the elevation where the monopile top meets the tower base — into one cantilever clamped at the mudline (hub_conn = 1), with the RNA lumped at the tower top via tip_mass. It is the WindIO analog of from_elastodyn_with_subdyn() (the ElastoDyn + SubDyn splice).

Parameters:

yaml_path (path to a WindIO ontology file carrying both a) – monopile and a tower component.
component_tower (component names to splice) – (defaults "tower" / "monopile").
component_monopile (component names to splice) – (defaults "tower" / "monopile").
thickness_interp ("linear" (default) or) – "piecewise_constant" — see from_windio().
tip_mass (optional tower-top RNA lump — a) – pybmodes.io.bmi.TipMassProps or a bare float (RNA mass in kg). None -> zero tip mass.
n_nodes (optional FE-mesh refinement applied per segment) – (each of the monopile and tower is re-gridded onto n_nodes evenly-spaced stations), mirroring from_windio()’s n_nodes. None keeps each component’s native WindIO grid.

Return type:

Notes

Requires the optional [windio] extra (PyYAML). This is the rigid fixed-base monopile path: the pile is clamped at the mudline with no soil flexibility, matching from_elastodyn_with_subdyn() and the bundled monopile samples. Distributed soil springs (a Winkler distr_k / hub_conn = 3 foundation) and Morison hydrodynamics are out of scope here and tracked separately. Raises ValueError if the monopile top and tower base do not meet at a common transition-piece elevation.

run(n_modes=20, *, check_model=True, on_error='raise')[source]

Solve the eigenvalue problem and return frequencies + mode shapes.

Parameters:

n_modes (number of modes to extract (must be >= 1; default 20).)
check_model (run pybmodes.checks.check_model() before the) – solve (default True). INFO findings are silent (call pybmodes.checks.check_model(model) explicitly to see those). Pass check_model=False to skip the pre-solve checks for scripted callers that have already validated their inputs.
on_error (how ERROR-severity findings are handled when) – check_model runs (default "raise", 1.14.0). ERROR findings flag non-physical input (NaN section properties, non-positive mass, a malformed support matrix), so the solve fails closed by raising pybmodes.checks.ModelValidationError rather than feeding the eigensolver garbage. Pass on_error="warn" to downgrade ERROR findings to UserWarning and continue, the pre-1.14.0 behaviour. WARN findings always emit as UserWarning regardless.

Return type:

ModalResult

Warning

n_modes affects the LAPACK solver path. For symmetric or nearly-symmetric towers (EI_FA ≈ EI_SS and small RNA c.m. offset), use n_modes >= 6. With n_modes <= 4, scipy.linalg.eigh invokes a subset eigenvalue routine that can artificially lift the degeneracy of the 1st FA / SS bending pair — the modes come back at slightly different frequencies and pre-separated, which prevents the degenerate-pair resolver in pybmodes.elastodyn.params from triggering. The polynomial fits still succeed, but the FA / SS classification may flip relative to a full solve, and downstream compute_tower_params_report may select different modes for TwFAM1Sh / TwSSM1Sh between runs at different n_modes.

Minimum recommended: n_modes >= 6 for reliable FA / SS classification on symmetric structures. The default of 20 is safely above this threshold.

ModalResult dataclass returned by RotatingBlade and Tower.

class ModalResult(frequencies, shapes, participation=None, fit_residuals=None, metadata=None, mode_labels=None, ignored_physics=(), diagnostics=None)[source]

Bases: object

Frequencies and mode shapes from a single FEM solve.

Variables:

frequencies ((n_modes,) array of natural frequencies in Hz.)
shapes (list of NodeModeShape, one per mode, ordered root) – to tip.
participation ((n_modes, 3) array of energy fractions in the) – per-mode (flap or FA, edge or SS, torsion) axes — populated by downstream classification code such as pybmodes.campbell._participation(). Each row sums to 1. None when not yet computed; included in the saved archive only when set.
fit_residuals (optional {block_name: rms_value} dict of) – polynomial-fit RMS residuals — populated by pybmodes.elastodyn.compute_tower_params() / compute_blade_params callers that want to embed the fit quality in the serialised result. None when not set.
metadata (optional metadata dict (pyBmodes version, source file,) – save timestamp, git hash) attached automatically by save() / to_json() if not already populated.
mode_labels (optional per-mode classification labels, one entry) – per mode (parallel to shapes / frequencies). For a floating model (hub_conn = 2 with a PlatformSupport) the platform rigid-body modes are named "surge" / "sway" / "heave" / "roll" / "pitch" / "yaw"; a mode the classifier can’t confidently attribute to a single platform DOF (a flexible tower mode, or a strongly coupled / rotated pair) is left as None. The whole list is None for a non-floating model (cantilever / monopile have no rigid-body modes to name). Added 1.3.0 — the last dataclass field, so the pre-1.3.0 positional constructor ABI is preserved (see the field-order note below); included in the saved archive only when set, like participation / fit_residuals.
ignored_physics (tuple of short labels naming any parsed-but-not-) – modelled physics the solve silently dropped, e.g. "distributed added mass (distr_m)". Empty when the model was solved at full fidelity. Lets a downstream consumer of a saved result see that it is an approximation rather than guess. Added 1.14.0; persisted when non-empty.
diagnostics (optional SolverDiagnostics) – for the solve that produced this result (path taken, sparse-to- dense fallback, mode-count guarantee, per-mode residuals, mass-matrix conditioning). Attached by run_fem(); not serialised (it describes the live solve run, not the persistent result data) and excluded from equality. Added 1.14.0.

Parameters:

frequencies (np.ndarray)
shapes (list[NodeModeShape])
participation (np.ndarray | None)
fit_residuals (dict[str, float] | None)
metadata (dict[str, Any] | None)
mode_labels (list[str | None] | None)
ignored_physics (tuple[str, ...])
diagnostics (SolverDiagnostics | None)

diagnostics: SolverDiagnostics | None = None

fit_residuals: dict[str, float] | None = None

frequencies: np.ndarray

classmethod from_json(path)[source]

Read a result back from a JSON file saved by to_json().

Parameters:: path (str | Path)
Return type:: ModalResult

ignored_physics: tuple[str, ...] = ()

classmethod load(path, *, allow_legacy_pickle=False)[source]

Read a result back from a .npz archive saved by save(). The reconstructed instance is value-equal to the original modulo numpy dtype promotion.

Archives written by current pyBmodes are pickle-free and load under allow_pickle=False. A legacy pre-1.0 archive whose __meta__ is a pickled object array is refused by default (object-array unpickling can execute arbitrary code); pass allow_legacy_pickle=True to opt in for a file you trust.

Parameters:

path (str | Path)
allow_legacy_pickle (bool)

Return type:

ModalResult

metadata: dict[str, Any] | None = None

mode_labels: list[str | None] | None = None

participation: np.ndarray | None = None

save(path, *, source_file=None)[source]

Write the result to a .npz archive.

source_file is recorded in the metadata when supplied (typically the BMI / ElastoDyn deck the solve came from).

Parameters:

path (str | Path)
source_file (str | Path | None)

Return type:

None

shapes: list[NodeModeShape]

to_json(path, *, source_file=None)[source]

Write the result to a JSON file. Arrays are emitted as nested lists; metadata is embedded under "metadata".

Parameters:

path (str | Path)
source_file (str | Path | None)

Return type:

None

Input / output

Parse errors (common across every `pybmodes.io.*` reader)

Unified parse-error base class for every pybmodes.io.* reader.

Every input-format parser pybmodes ships (BModes .bmi / section-properties .dat / OpenFAST ElastoDyn / SubDyn / HydroDyn / MoorDyn .dat / WAMIT .1 .hst / BModes .out reference output / WISDEM / WindIO ontology .yaml) previously raised either a bare ValueError with a format-specific message or, in the case of pybmodes.io.out_parser, the legacy BModeOutParseError. That made downstream try / except ValueError callers correct but indiscriminate, and the per-error file / line context was unstructured prose.

ParseError is the unified base. It inherits ValueError so existing except ValueError callers continue to catch every parse error untouched — the inheritance addition is non-breaking.

Each parser will (incrementally) start raising ParseError subclasses with structured file / line / column / context fields. Callers that want file / line context can switch to except ParseError and read the typed fields; callers that only need “the file is broken” stay on except ValueError.

Examples

Catch a parse error and pull out the file / line / column context:

from pybmodes.io.errors import ParseError
from pybmodes.io.bmi import read_bmi

try:
    bmi = read_bmi("malformed.bmi")
except ParseError as err:
    print(f"parse failed in {err.file}:{err.line}: {err}")
    print(err.format_diagnostic())   # one-line diagnostic string

Keep an older except ValueError pattern working — ParseError is a subclass so the catch still triggers:

try:
    bmi = read_bmi("malformed.bmi")
except ValueError as err:
    print(f"oops: {err}")   # also catches every ParseError

exception BMIParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ParseError

Raised by pybmodes.io.bmi.read_bmi() and the companion section-properties parser when the input deck is malformed.

The parser is line-oriented; line is the 1-based row in the source file. context is the offending line (truncated).

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

exception ElastoDynParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ParseError

Raised by the ElastoDyn deck reader (pybmodes.io.elastodyn_reader / the private pybmodes.io._elastodyn sub-package) on malformed input.

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

exception MoorDynParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ParseError

Raised by pybmodes.mooring.MooringSystem.from_moordyn() when a MoorDyn .dat carries an unrecognised layout or a point-ID column ordering the parser can’t auto-detect.

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

exception ParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ValueError

Base class for every pybmodes.io.* parser exception.

Inherits ValueError so except ValueError catches it unchanged. Inheriting parsers add structured context fields below.

Variables:

message (str) – The human-readable description of the failure. Same as str(err).
file (pathlib.Path-compatible str | None) – The source file the parser was reading from, when known. None for in-memory parses (e.g. yaml from a string).
line (int | None) – 1-based line number where the error was detected, when known. The parser is encouraged to populate this even for token-level errors — fall back to the start of the containing block.
column (int | None) – 1-based column number where the error was detected, when known. Optional; many parsers don’t track column.
context (str | None) – A short snippet — typically the offending line or token — that gives the reader visual confirmation of what the parser tripped on. Up to one or two lines; the parser should truncate large blobs.

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

column: int | None = None

context: str | None = None

file: str | None = None

format_diagnostic()[source]

Render a one-line file:line:column — message diagnostic.

Useful as a fallback when callers want a uniform error format regardless of which parser raised. Falls back to just the message when no file / line context is available.

Return type:: str

line: int | None = None

message: str

exception SubDynParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ParseError

Raised by pybmodes.io.subdyn_reader when the SubDyn joints / members / reaction-joint block can’t be parsed.

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

exception WAMITParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ParseError

Raised by pybmodes.io.wamit_reader.HydroDynReader and the underlying .1 / .hst readers on malformed WAMIT output (bad re-dimensionalisation, missing files behind PotFile, asymmetric matrices that can’t be mirrored).

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

exception WindIOParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ParseError

Raised by the WindIO ontology readers (pybmodes.io.windio / pybmodes.io.windio_blade / pybmodes.io.windio_floating) on schema-drift or malformed published blocks.

Includes both components.<component> lookup failures and the “elastic_properties block is present but unparseable” path that elastic="auto" / "file" distinguish.

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

BModes `.bmi`

Parser for .bmi main input files.

The .bmi format is line-ordered: values precede their labels. The reader follows the format conventions:

ReadCom  -> consume one line verbatim (section headers / blank lines)
ReadStr  -> consume one line, return it as a string
ReadVar  -> skip blanks, return first whitespace token of next non-blank line
ReadAry  -> skip blanks, return first N tokens of next non-blank line

class BMIFile(title, echo, beam_type, rot_rpm, rpm_mult, radius, hub_rad, precone, bl_thp, hub_conn, n_modes_print, tab_delim, mid_node_tw, tip_mass, id_mat, sec_props_file, scaling, n_elements, el_loc, tow_support=0, support=None, source_file=None)[source]

Bases: object

All parameters parsed from a .bmi main input file.

See Coordinate systems, origins & datums for the coordinate system, origin (tower base) and boundary-condition definitions.

Key fields

beam_typeint: 1 = blade (rotating), 2 = tower.
radiusfloat: Span length of the finite-element beam, metres — the blade radius or the flexible tower length (TowerHt - TowerBsHt). The beam runs from base (z = 0) to this length; the absolute base elevation is conveyed separately (via the PlatformSupport.draft for a floater).
hub_radfloat: Hub radius (blade root radial offset from the rotation axis), m.
hub_connint: Base boundary condition: 1 cantilever (clamped), 2 free-free floating (reactions from support), 3 soft monopile (lateral + rocking free), 4 pinned-free (bending slopes free). See the hub_conn table in Coordinate systems, origins & datums.
rot_rpm, preconefloat: Rotor speed (rpm) and coning angle (deg) for a rotating blade.
tip_massTipMassProps: Tower-top RNA / blade-tip lumped mass.
tow_supportint: 0 none / wires-only, 1 inline platform block, normalised to a PlatformSupport by the parser.
supportTensionWireSupport | PlatformSupport | None: The offshore platform or guy-wire support, when present.

beam_type: int

bl_thp: float

echo: bool

el_loc: ndarray

hub_conn: int

hub_rad: float

id_mat: int

mid_node_tw: bool

n_elements: int

n_modes_print: int

precone: float

radius: float

resolve_sec_props_path()[source]

Return absolute path to the section properties file.

Return type:: Path

rot_rpm: float

rpm_mult: float

scaling: ScalingFactors

sec_props_file: str

source_file: Path | None = None

support: TensionWireSupport | PlatformSupport | None = None

tab_delim: bool

tip_mass: TipMassProps

title: str

tow_support: int = 0

Parameters:

title (str)
echo (bool)
beam_type (int)
rot_rpm (float)
rpm_mult (float)
radius (float)
hub_rad (float)
precone (float)
bl_thp (float)
hub_conn (int)
n_modes_print (int)
tab_delim (bool)
mid_node_tw (bool)
tip_mass (TipMassProps)
id_mat (int)
sec_props_file (str)
scaling (ScalingFactors)
n_elements (int)
el_loc (ndarray)
tow_support (int)
support (TensionWireSupport | PlatformSupport | None)
source_file (Path | None)

class PlatformSupport(draft, cm_pform, mass_pform, i_matrix, ref_msl, hydro_M, hydro_K, mooring_K, distr_m_z, distr_m, distr_k_z, distr_k, wires=None, cm_pform_x=0.0, cm_pform_y=0.0, ref_x=0.0, ref_y=0.0)[source]

Bases: object

Offshore floating-platform or monopile support.

All quantities follow the conventions in Coordinate systems, origins & datums — read that page if anything below is unclear. The single origin is the tower base; the vertical datum is mean sea level (MSL, z = 0); the 6×6 matrices use OpenFAST DOF order (surge, sway, heave, roll, pitch, yaw — see pybmodes.coords).

Variables:

draft (float) –
Signed elevation of the tower base relative to MSL, in metres, negative = above MSL (a tower base 15 m above the waterline is draft = -15).

Note

Despite the name, this is not the naval-architecture draft (the keel depth below the waterline). The field name is inherited verbatim from the BModes .bmi platform block; here it means “tower-base z relative to MSL, negative up”. pyBmodes forms the CM→tower-base vertical lever internally as cm_pform - draft, so the base elevation enters the maths exactly once — do not also add it to cm_pform / ref_msl.
cm_pform (float) – Platform centre-of-mass depth below MSL, metres, positive downward. The vertical reference of i_matrix.
mass_pform (float) – Platform (substructure) mass, kg.
i_matrix (np.ndarray) – Platform 6×6 rigid-body inertia (kg, kg·m²) about the platform CM, in OpenFAST DOF order. Transferred to the tower base by the rigid-arm transform using cm_pform (vertical) and cm_pform_x / cm_pform_y (horizontal).
ref_msl (float) – Depth below MSL of the reference point for the hydro and mooring matrices, metres, positive downward (the WAMIT/HydroDyn PtfmRefzt; usually 0, i.e. at the waterline). This reference is assumed to lie on the tower axis (PtfmRefxt = PtfmRefyt = 0); hydro_* / mooring_K receive no horizontal arm (see cm_pform_x).
hydro_M (np.ndarray) – Infinite-frequency hydrodynamic added mass A_inf, 6×6, about the ref_msl reference, OpenFAST DOF order.
hydro_K (np.ndarray) – Hydrostatic restoring C_hst, 6×6, about the ref_msl reference. May be (legitimately) negative-definite in roll/pitch for an unstable spar stabilised by mooring.
mooring_K (np.ndarray) – Linearised mooring stiffness, 6×6, about the ref_msl reference, OpenFAST DOF order.
distr_m (distr_m_z,) – Optional distributed added mass vs depth (Morison), if used.
distr_k (distr_k_z,) – Optional distributed Winkler soil/foundation stiffness vs depth (the hub_conn = 3 soft-monopile path).
wires (TensionWireSupport | None) – Optional guy-wire support.
cm_pform_y (cm_pform_x,) – Horizontal offset of the platform CM from the tower axis (applied to the inertia transform; see the inline note below).
ref_y (ref_x,) – Horizontal position of the hydro/mooring reference point (PtfmRefxt / PtfmRefyt) from the tower axis, metres. Default 0.0 (reference on the tower axis — every standard deck). Set non-zero for an off-axis floater so hydro_* / mooring_K are carried horizontally to the tower base too (issue #100, 1.13.0).
tower_base_z (float (property)) – Intuitive positive-up alias for draft (tower_base_z == -draft).

Parameters:

draft (float)
cm_pform (float)
mass_pform (float)
i_matrix (ndarray)
ref_msl (float)
hydro_M (ndarray)
hydro_K (ndarray)
mooring_K (ndarray)
distr_m_z (ndarray)
distr_m (ndarray)
distr_k_z (ndarray)
distr_k (ndarray)
wires (TensionWireSupport | None)
cm_pform_x (float)
cm_pform_y (float)
ref_x (float)
ref_y (float)

cm_pform: float

cm_pform_x: float = 0.0

cm_pform_y: float = 0.0

distr_k: ndarray

distr_k_z: ndarray

distr_m: ndarray

distr_m_z: ndarray

draft: float

hydro_K: ndarray

hydro_M: ndarray

i_matrix: ndarray

mass_pform: float

mooring_K: ndarray

ref_msl: float

ref_x: float = 0.0

ref_y: float = 0.0

property tower_base_z: float

Elevation of the tower base above MSL, metres, positive up.

Intuitive alias for draft, which is the BModes-inherited signed, negative-up spelling: tower_base_z == -draft. A tower base 15 m above the waterline is tower_base_z = 15 (equivalently draft = -15). Reading or assigning one keeps the other in sync. Added 1.13.0 (issue #100).

wires: TensionWireSupport | None = None

class ScalingFactors(sec_mass=1.0, flp_iner=1.0, lag_iner=1.0, flp_stff=1.0, edge_stff=1.0, tor_stff=1.0, axial_stff=1.0, cg_offst=1.0, sc_offst=1.0, tc_offst=1.0)[source]

Bases: object

Multiplicative scaling factors applied to all section properties.

Parameters:

sec_mass (float)
flp_iner (float)
lag_iner (float)
flp_stff (float)
edge_stff (float)
tor_stff (float)
axial_stff (float)
cg_offst (float)
sc_offst (float)
tc_offst (float)

axial_stff: float = 1.0

cg_offst: float = 1.0

edge_stff: float = 1.0

flp_iner: float = 1.0

flp_stff: float = 1.0

lag_iner: float = 1.0

sc_offst: float = 1.0

sec_mass: float = 1.0

tc_offst: float = 1.0

tor_stff: float = 1.0

class TensionWireSupport(n_attachments, n_wires, node_attach, wire_stiffness, th_wire)[source]

Bases: object

Tension wire (guy wire) support for land-based towers.

Parameters:

n_attachments (int)
n_wires (list[int])
node_attach (list[int])
wire_stiffness (list[float])
th_wire (list[float])

n_attachments: int

n_wires: list[int]

node_attach: list[int]

th_wire: list[float]

wire_stiffness: list[float]

class TipMassProps(mass, cm_offset, cm_axial, ixx, iyy, izz, ixy, izx, iyz)[source]

Bases: object

Mass and inertia of the blade tip or tower-top concentrated mass.

Lumps the rotor-nacelle assembly (RNA) at the tower top, or a tip mass at the blade tip. Offsets are in the span-end section frame: z along the span (positive toward / beyond the top), the transverse axis aligned with the lateral/sway direction. See Coordinate systems, origins & datums.

Variables:

mass (float) – Concentrated mass, kg. Omitting the RNA (mass = 0) on a tower makes the 1st fore-aft frequency ~10–30 % too high.
cm_offset (float) – CM offset transverse to the span (along the tip-section axis aligned with lateral/sway), metres.
cm_axial (float) – CM offset along the span (z), metres — how far the lumped CM sits beyond the span end (e.g. the RNA CM above the tower top).
iyz (ixx, iyy, izz, ixy, izx,) – Mass moments / products of inertia about the CM, kg·m².

Parameters:

mass (float)
cm_offset (float)
cm_axial (float)
ixx (float)
iyy (float)
izz (float)
ixy (float)
izx (float)
iyz (float)

cm_axial: float

cm_offset: float

ixx: float

ixy: float

iyy: float

iyz: float

izx: float

izz: float

mass: float

read_bmi(path)[source]

Parse a .bmi main input file and return a BMIFile.

Parameters:: path (str | Path)
Return type:: BMIFile

ElastoDyn deck reader

Parser, writer, and pyBmodes adapter for OpenFAST ElastoDyn .dat files.

ElastoDyn .dat files come in three flavours, all of which share a common line-ordered convention: each significant line carries a value (or short list) followed by a label keyword and an optional - comment tail. The parsers below are label-based rather than position-based, which keeps them robust across the FAST v8 → OpenFAST v3+ format drift (renamed/added/removed scalars between versions; BldFile(1) vs BldFile1).

Three entry points cover the three file flavours:

read_elastodyn_main() — the top-level ElastoDyn input file.
read_elastodyn_tower() — the tower-properties file referenced via TwrFile.
read_elastodyn_blade() — the blade-properties file referenced via BldFile(1..3).

Each returns a dataclass holding the parsed values plus enough raw- line metadata to re-emit a semantically identical file via the matching write_* function.

Two adapter helpers turn a parsed ElastoDyn bundle into pyBmodes BMIFile + SectionProperties records ready for the FEM core:

to_pybmodes_tower(main, tower, blade=None)()
to_pybmodes_blade(main, blade)()

The implementation lives in the private sub-package pybmodes.io._elastodyn:

types — dataclasses
lex — line/token scanning helpers
parser — line-driven flavour parsers
writer — canonical re-emitters
adapter — pyBmodes BMI/SectionProperties synthesisers

This module re-exports the public names from those sub-modules so existing from pybmodes.io.elastodyn_reader import … imports keep working unchanged. Field-set discrepancies vs. the user-specified spec are documented under each dataclass.

Round-trip contract

write_* functions emit a canonically formatted file that parses back to an equal dataclass but is not byte-identical to the original. Whitespace, label column position, and comment text are normalised. The test suite compares the parse-emit-reparse fixed point (see tests/test_elastodyn_reader.py).

class ElastoDynBlade(header, title, source_file=None, n_bl_inp_st=0, bld_fl_dmp=<factory>, bld_ed_dmp=<factory>, fl_st_tunr=<factory>, adj_bl_ms=1.0, adj_fl_st=1.0, adj_ed_st=1.0, bl_fract=<factory>, pitch_axis=None, strc_twst=<factory>, b_mass_den=<factory>, flp_stff=<factory>, edg_stff=<factory>, rotary_inertia_available=False, bld_fl1_sh=<factory>, bld_fl2_sh=<factory>, bld_edg_sh=<factory>, distr_header_lines=<factory>, section_dividers=<factory>)[source]

Bases: object

Parsed ElastoDyn blade input file.

ElastoDyn ships only translational mass density (BMassDen) and bending stiffnesses (FlpStff, EdgStff) per spanwise station; it has no per-section rotary mass moments of inertia. Those live in BeamDyn or come from a cross-section pre-processor (VABS, PreComp).

The rotary_inertia_available flag is therefore always False after parsing an ElastoDyn blade file. Downstream code (to_pybmodes_blade) treats the rotary inertia contributions as zero, which is the correct Euler-Bernoulli limit for slender blades and is sufficient for the bending modes (1–4 flap/edge) pyBmodes targets. A tiny regularisation floor is added in the section-property builder to keep the global mass matrix positive-definite without fabricating physically meaningful rotary terms.

Parameters:

header (str)
title (str)
source_file (Path | None)
n_bl_inp_st (int)
bld_fl_dmp (list[float])
bld_ed_dmp (list[float])
fl_st_tunr (list[float])
adj_bl_ms (float)
adj_fl_st (float)
adj_ed_st (float)
bl_fract (ndarray)
pitch_axis (ndarray | None)
strc_twst (ndarray)
b_mass_den (ndarray)
flp_stff (ndarray)
edg_stff (ndarray)
rotary_inertia_available (bool)
bld_fl1_sh (ndarray)
bld_fl2_sh (ndarray)
bld_edg_sh (ndarray)
distr_header_lines (list[str])
section_dividers (list[str])

adj_bl_ms: float = 1.0

adj_ed_st: float = 1.0

adj_fl_st: float = 1.0

b_mass_den: ndarray

bl_fract: ndarray

bld_ed_dmp: list[float]

bld_edg_sh: ndarray

bld_fl1_sh: ndarray

bld_fl2_sh: ndarray

bld_fl_dmp: list[float]

distr_header_lines: list[str]

edg_stff: ndarray

fl_st_tunr: list[float]

flp_stff: ndarray

header: str

n_bl_inp_st: int = 0

pitch_axis: ndarray | None = None

rotary_inertia_available: bool = False

section_dividers: list[str]

source_file: Path | None = None

strc_twst: ndarray

title: str

class ElastoDynMain(header, title, source_file=None, num_bl=3, tip_rad=0.0, hub_rad=0.0, pre_cone=<factory>, hub_cm=0.0, overhang=0.0, shft_tilt=0.0, twr2shft=0.0, tower_ht=0.0, tower_bs_ht=0.0, nac_cm_xn=0.0, nac_cm_yn=0.0, nac_cm_zn=0.0, rot_speed_rpm=0.0, tip_mass=<factory>, hub_mass=0.0, hub_iner=0.0, gen_iner=0.0, nac_mass=0.0, nac_y_iner=0.0, yaw_br_mass=0.0, bld_file=<factory>, twr_file='', scalars=<factory>, out_list=<factory>, nodal_out_list=<factory>, section_dividers=<factory>)[source]

Bases: object

Parsed top-level ElastoDyn input file.

Parameters:

header (str)
title (str)
source_file (Path | None)
num_bl (int)
tip_rad (float)
hub_rad (float)
pre_cone (list[float])
hub_cm (float)
overhang (float)
shft_tilt (float)
twr2shft (float)
tower_ht (float)
tower_bs_ht (float)
nac_cm_xn (float)
nac_cm_yn (float)
nac_cm_zn (float)
rot_speed_rpm (float)
tip_mass (list[float])
hub_mass (float)
hub_iner (float)
gen_iner (float)
nac_mass (float)
nac_y_iner (float)
yaw_br_mass (float)
bld_file (list[str])
twr_file (str)
scalars (dict[str, str])
out_list (list[str])
nodal_out_list (list[str])
section_dividers (list[str])

bld_file: list[str]

compute_rot_mass(blade)[source]

Total rotor mass = hub + N · AdjBlMs · ∫ BMassDen ds along the blade.

Requires the blade file to integrate the distributed mass density. The AdjBlMs scalar from the blade file is an ElastoDyn-side multiplier on the entire mass distribution (the “blade-mass adjustment factor”); ignoring it under- / over-reports rotor mass on any deck where it deviates from 1.0. The adapter at to_pybmodes_blade() already applies it; this method previously did not.

Parameters:: blade (ElastoDynBlade)
Return type:: float

gen_iner: float = 0.0

header: str

hub_cm: float = 0.0

property hub_ht: float: Hub height above tower base, derived from geometry.

hub_iner: float = 0.0

hub_mass: float = 0.0

hub_rad: float = 0.0

nac_cm_xn: float = 0.0

nac_cm_yn: float = 0.0

nac_cm_zn: float = 0.0

nac_mass: float = 0.0

nac_y_iner: float = 0.0

nodal_out_list: list[str]

num_bl: int = 3

out_list: list[str]

overhang: float = 0.0

pre_cone: list[float]

rot_speed_rpm: float = 0.0

scalars: dict[str, str]

section_dividers: list[str]

shft_tilt: float = 0.0

source_file: Path | None = None

tip_mass: list[float]

tip_rad: float = 0.0

title: str

tower_bs_ht: float = 0.0

tower_ht: float = 0.0

twr2shft: float = 0.0

twr_file: str = ''

yaw_br_mass: float = 0.0

class ElastoDynTower(header, title, source_file=None, n_tw_inp_st=0, twr_fa_dmp=<factory>, twr_ss_dmp=<factory>, fa_st_tunr=<factory>, ss_st_tunr=<factory>, adj_tw_ma=1.0, adj_fa_st=1.0, adj_ss_st=1.0, ht_fract=<factory>, t_mass_den=<factory>, tw_fa_stif=<factory>, tw_ss_stif=<factory>, tw_fa_iner=<factory>, tw_ss_iner=<factory>, tw_fa_cg_of=<factory>, tw_ss_cg_of=<factory>, tw_fa_m1_sh=<factory>, tw_fa_m2_sh=<factory>, tw_ss_m1_sh=<factory>, tw_ss_m2_sh=<factory>, distr_header_lines=<factory>, section_dividers=<factory>)[source]

Bases: object

Parsed ElastoDyn tower input file.

Parameters:

header (str)
title (str)
source_file (Path | None)
n_tw_inp_st (int)
twr_fa_dmp (list[float])
twr_ss_dmp (list[float])
fa_st_tunr (list[float])
ss_st_tunr (list[float])
adj_tw_ma (float)
adj_fa_st (float)
adj_ss_st (float)
ht_fract (ndarray)
t_mass_den (ndarray)
tw_fa_stif (ndarray)
tw_ss_stif (ndarray)
tw_fa_iner (ndarray)
tw_ss_iner (ndarray)
tw_fa_cg_of (ndarray)
tw_ss_cg_of (ndarray)
tw_fa_m1_sh (ndarray)
tw_fa_m2_sh (ndarray)
tw_ss_m1_sh (ndarray)
tw_ss_m2_sh (ndarray)
distr_header_lines (list[str])
section_dividers (list[str])

adj_fa_st: float = 1.0

adj_ss_st: float = 1.0

adj_tw_ma: float = 1.0

distr_header_lines: list[str]

fa_st_tunr: list[float]

header: str

ht_fract: ndarray

n_tw_inp_st: int = 0

section_dividers: list[str]

source_file: Path | None = None

ss_st_tunr: list[float]

t_mass_den: ndarray

title: str

tw_fa_cg_of: ndarray

tw_fa_iner: ndarray

tw_fa_m1_sh: ndarray

tw_fa_m2_sh: ndarray

tw_fa_stif: ndarray

tw_ss_cg_of: ndarray

tw_ss_iner: ndarray

tw_ss_m1_sh: ndarray

tw_ss_m2_sh: ndarray

tw_ss_stif: ndarray

twr_fa_dmp: list[float]

twr_ss_dmp: list[float]

read_elastodyn_blade(path)[source]

Parameters:: path (str | Path)
Return type:: ElastoDynBlade

read_elastodyn_main(path)[source]

Parse a top-level ElastoDyn .dat file.

Parameters:: path (str | Path)
Return type:: ElastoDynMain

read_elastodyn_tower(path)[source]

Parameters:: path (str | Path)
Return type:: ElastoDynTower

to_pybmodes_blade(main, blade)[source]

Build pyBmodes BMIFile and SectionProperties for blade modal analysis at the operating RotSpeed from the main file.

Parameters:

main (ElastoDynMain)
blade (ElastoDynBlade)

Return type:

tuple[BMIFile, SectionProperties]

to_pybmodes_tower(main, tower, blade=None, *, physical_sec_props=False)[source]

Build pyBmodes BMIFile and SectionProperties for tower modal analysis from a parsed ElastoDyn bundle. blade is optional; when omitted, the rotor mass is approximated as HubMass only.

physical_sec_props selects the section-property synthesis (see _stack_tower_section_props()): False (default) for the clamped-base cantilever / monopile path; True for the free-base floating path, where the cantilever proxies wreck conditioning.

Parameters:

main (ElastoDynMain)
tower (ElastoDynTower)
blade (ElastoDynBlade | None)
physical_sec_props (bool)

Return type:

tuple[BMIFile, SectionProperties]

write_elastodyn_blade(obj, path=None)[source]

Parameters:

obj (ElastoDynBlade)
path (str | Path | None)

Return type:

str

write_elastodyn_main(obj, path=None)[source]

Re-emit a main ElastoDyn file. Returns the text and optionally writes it to path.

Parameters:

obj (ElastoDynMain)
path (str | Path | None)

Return type:

str

write_elastodyn_tower(obj, path=None)[source]

Parameters:

obj (ElastoDynTower)
path (str | Path | None)

Return type:

str

SubDyn reader

Minimal parser and adapter for OpenFAST SubDyn .dat substructure files.

Scope is intentionally narrow — only what’s needed to drive a clamped-base monopile model from the existing reference decks shipped with the OpenFAST regression-test corpus. Specifically:

JOINTS — XYZ positions and joint type.
BASE REACTION JOINTS — read the single clamped-base joint id and its DOF flags; SSI files are not handled (the 5MW OC3 deck has none).
INTERFACE JOINTS — read the transition-piece joint id.
MEMBERS — beam connectivity (joint pairs and property-set ids).
CIRCULAR BEAM CROSS-SECTION PROPERTIES — E, G, MatDens, outer diameter, wall thickness.

Sections we don’t parse: rectangular / arbitrary cross-sections, cables, rigid links, spring elements, cosine matrices, concentrated masses, output settings. These are zero in the bundled OC3 monopile and SubDyn raises an explicit error if asked to read tables we don’t know about, so silently skipping them is fine; we just consume their rows until the next section divider. Extending the parser to cover one of those sections is a localised change — drop another _consume_table call into _parse() and add a typed list field on SubDynFile.

The companion to_pybmodes_pile_tower() adapter takes the parsed SubDyn together with the ElastoDyn main + tower files and synthesises a single combined-cantilever BMIFile + SectionProperties for pyBmodes’ tower modal solver — a “rigid base + flexible pile + flexible tower” model with no soil flexibility (matching the OC3 reference design, which clamps the pile rigidly at the seabed).

class SubDynCircProp(prop_set_id, E, G, rho, D, t)[source]

Bases: object

Circular beam cross-section property set.

Parameters:

prop_set_id (int)
E (float)
G (float)
rho (float)
D (float)
t (float)

D: float

E: float

property EA: float: Axial stiffness, N.

property EI: float: Bending stiffness about either transverse axis, N·m².

G: float

property GJ: float: Torsional stiffness, N·m².

property I: float: Second moment of area, m⁴ (same about either bending axis).

property J: float: Polar moment of area, m⁴.

property area: float: Cross-sectional area of the thin-walled tube, m².

property mass_per_length: float: Distributed mass density, kg/m.

prop_set_id: int

rho: float

t: float

class SubDynFile(header='', title='', source_file=None, joints=<factory>, members=<factory>, circ_props=<factory>, reaction_joint_id=0, interface_joint_id=0)[source]

Bases: object

All parameters parsed from a SubDyn .dat file.

Parameters:

header (str)
title (str)
source_file (Path | None)
joints (list[SubDynJoint])
members (list[SubDynMember])
circ_props (list[SubDynCircProp])
reaction_joint_id (int)
interface_joint_id (int)

circ_props: list[SubDynCircProp]

header: str = ''

interface_joint_id: int = 0

joints: list[SubDynJoint]

members: list[SubDynMember]

reaction_joint_id: int = 0

source_file: Path | None = None

title: str = ''

class SubDynJoint(joint_id, x, y, z, joint_type=1)[source]

Bases: object

One node in the substructure graph (SS-coordinate system, metres).

Parameters:

joint_id (int)
x (float)
y (float)
z (float)
joint_type (int)

joint_id: int

joint_type: int = 1

x: float

y: float

z: float

class SubDynMember(member_id, joint_a, joint_b, prop_set_a, prop_set_b, member_type='1c')[source]

Bases: object

One beam member connecting two joints.

Parameters:

member_id (int)
joint_a (int)
joint_b (int)
prop_set_a (int)
prop_set_b (int)
member_type (str)

joint_a: int

joint_b: int

member_id: int

member_type: str = '1c'

prop_set_a: int

prop_set_b: int

read_subdyn(path)[source]

Parse a SubDyn .dat file. See module docstring for the supported subset of sections.

Parameters:: path (str | Path)
Return type:: SubDynFile

to_pybmodes_pile_tower(main, tower, subdyn, blade=None)[source]

Build a combined pile + tower BMI + SectionProperties for OC3-style rigid-base monopiles.

Layout (axial coordinate z increases upward):

z = z_seabed   --- rigid clamped base (SubDyn reaction joint)
...                pile section, properties from SubDyn members
z = z_TP       --- transition piece (SubDyn interface joint
                    / ElastoDyn ``TowerBsHt``)
...                tower section, properties from ElastoDyn tower
z = z_top      --- tower top, lumped RNA tip mass

The returned BMIFile describes the full beam from z_seabed to z_top as a single cantilever; the SectionProperties array splices the pile and tower property tables at the transition piece with two stations very close together so the FE interpolant captures the cross-section discontinuity correctly.

No soil flexibility, no hydrodynamic added mass — the OC3 design fixes the pile rigidly at the seabed and the user-selected scope here excludes hydro coupling. See the case-study script in cases/nrel5mw_monopile/run.py for context.

Parameters:: subdyn (SubDynFile)

WAMIT / HydroDyn reader

WAMIT v7 output file reader for OpenFAST / HydroDyn floating-platform decks.

Parses .1 (added-mass + radiation-damping) and .hst (hydrostatic restoring) files written by WAMIT, redimensionalises them per the WAMIT v7 convention (ρ·L^k and ρ·g·L^k factors), and returns the 6 × 6 SI matrices a pyBmodes PlatformSupport block consumes.

The reader is read-only and pure-Python; pairs with a thin HydroDynReader that picks WAMITULEN and the PotFile root from an OpenFAST HydroDyn .dat so callers can chain HydroDynReader(...).read_platform_matrices() in a single step.

WAMIT file formats

.1 (added mass / radiation damping) — each line carries:

period  i  j  A(i,j)  [B(i,j)]

period = -1.0 rows are the infinite-frequency added mass A_inf (only the A column is present; no damping at infinite frequency).
period = 0.0 rows are the zero-frequency added mass A_0 (also A only).
period > 0 rows are frequency-dependent A(ω) + B(ω); this reader currently extracts only A_inf and A_0.

Indices i, j are 1-indexed over the six rigid-body DOFs {1: surge, 2: sway, 3: heave, 4: roll, 5: pitch, 6: yaw}. The file is sparse — only non-trivial entries are listed; missing entries are zero.

.hst (hydrostatic restoring) — each line carries:

i  j  C(i,j)

Same 1-indexed convention. Some WAMIT runs write only the upper triangle; others write the full 6 × 6 including explicit zeros.

WAMIT v7 nondimensionalisation

All WAMIT outputs are dimensionless. Redimensionalisation factors depend on the DOF-pair type:

Added mass (.1): trans-trans → ρ·L³, trans-rot → ρ·L⁴, rot-rot → ρ·L⁵.
Hydrostatic stiffness (.hst): trans-trans → ρ·g·L², trans-rot → ρ·g·L³, rot-rot → ρ·g·L⁴.

Here L = WAMITULEN (declared in the HydroDyn .dat), ρ is the sea-water density, and g is gravitational acceleration. The exponent table is captured in WamitReader._dim_factor() as base + n_rot_dofs where base ∈ {3, 2} and n_rot_dofs ∈ {0, 1, 2}.

class HydroDynReader(dat_path)[source]

Bases: object

Minimal reader for the floating-platform section of a HydroDyn .dat.

Surfaces WAMITULEN, PotMod, PotFile, and PtfmRefzt — the four values needed to drive WamitReader. WtrDens and Gravity are NOT in the modern HydroDyn .dat (HydroDyn ≥ v2.03 delegates those to the paired SeaState input file), so the corresponding properties fall back to standard sea-water defaults.

The parser is deliberately loose: it scans each non-blank line for the pattern <value> <label> ... and stores the first occurrence of each label. Matrix continuation rows (e.g. the five extra rows of AddCLin) and wrapped comment text get parsed too but produce harmless _values[<numeric_string>] = ... entries that no caller asks for.

Parameters:: dat_path (str | pathlib.Path)

property gravity: float

Gravitational acceleration (m/s²).

OpenFAST stores Gravity at the top-level .fst file, not in HydroDyn. We fall back to ISO standard gravity (9.80665) so the reader is self-contained.

property pot_file: str: PotFile — WAMIT root path (verbatim, may carry quotes).

property pot_mod: int: PotMod — 0 = no potential-flow model, 1 = WAMIT.

property ptfm_ref_zt: float: PtfmRefzt — vertical offset of the body reference point (m).

read_platform_matrices()[source]

Resolve PotFile alongside this deck and read its WAMIT outputs.

Raises ValueError if PotMod == 0 (no WAMIT data attached).

Return type:: WamitData

property rho_water: float

Water density (kg/m³).

HydroDyn v2.03+ moves WtrDens to the paired SeaState file; if the legacy inline value isn’t present, or carries Fortran’s "DEFAULT" keyword, we fall back to a standard sea-water default of 1025 kg/m³.

property ulen: float: WAMITULEN — WAMIT reference length used for re-dimensionalisation.

class WamitData(A_inf, A_0, C_hst, rho, g, ulen, pot_file_root)[source]

Bases: object

Dimensionalised WAMIT output for a single floating body.

All matrices are in SI units. Added-mass entries are kg for trans-trans, kg·m for trans-rot, and kg·m² for rot-rot. Hydrostatic-stiffness entries are N/m for trans-trans, N (or N·m/m, same thing) for trans-rot, and N·m/rad for rot-rot.

Variables:

A_inf (ndarray, shape (6, 6)) – Infinite-frequency added mass (period = -1 rows of .1).
A_0 (ndarray, shape (6, 6)) – Zero-frequency added mass (period = 0 rows of .1).
C_hst (ndarray, shape (6, 6)) – Hydrostatic restoring matrix (all rows of .hst).
ulen (rho, g,) – Re-dimensionalisation constants applied; stored so downstream code can audit the SI conversion that was used.
pot_file_root (pathlib.Path) – Absolute path of the resolved PotFile (no extension).

Parameters:

A_inf (ndarray)
A_0 (ndarray)
C_hst (ndarray)
rho (float)
g (float)
ulen (float)
pot_file_root (Path)

A_0: ndarray

A_inf: ndarray

C_hst: ndarray

g: float

pot_file_root: Path

rho: float

ulen: float

class WamitReader(pot_file_root, dat_dir, rho=1025.0, g=9.80665, ulen=1.0)[source]

Bases: object

Read a WAMIT v7 .1 / .hst pair and return SI 6 × 6 matrices.

The constructor only normalises the inputs; on-disk reads happen inside read() so callers can catch FileNotFoundError at a single well-defined point.

Parameters:

pot_file_root (str or Path) – PotFile value taken verbatim from the HydroDyn .dat. May carry surrounding double or single quotes, may use backslashes on Windows, and may be relative (resolved against dat_dir) or absolute.
dat_dir (pathlib.Path) – Directory of the HydroDyn .dat whose PotFile is being followed. Used to resolve relative PotFile values.
rho (float) – Sea-water density (kg/m³), gravitational acceleration (m/s²), and WAMIT reference length (WAMITULEN, m). Default to standard sea-water values so simple smoke tests don’t need an explicit HydroDynReader.
g (float) – Sea-water density (kg/m³), gravitational acceleration (m/s²), and WAMIT reference length (WAMITULEN, m). Default to standard sea-water values so simple smoke tests don’t need an explicit HydroDynReader.
ulen (float) – Sea-water density (kg/m³), gravitational acceleration (m/s²), and WAMIT reference length (WAMITULEN, m). Default to standard sea-water values so simple smoke tests don’t need an explicit HydroDynReader.

read()[source]

Parse the .1 and .hst files alongside the resolved root.

Raises:: FileNotFoundError – If either <root>.1 or <root>.hst is missing; the message names the expected absolute path and the verbatim PotFile value that produced it.
Return type:: WamitData

BModes `.out` output

Parser for modal analysis output files (.out).

Output format:

Header / title
“rotating blade frequencies & mode shapes” OR “tower frequencies & mode shapes”
For each mode:
——– Mode No. N (freq = X.XXXE+XX Hz) <column-header line> <blank line> <data rows: span_loc + 5 columns> <blank line>

Column order:

Blade : span_loc flap_disp flap_slope lag_disp lag_slope twist Tower : span_loc ss_disp ss_slope fa_disp fa_slope twist

exception BModeOutParseError(message, file=None, line=None, column=None, context=None)[source]

Bases: ParseError

Raised by read_out() under strict=True when a .out file is not cleanly parseable.

The default (strict=False) parser is intentionally tolerant — it skips malformed rows so a partially-corrupt file still yields whatever modes it can. Validation / cross-solver-comparison workflows that must trust every value should pass strict=True; the message then carries the source file, the 1-based line number, and the mode context so the offending location is unambiguous.

Inherits pybmodes.io.errors.ParseError (which in turn inherits ValueError), so existing except ValueError callers still catch this exception unchanged — the inheritance addition is backward-compatible. New callers can except ParseError to get a typed handler that works across every pybmodes.io.* parser.

Parameters:

message (str)
file (str | None)
line (int | None)
column (int | None)
context (str | None)

Return type:

None

class BModeOutput(title, beam_type, modes, source_file=None)[source]

Bases: object

All mode shapes parsed from a .out file.

Parameters:

title (str)
beam_type (str)
modes (list[ModeShape])
source_file (Path | None)

beam_type: str

frequencies()[source]

Return array of natural frequencies (Hz) in mode order.

Return type:: ndarray

modes: list[ModeShape]

source_file: Path | None = None

title: str

class ModeShape(mode_number, frequency, span_loc, col1, col2, col3, col4, twist, col_names=<factory>)[source]

Bases: object

Mode shape data for a single mode.

Parameters:

mode_number (int)
frequency (float)
span_loc (ndarray)
col1 (ndarray)
col2 (ndarray)
col3 (ndarray)
col4 (ndarray)
twist (ndarray)
col_names (list[str])

col1: ndarray

col2: ndarray

col3: ndarray

col4: ndarray

col_names: list[str]

property fa_disp: ndarray

property fa_slope: ndarray

property flap_disp: ndarray

property flap_slope: ndarray

frequency: float

property lag_disp: ndarray

property lag_slope: ndarray

mode_number: int

span_loc: ndarray

property ss_disp: ndarray

property ss_slope: ndarray

twist: ndarray

read_out(path, *, strict=False)[source]

Parse a .out file and return a BModeOutput.

With strict=False (default, unchanged behaviour) malformed data rows are silently skipped so a partially-corrupt file still yields whatever modes it can. With strict=True the parser raises BModeOutParseError — with the source file, the 1-based line number, and the mode context — on a short data row, an unparseable number, a non-finite (NaN / inf) value, a duplicate mode number, or a file that yields no modes at all. Use strict=True for validation / cross-solver comparison.

Parameters:

path (str | Path)
strict (bool)

Return type:

BModeOutput

WindIO ontology — tubular tower / monopile

Read the structural subset of a WindIO ontology .yaml for a tubular tower or monopile (issue #35).

WindIO describes a tower / monopile as a circular tube via:

components.<component>.outer_shape.outer_diameter.{grid, values}
components.<component>.structure.layers[] — each {material, thickness.{grid, values}} (summed for the wall)
components.<component>.structure.outfitting_factor — the non-structural mass multiplier (internals / flanges / paint)
components.<component>.reference_axis.z.{grid, values} — physical station heights (m); the span = |z[-1] - z[0]|
top-level materials[] — the layer’s material {E, rho, nu}

That is exactly what pybmodes.io.geometry.tubular_section_props() needs, so pybmodes.models.Tower.from_windio() is a thin wrapper.

This module is the tubular (tower / monopile) reader only. A WindIO blade is a composite layup whose beam properties need a PreComp-class thin-wall cross-section reduction — that lives in pybmodes.io.windio_blade (read_windio_blade() / pybmodes.models.RotatingBlade.from_windio()), not here.

Requires the optional [windio] extra (PyYAML); the runtime core stays numpy + scipy only, mirroring the [plots] / [notebook] extras.

class WindIOMonopileTower(section_props, combined_length, el_loc, transition_frac, z_base, z_transition, z_top)[source]

Bases: object

A monopile + tower spliced into a single fixed-bottom cantilever.

The combined cantilever runs from the monopile base (mudline, z_base) up through the transition piece (z_transition, where the monopile meets the tower) to the tower top (z_top). The section-property table carries both segments’ own wall schedules and materials, joined at the transition with a near-coincident station pair so the FE interpolant captures the cross-section step.

Parameters:

section_props (SectionProperties)
combined_length (float)
el_loc (ndarray)
transition_frac (float)
z_base (float)
z_transition (float)
z_top (float)

combined_length: float

el_loc: ndarray

section_props: SectionProperties

transition_frac: float

z_base: float

z_top: float

z_transition: float

class WindIOTubular(station_grid, outer_diameter, wall_thickness, flexible_length, E, rho, nu, outfitting_factor, z_base=0.0, z_top=0.0)[source]

Bases: object

Geometry + material extracted from a WindIO tower / monopile.

Parameters:

station_grid (ndarray)
outer_diameter (ndarray)
wall_thickness (ndarray)
flexible_length (float)
E (float)
rho (float)
nu (float)
outfitting_factor (float)
z_base (float)
z_top (float)

E: float

flexible_length: float

nu: float

outer_diameter: ndarray

outfitting_factor: float

rho: float

station_grid: ndarray

wall_thickness: ndarray

z_base: float = 0.0

z_top: float = 0.0

read_windio_monopile_tower(yaml_path, *, component_tower='tower', component_monopile='monopile', thickness_interp='linear', n_nodes=None)[source]

Reduce the monopile and tower components and splice them into one fixed-bottom cantilever (issue #92).

Each component is reduced independently through the closed-form circular-tube relations (so they keep their own wall schedule and steel grade), then concatenated bottom-to-top at the transition piece — the elevation where the monopile top meets the tower base. The result is the WindIO analog of pybmodes.io.subdyn_reader.to_pybmodes_pile_tower() (the ElastoDyn + SubDyn splice).

Parameters:

yaml_path (path to a WindIO ontology file carrying both a) – monopile and a tower component.
component_tower (component names to splice) – (defaults "tower" / "monopile").
component_monopile (component names to splice) – (defaults "tower" / "monopile").
thickness_interp ("linear" or "piecewise_constant" — passed) – through to each component’s reduction (see read_windio_tubular()).
n_nodes (optional FE-mesh refinement, applied per segment: each) – of the monopile and tower is re-gridded onto n_nodes evenly- spaced stations (geometry interpolated, tube properties recomputed exactly), mirroring Tower.from_windio()’s n_nodes. None keeps each component’s native WindIO grid.

Return type:

WindIOMonopileTower

:raises ValueError : when the monopile top and tower base do not meet at a: common transition-piece elevation (a gap or overlap of more than 1 mm), since a non-contiguous pair cannot be spliced into one beam.

read_windio_tubular(yaml_path, *, component='tower', thickness_interp='linear')[source]

Parse the structural subset of component from a WindIO file.

Handles both WindIO key dialects (modern outer_shape/structure and older outer_shape_bem/internal_structure_2d_fem); see _shape_and_structure().

Parameters:

yaml_path (str | Path)
component (str)
thickness_interp (str)

Return type:

WindIOTubular

WindIO ontology — composite blade

Read a WindIO ontology .yaml blade and reduce it to the FEM section-property table (issue #35).

This is the public glue that ties Phase-2 together, mirroring pybmodes.io.windio (tower) / pybmodes.io.geometry. tubular_section_props():

read_windio_blade() — parse the blade component (dialect-robust, reusing the duplicate-anchor-tolerant loader from pybmodes.io.windio): span axis, chord, twist, reference-axis chordwise location, the spanwise airfoil set, the resolved web / layer nd_arc bands (pybmodes.io._precomp.arc_resolver), and the material table.
windio_blade_section_props() — walk the span, blend the airfoil, build each station’s shell-layer / web stacks, run the thin-wall reduction (pybmodes.io._precomp.reduction), and assemble a pybmodes.io.sec_props.SectionProperties ready for pybmodes.models.RotatingBlade.

Both WindIO key dialects are handled — modern outer_shape / structure (IEA-15 WT_Ontology, every WISDEM example incl. the floating ones) and older outer_shape_bem / internal_structure_2d_fem (IEA-3.4 / 10 / 22). Needs the optional [windio] extra (PyYAML); the runtime core stays numpy+scipy.

class WindIOBlade(span_grid, flexible_length, chord, twist_deg, ref_axis_xc, profiles, resolved, materials, elastic=None, elastic_parse_error=None)[source]

Bases: object

Geometry + layup of a WindIO blade, resolved onto a span grid.

Parameters:

span_grid (ndarray)
flexible_length (float)
chord (ndarray)
twist_deg (ndarray)
ref_axis_xc (ndarray)
profiles (list[Profile])
resolved (ResolvedBladeStructure)
materials (dict)
elastic (dict | None)
elastic_parse_error (str | None)

chord: ndarray

elastic: dict | None = None: Pre-computed distributed beam properties parsed straight from the WindIO elastic_properties / elastic_properties_mb block (the published reference), interpolated onto span_grid; None when the file carries only the layup. Keys: mass_den/flp_iner/edge_iner/flp_stff/ edge_stff/tor_stff/axial_stff/cg_offst.

elastic_parse_error: str | None = None: Non-None when a published block was present but could not be parsed (schema drift / malformed). elastic is then None too, but this distinguishes “absent” (silent PreComp fallback is correct) from “present-but-broken” ("auto" warns; "file" raises) — so a typo can’t hide behind a plausible lower-fidelity result.

flexible_length: float

materials: dict

profiles: list[Profile]

ref_axis_xc: ndarray

resolved: ResolvedBladeStructure

span_grid: ndarray

twist_deg: ndarray

read_windio_blade(yaml_path, *, component='blade', n_span=30)[source]

Parse the structural subset of a WindIO blade component.

Parameters:

yaml_path (str | Path)
component (str)
n_span (int)

Return type:

WindIOBlade

windio_blade_section_props(blade, *, n_perim=300, title='WindIO composite-blade section properties', elastic='auto')[source]

Reduce every span station to the FEM section-property table.

elastic selects the property source (issue #48 — keep deltas to the reference model small):

"auto" (default) — use the WindIO published distributed beam properties (elastic_properties / elastic_properties_mb) when the file carries them, so pyBmodes matches the reference model’s stiffness/inertia exactly; fall back to the PreComp thin-wall reduction of the layup only when they are absent.
"precomp" — always run the PreComp reduction (the pre-1.5 behaviour), even when published properties exist.
"file" — require the published properties; raise ValueError if the file has only the layup or carries a published block that could not be parsed.

If a published block is present but unparseable (schema drift / malformed), "auto" does not silently fall back to the lower-fidelity PreComp result — it emits a UserWarning naming the parse problem before reducing the layup, and "file" raises (issue #47 follow-up — static review).

Parameters:

blade (WindIOBlade)
n_perim (int)
title (str)
elastic (str)

Return type:

SectionProperties

WindIO ontology — floating substructure

Read the WindIO components.floating_platform + mooring blocks (issue #35).

A WindIO floating substructure is a set of named joints (3-D points, MSL datum, z up) connected by slender circular members (each a wall layup + optional bulkhead + ballast, plus Morison Ca/Cd). The downstream physics — hydrostatic restoring, Morison added mass + buoyancy + rigid-body inertia, and the catenary mooring stiffness (pybmodes.mooring) — all build on the geometry parsed here.

This module is only the parser + geometry primitives. Conventions:

cylindrical: true joints give location = [r, theta_deg, z] (WISDEM windIO convention) → x = r cos θ, y = r sin θ.
axial_joints (named fractions along a member, e.g. a fairlead) are resolved into the joint table so mooring can reference them.
the joint flagged transition: true is where the tower foots.

Needs the optional [windio] extra (PyYAML); runtime core stays numpy+scipy. Dialect-robust via the duplicate-anchor-tolerant loader shared with pybmodes.io.windio.

class FloatingMember(name, end1, end2, od_grid, od_values, wall_material, wall_t_grid, wall_t_values, ca=1.0, bulkhead_material=None, bulkhead_t=0.0, ballast=<factory>)[source]

Bases: object

One circular member: end points (m, MSL datum), the spanwise outer-diameter curve, the wall + bulkhead layup, ballast, and the Morison coefficients.

Parameters:

name (str)
end1 (ndarray)
end2 (ndarray)
od_grid (ndarray)
od_values (ndarray)
wall_material (str)
wall_t_grid (ndarray)
wall_t_values (ndarray)
ca (float)
bulkhead_material (str | None)
bulkhead_t (float)
ballast (list)

property axis: ndarray: Unit vector end1 → end2.

ballast: list

bulkhead_material: str | None = None

bulkhead_t: float = 0.0

ca: float = 1.0

diameter_at(frac)[source]: Outer diameter at member fraction(s) frac ∈ [0, 1] (scalar or array).

end1: ndarray

end2: ndarray

property length: float

name: str

od_grid: ndarray

od_values: ndarray

point_at(frac)[source]

3-D location at member fraction frac ∈ [0, 1].

Parameters:: frac (float)
Return type:: ndarray

wall_material: str

wall_t_at(frac)[source]: Wall thickness at member fraction(s) frac (scalar/array).

wall_t_grid: ndarray

wall_t_values: ndarray

class WindIOFloating(members, joints, transition_joint, transition_piece_mass, mooring, materials)[source]

Bases: object

Parsed floating substructure + raw mooring block.

Parameters:

members (list[FloatingMember])
joints (dict)
transition_joint (str | None)
transition_piece_mass (float)
mooring (dict)
materials (dict)

joints: dict

materials: dict

members: list[FloatingMember]

mooring: dict

transition_joint: str | None

transition_piece_mass: float

added_mass(floating, ref_point=None, *, rho=1025.0, z_msl=0.0, n_seg=200, ca_end=0.6)[source]

6×6 infinite-frequency added mass (Morison + member-end caps).

Following RAFT (raft_member): each submerged member element contributes a transverse added mass ρ·Ca·(πD²/4) ⟂ to its axis (M3 = a'(I − n nᵀ)), and each submerged member end contributes an axial end-cap added mass ρ·Ca_End·(2/3)πr³ along the axis (n nᵀ) — the heave-plate / end effect that closes most of the strip-only heave gap (Ca_End default 0.6, RAFT’s default). Both are kinematically transformed to the platform reference. Still a Morison proxy (no radiation diffraction / member interaction), so a documented approximation to a potential-flow A_inf; the WAMIT deck-fallback supplies the exact matrix when present.

Parameters:

floating (WindIOFloating)
ref_point (ndarray | None)
rho (float)
z_msl (float)
n_seg (int)
ca_end (float)

Return type:

hydrostatic_restoring(floating, *, rho=1025.0, g=9.80665, z_msl=0.0, n_seg=200)[source]

6×6 hydrostatic restoring (DOF order surge, sway, heave, roll, pitch, yaw) from member geometry — the WAMIT/.hst buoyancy + waterplane convention (no gravitational term; that enters via the body mass elsewhere). For a freely-floating semi the heave / roll / pitch entries are geometry-exact, so this matches a potential-flow .hst closely (integration anchor).

Parameters:

floating (WindIOFloating)
rho (float)
g (float)
z_msl (float)
n_seg (int)

Return type:

read_windio_floating(yaml_path, *, component='floating_platform')[source]

Parse the floating substructure + mooring from a WindIO file.

Parameters:

yaml_path (str | Path)
component (str)

Return type:

WindIOFloating

rigid_body_inertia(floating, ref_point=None, *, n_seg=200)[source]

Structural rigid-body mass of the substructure about ref_point: returns (mass, M6x6, cg).

Counts the member steel wall (thin-wall ρ·πD·t per length), end bulkheads, fixed ballast (explicit volume × material density) and the transition-piece mass. Variable (trim) ballast is intentionally excluded — it is an equilibrium quantity of the fully-assembled turbine, not derivable from the floating component alone (the validated total is taken from the companion ElastoDyn PtfmMass when available).

Parameters:

floating (WindIOFloating)
ref_point (ndarray | None)
n_seg (int)

Return type:

tuple[float, ndarray, ndarray]

Geometry & section properties

Build pyBmodes section properties from tubular geometry instead of pre-computed structural properties.

Wind-turbine towers and monopiles are circular tubes: given the outer diameter, wall thickness, and material, every structural property the FEM needs is an exact closed-form expression — so the user supplies only what they actually know (geometry) and pyBmodes derives mass / EI / GJ / EA, eliminating the hand-computation error class (issue #35).

For a circular tube of outer radius Ro and inner Ri made of an isotropic material (E, rho, nu):

A   = pi (Ro^2 - Ri^2)                      cross-section area
I   = (pi/4) (Ro^4 - Ri^4)                  area 2nd moment (FA == SS)
J   = 2 I                                   polar 2nd moment (tube)
G   = E / (2 (1 + nu))                      shear modulus

mass_den   = rho * A * outfitting_factor    kg/m   (outfitting lumps
                                             internals / flanges /
                                             paint into the mass)
flp_stff   = edge_stff = E * I              N*m^2
tor_stff   = G * J = E * I / (1 + nu)       N*m^2
axial_stff = E * A                          N
flp_iner   = edge_iner = rho * I            kg*m   (rotary inertia)

These are the same homogeneous-material identities the floating section-property path uses (axial = mass.E/rho, tor = EI/(1+nu), rho.I = EI.rho/E) — here derived forward from geometry rather than back-solved from stiffness. outfitting_factor multiplies only the mass terms (it is non-structural mass), never the stiffness — the same separation the AdjTwMa fix established.

tubular_section_props(span_loc, outer_diameter, wall_thickness, *, E, rho, nu=0.3, outfitting_factor=1.0, title='geometry-derived tubular section properties')[source]

Exact circular-tube section properties for a steel/iso tower.

Parameters:

span_loc ((n,) normalised station locations in [0, 1] (root) – -> tip), strictly the same convention the solver expects.
outer_diameter ((n,) metres, per station.)
wall_thickness ((n,) metres, per station.)
E (isotropic material — Young's modulus (Pa), density) – (kg/m^3), Poisson’s ratio.
rho (isotropic material — Young's modulus (Pa), density) – (kg/m^3), Poisson’s ratio.
nu (isotropic material — Young's modulus (Pa), density) – (kg/m^3), Poisson’s ratio.
outfitting_factor (non-structural mass multiplier (internals,) – flanges, paint, bolts). Multiplies the distributed mass density ONLY. Rotary inertia is treated as a structural section property and is not scaled (it stays rho * i_area); stiffness is never scaled.
title (str)

Returns:

SectionProperties (FA == SS, no twist / offsets — an axisymmetric
tube has none), ready for the FEM pipeline.

Return type:

SectionProperties

Parser for beam section-properties files (.dat).

File structure (all units SI):

Line 1 : title string
Line 2 : n_secs  label  description
Line 3 : blank
Line 4 : column header
Line 5 : column units
Lines 6+: one row per spanwise station (13 space-separated values)
          trailing notes / blank lines after the data are ignored

Column order:

span_loc  str_tw  tw_iner  mass_den  flp_iner  edge_iner
flp_stff  edge_stff  tor_stff  axial_stff  cg_offst  sc_offst  tc_offst

class SectionProperties(title, n_secs, span_loc, str_tw, tw_iner, mass_den, flp_iner, edge_iner, flp_stff, edge_stff, tor_stff, axial_stff, cg_offst, sc_offst, tc_offst, source_file=None)[source]

Bases: object

Spanwise section property table.

Parameters:

title (str)
n_secs (int)
span_loc (ndarray)
str_tw (ndarray)
tw_iner (ndarray)
mass_den (ndarray)
flp_iner (ndarray)
edge_iner (ndarray)
flp_stff (ndarray)
edge_stff (ndarray)
tor_stff (ndarray)
axial_stff (ndarray)
cg_offst (ndarray)
sc_offst (ndarray)
tc_offst (ndarray)
source_file (Path | None)

axial_stff: ndarray

cg_offst: ndarray

edge_iner: ndarray

edge_stff: ndarray

flp_iner: ndarray

flp_stff: ndarray

mass_den: ndarray

n_secs: int

sc_offst: ndarray

source_file: Path | None = None

span_loc: ndarray

str_tw: ndarray

tc_offst: ndarray

title: str

tor_stff: ndarray

tw_iner: ndarray

read_sec_props(path)[source]

Parse a section-properties .dat file.

Parameters:: path (str | Path)
Return type:: SectionProperties

Polynomial fitting + ElastoDyn helpers

Constrained 6th-order polynomial fit for ElastoDyn mode shapes.

ElastoDyn requires mode shapes of the form phi(x) = C2*x^2 + C3*x^3 + C4*x^4 + C5*x^5 + C6*x^6 with constraint phi(1) = 1, i.e. C2+C3+C4+C5+C6 = 1.

The constraint is enforced analytically by substituting C6 = 1 - C2 - C3 - C4 - C5 and solving the reduced least-squares system.

class PolyFitResult(c2, c3, c4, c5, c6, rms_residual, tip_slope, cond_number)[source]

Bases: object

Polynomial fit coefficients and quality metrics for one mode component.

Parameters:

c2 (float)
c3 (float)
c4 (float)
c5 (float)
c6 (float)
rms_residual (float)
tip_slope (float)
cond_number (float)

c2: float

c3: float

c4: float

c5: float

c6: float

coefficients()[source]

Return [C2, C3, C4, C5, C6] as a length-5 array.

Return type:: ndarray

cond_number: float

evaluate(x)[source]

Evaluate the polynomial at positions x in [0, 1].

Parameters:: x (ndarray)
Return type:: ndarray

rms_residual: float

tip_slope: float

fit_mode_shape(span_loc, displacement)[source]

Fit a constrained 6th-order polynomial to a normalised mode shape.

Parameters:

span_loc (1-D array of normalised span positions, x in [0, 1].)
displacement (1-D array of mode-shape displacements at each station.)

Raises:

ValueError – If the inputs are not 1-D, lengths differ, the array is too short to fit five free coefficients, any element is non-finite, span_loc is not strictly increasing, or the tip displacement is effectively zero.

Return type:

PolyFitResult

Map ModalResult + poly fit to named ElastoDyn input parameters.

class BladeElastoDynParams(BldFl1Sh, BldFl2Sh, BldEdgSh)[source]

Bases: object

Polynomial fits for the three blade mode shapes required by ElastoDyn.

Parameters:

BldFl1Sh (PolyFitResult)
BldFl2Sh (PolyFitResult)
BldEdgSh (PolyFitResult)

BldEdgSh: PolyFitResult

BldFl1Sh: PolyFitResult

BldFl2Sh: PolyFitResult

as_dict()[source]

Flat {ElastoDyn_param_name: coefficient_value} for writer.

Return type:: dict[str, float]

class TowerElastoDynParams(TwFAM1Sh, TwFAM2Sh, TwSSM1Sh, TwSSM2Sh)[source]

Bases: object

Polynomial fits for the four tower mode shapes required by ElastoDyn.

Field names match the OpenFAST tower sub-file <turbine>_ElastoDyn_tower.dat:

TwFAM1Sh / TwFAM2Sh — fore-aft modes 1 and 2
TwSSM1Sh / TwSSM2Sh — side-side modes 1 and 2

Parameters:

TwFAM1Sh (PolyFitResult)
TwFAM2Sh (PolyFitResult)
TwSSM1Sh (PolyFitResult)
TwSSM2Sh (PolyFitResult)

TwFAM1Sh: PolyFitResult

TwFAM2Sh: PolyFitResult

TwSSM1Sh: PolyFitResult

TwSSM2Sh: PolyFitResult

as_dict()[source]

Return type:: dict[str, float]

class TowerFamilyMemberReport(mode_number, frequency_hz, family_rank, is_fa, fa_rms, ss_rms, direction_ratio, fit_rms, fit_is_good, selected, fa_participation=0.0, ss_participation=0.0, torsion_participation=0.0, torsion_rejected=False)[source]

Bases: object

Diagnostic view of one scored FA/SS tower family candidate.

fa_participation / ss_participation / torsion_participation are normalised modal kinetic-energy fractions (unit-mass approximation: Σ φ_axis² / Σ φ_total² over all FEM nodes). Sum to 1 for every mode. torsion_rejected is True when the torsion fraction crosses _TORSION_REJECT_THRESHOLD (10 %); such modes are dropped from the family selection because they are no longer pure bending modes ElastoDyn’s polynomial ansatz can faithfully represent.

Parameters:

mode_number (int)
frequency_hz (float)
family_rank (int)
is_fa (bool)
fa_rms (float)
ss_rms (float)
direction_ratio (float)
fit_rms (float)
fit_is_good (bool)
selected (bool)
fa_participation (float)
ss_participation (float)
torsion_participation (float)
torsion_rejected (bool)

direction_ratio: float

fa_participation: float = 0.0

fa_rms: float

family_rank: int

fit_is_good: bool

fit_rms: float

frequency_hz: float

is_fa: bool

mode_number: int

selected: bool

ss_participation: float = 0.0

ss_rms: float

torsion_participation: float = 0.0

torsion_rejected: bool = False

class TowerSelectionReport(fa_family, ss_family, selected_fa_modes, selected_ss_modes, rejected_fa_modes=(), rejected_ss_modes=())[source]

Bases: object

Structured report of tower-mode family scoring and final selection.

rejected_fa_modes / rejected_ss_modes carry the mode numbers of family candidates that were dropped because their torsion-participation crossed _TORSION_REJECT_THRESHOLD. They’re empty when every candidate passed the gate.

Parameters:

fa_family (tuple[TowerFamilyMemberReport, ...])
ss_family (tuple[TowerFamilyMemberReport, ...])
selected_fa_modes (tuple[int, int])
selected_ss_modes (tuple[int, int])
rejected_fa_modes (tuple[int, ...])
rejected_ss_modes (tuple[int, ...])

fa_family: tuple[TowerFamilyMemberReport, ...]

rejected_fa_modes: tuple[int, ...] = ()

rejected_ss_modes: tuple[int, ...] = ()

selected_fa_modes: tuple[int, int]

selected_ss_modes: tuple[int, int]

ss_family: tuple[TowerFamilyMemberReport, ...]

compute_blade_params(modal)[source]

Fit polynomials to the 1st/2nd flap and 1st edge modes.

Parameters:: modal (ModalResult)
Return type:: BladeElastoDynParams

compute_tower_params(modal)[source]

Fit polynomials to the 1st/2nd FA and 1st/2nd SS tower modes.

Parameters:: modal (ModalResult)
Return type:: TowerElastoDynParams

compute_tower_params_report(modal)[source]

Compute tower ElastoDyn parameters and return a selection report.

Pre-rotates any degenerate FA/SS eigenpair (consecutive modes whose relative frequency gap is below _DEGENERATE_FREQ_RTOL) into clean direction-aligned shapes before classification. The original modal.shapes is not mutated; the rotation only feeds the candidate builder used here. See _rotate_degenerate_pairs().

Emits a RuntimeWarning if the polynomial-fit design-matrix condition number exceeds _FIT_COND_WARN (1e4) — the fit may be unreliable because the reduced Vandermonde basis is poorly conditioned at the mesh stations supplied. Emits a stronger warning above _FIT_COND_FAIL (1e6); the reconstructed shape may still be accurate but individual C2..C6 coefficients are likely to swing wildly under small perturbations of the input data.

Parameters:: modal (ModalResult)
Return type:: tuple[TowerElastoDynParams, TowerSelectionReport]

Coefficient-consistency validation for OpenFAST ElastoDyn decks.

The polynomial coefficient blocks shipped in industry ElastoDyn .dat files (NREL 5MW Reference Turbine, IEA-3.4-130-RWT, and others) are demonstrably inconsistent with the structural-property blocks in the same files — see cases/ECOSYSTEM_FINDING.md for the longer write-up. This module exposes a programmatic way to surface that inconsistency:

validate_dat_coefficients() parses the deck, runs pyBmodes on the same structural inputs, fits its own polynomials, and computes per-block RMS residuals for both the file’s polynomial and pyBmodes’ polynomial against the FEM-computed mode shape.
ValidationResult and CoeffBlockResult carry the results in a form the CLI (pybmodes validate) and the Tower.from_elastodyn(validate_coeffs=True) path both consume.

class CoeffBlockResult(name, file_rms, pybmodes_rms, ratio, verdict, file_coeffs, pybmodes_coeffs, fa_participation=nan, ss_participation=nan, torsion_participation=nan, rejected_modes=<factory>)[source]

Bases: object

Validation result for a single ElastoDyn coefficient block.

Both file_rms and pybmodes_rms are RMS residuals of the respective polynomial evaluated at the FEM stations against the same tip-normalised mode-shape data. They are directly comparable. ratio = file_rms / pybmodes_rms quantifies how much worse the file polynomial fits the structural model than pyBmodes’ fit does; a ratio above ~50× indicates the file polynomial was generated against a different structural model than the one in the deck.

Tower blocks (TwFAM1Sh / TwFAM2Sh / TwSSM1Sh / TwSSM2Sh) carry four extra fields populated from the family- selection report:

fa_participation / ss_participation / torsion_participation — modal kinetic-energy fractions (unit-mass approximation) of the selected mode for this block. Sum to 1.
rejected_modes — mode numbers that were dropped from this block’s family during selection because their torsion fraction exceeded 10 %. Empty when every candidate passed the gate.

Blade blocks leave all four fields at their NaN / empty defaults because the blade adapter doesn’t run a torsion-contamination filter (blade torsion is part of the structural model and not handled the same way as tower torsion).

Parameters:

name (str)
file_rms (float)
pybmodes_rms (float)
ratio (float)
verdict (Literal['PASS', 'WARN', 'FAIL'])
file_coeffs (list[float])
pybmodes_coeffs (list[float])
fa_participation (float)
ss_participation (float)
torsion_participation (float)
rejected_modes (tuple[int, ...])

fa_participation: float = nan

file_coeffs: list[float]

file_rms: float

name: str

pybmodes_coeffs: list[float]

pybmodes_rms: float

ratio: float

rejected_modes: tuple[int, ...]

ss_participation: float = nan

torsion_participation: float = nan

verdict: Literal['PASS', 'WARN', 'FAIL']

class ValidationResult(dat_path, tower_results=<factory>, blade_results=<factory>, overall='PASS', summary='')[source]

Bases: object

Aggregated coefficient-validation result for an ElastoDyn deck.

Parameters:

dat_path (Path)
tower_results (dict[str, CoeffBlockResult])
blade_results (dict[str, CoeffBlockResult])
overall (Literal['PASS', 'WARN', 'FAIL'])
summary (str)

all_blocks()[source]

Iterate tower-then-blade blocks in canonical order.

Return type:: dict[str, CoeffBlockResult]

blade_results: dict[str, CoeffBlockResult]

dat_path: Path

failing_blocks()[source]

Return type:: list[CoeffBlockResult]

overall: Literal['PASS', 'WARN', 'FAIL'] = 'PASS'

summary: str = ''

tower_results: dict[str, CoeffBlockResult]

warning_blocks()[source]

Return type:: list[CoeffBlockResult]

validate_dat_coefficients(dat_path, *, verbose=False, n_modes=10)[source]

Validate the polynomial coefficient blocks in an ElastoDyn deck.

Parses the main .dat file plus the tower and blade files it references, builds pyBmodes Tower and RotatingBlade models from the structural-property blocks (NOT the polynomial blocks), runs the eigensolver, fits 6th-order polynomials, and compares those fits against the polynomial coefficients embedded in the deck.

Parameters:

dat_path (Path | str) – Path to the ElastoDyn main .dat file.
verbose (bool) – Reserved for future diagnostic output. The function currently emits no prints regardless of this flag — the CLI front-end is responsible for stdout formatting.
n_modes (int) – Number of FEM modes to extract per model (default 10; must be large enough to cover the four tower bending modes plus the three blade bending modes — pyBmodes’ family selectors warn if any required mode falls outside the requested range).

Returns:

Carries seven CoeffBlockResult entries (4 tower + 3 blade) and an overall PASS/WARN/FAIL verdict (worst across all blocks).

Return type:

ValidationResult

Patch ElastoDyn .dat files with computed polynomial coefficients.

Each coefficient occupies a separate line of the form <value> <ParamName(k)> - comment text. The writer finds each line by matching the parameter name token, then replaces the leading value in-place, preserving indentation and all trailing text.

patch_dat(path, params)[source]

Patch named mode-shape coefficient lines in an ElastoDyn .dat file.

Each parameter in params is located by searching for its exact name as a whitespace-delimited token. The value on the same line (the token before the name) is replaced with the computed coefficient. All other content (indentation, comment text) is left unchanged.

Parameters:

path (path to the ElastoDyn .dat file (modified in place).)
params (BladeElastoDynParams or TowerElastoDynParams.)

Raises:

KeyError if a required parameter name is not found in the file. –

Return type:

None

Floating-tower frequency-gap diagnostic.

A floating ElastoDyn deck has two natural tower bending frequencies that differ by design. The polynomial coefficient blocks in the deck encode the cantilever tower mode shape, which is the only basis ElastoDyn’s SHP = sum_i c_i * (h/H)^(i+1) ansatz can represent (the source-code citations live in src/pybmodes/_examples/reference_decks/FLOATING_CASES.md). The coupled-system frequency that OpenFAST linearisation reports includes platform 6-DOF participation, mooring restoring, and hydrostatic restoring. The two can differ by 20-30 percent on floating platforms, and the gap is expected rather than a bug.

report_floating_frequency_gap runs both pyBmodes solves on the same deck and reports the gap as a short text block, so users reconciling pyBmodes-generated polynomials against OpenFAST linearisation output do not have to re-derive the cantilever-vs-coupled architecture from scratch.

class FloatingFrequencyGap(cantilever_fa_1, cantilever_ss_1, coupled_fa_1, coupled_ss_1)[source]

Bases: object

Cantilever vs coupled tower bending frequencies for one deck.

Cantilever frequencies come from a clamped-base from_elastodyn() solve. This is the modal basis ElastoDyn integrates into MTFA/KTFA at runtime and the basis the shipped polynomial coefficients describe.

Coupled frequencies come from a free-base from_elastodyn_with_mooring() solve with mooring stiffness, hydrostatic restoring, and platform inertia engaged. These are the numbers an OpenFAST linearisation reports.

Parameters:

cantilever_fa_1 (float)
cantilever_ss_1 (float)
coupled_fa_1 (float)
coupled_ss_1 (float)

cantilever_fa_1: float

cantilever_ss_1: float

coupled_fa_1: float

coupled_ss_1: float

format_report()[source]

Return a short text block summarising the gap.

Return type:: str

property gap_fa_1_pct: float

property gap_ss_1_pct: float

report_floating_frequency_gap(main_dat_path, moordyn_dat_path, hydrodyn_dat_path=None, *, n_modes=10)[source]

Run cantilever and coupled solves on the same floating deck.

The cantilever solve is the modal basis ElastoDyn consumes for its runtime tower-bending DOFs. The coupled solve is what OpenFAST linearisation produces when platform 6-DOF, mooring, and hydrostatic restoring are all engaged. The returned FloatingFrequencyGap lets a user reconcile pyBmodes-generated polynomial coefficients against OpenFAST linearisation output without re-deriving the architectural reason they differ.

On the coupled solve, n_modes + 6 modes are requested so that after the six platform rigid-body modes (surge / sway / heave / roll / pitch / yaw) are filtered out, the tower-bending family selector still sees n_modes candidates.

Parameters:

main_dat_path (str | Path)
moordyn_dat_path (str | Path)
hydrodyn_dat_path (str | Path | None)
n_modes (int)

Return type:

FloatingFrequencyGap

Campbell

Campbell-diagram support: rotor-speed sweep with MAC-tracked blade modes and constant-frequency tower modes overlaid on the same plot.

A Campbell diagram plots a turbine’s natural frequencies against rotor speed and overlays the per-revolution excitation lines (1P, 2P, 3P, …); crossings between excitation lines and structural-mode lines flag resonance risks. For a wind-turbine blade the centrifugal-stiffening contribution to the FEM stiffness matrix raises flap-dominated frequencies markedly with rotor speed while edgewise (lag-dominated) modes barely move. The tower lives in an Earth-fixed frame, so its fore-aft / side-to-side bending frequencies don’t depend on rotor speed at all and show up as horizontal lines on the diagram. The NREL 5MW turbine’s canonical resonance call-out — 3P crossing the 1st tower fore-aft mode near ~6.4 rpm — sits right where the cut-in operating envelope begins, which is exactly the kind of constraint this diagram is designed to surface.

Module layout

Phase 3 PR C1 of the v1.x architecture refactor split this from a single 1301-line module into a sub-package. The public API is unchanged; internal helpers live in private sub-modules so each file covers one concern and stays under a few hundred lines:

pybmodes.campbell.result — CampbellResult dataclass plus its NPZ / CSV round-trip.
_models — input dispatcher: path-vs-loaded model, .dat vs. .bmi, optional tower_input keyword.
_classify — mode-naming heuristics (1st flap / 1st tower FA / platform DOFs).
_mac — MAC matrix helpers + Hungarian assignment used by the blade-sweep tracker.
_sweep — rotor-speed sweep drivers and the public campbell_sweep() entry point.
_plot — plot_campbell().

Public API

campbell_sweep() — given an OpenFAST ElastoDyn main .dat, loads the blade and tower from the same deck, sweeps the blade across omega_rpm (with MAC-based mode tracking), solves the tower once, and packs both into a single CampbellResult. .bmi inputs are also accepted and route to blade-only or tower-only sweeps based on beam_type; an explicit tower_input=... keyword adds a tower file alongside a blade .bmi.
plot_campbell() — renders the result with blade modes as solid coloured lines, tower modes as horizontal dashed dark-grey lines, and the per-rev excitation family as light grey rays from the origin. Optional vertical marker at the rated rotor speed.

Defaults are deliberately spare (n_blade_modes=4, n_tower_modes=4) so the diagram shows the modes that actually drive resonance design — 1st/2nd flap, 1st/2nd edge, 1st/2nd tower FA, 1st/2nd tower SS — without crowding the plot with high-order modes that the per-rev family doesn’t reach inside any realistic operating envelope.

class CampbellResult(omega_rpm, frequencies, labels, participation, n_blade_modes, n_tower_modes, mac_to_previous=<factory>)[source]

Bases: object

Frequencies and labels from a Campbell sweep — blade + tower combined.

Variables:

omega_rpm ((N,) array of rotor speeds in rpm.)
frequencies ((N, n_total_modes) array of natural frequencies in Hz.) – Columns are ordered blade modes first, then tower modes. With MAC tracking enabled, blade columns hold the same physical mode across all rotor speeds. Tower columns are constant across rows (tower frequencies don’t depend on rotor speed).
labels (list of length n_total_modes with human-readable mode) – names — blade modes look like "1st flap" / "2nd edge", tower modes are prefixed with "tower" (e.g. "1st tower FA", "1st tower SS") so callers can split the two by string match if needed.
participation ((N, n_total_modes, 3) array of energy fractions in) – the FEM’s per-mode (flap or FA, edge or SS, torsion) axes. Each row sums to 1. Note the axis interpretation is beam-type-specific: blade columns use flap/edge/torsion, tower columns use FA/SS/torsion.
mac_to_previous ((N, n_total_modes) array of per-step MAC values) – between each output slot’s mode shape at step k and the same slot at step k - 1 (i.e. the tracking confidence). Row 0 is filled with NaN (no previous step). Tower columns are also NaN (tower modes don’t change with rotor speed, so a MAC confidence is not meaningful for them).
n_blade_modes (how many of the leading columns are blade modes.)
n_tower_modes (how many of the trailing columns are tower modes.)

Parameters:

omega_rpm (ndarray)
frequencies (ndarray)
labels (list[str])
participation (ndarray)
n_blade_modes (int)
n_tower_modes (int)
mac_to_previous (ndarray)

frequencies: ndarray

labels: list[str]

classmethod load(path, *, allow_legacy_pickle=False)[source]

Read a sweep result back from a .npz archive saved by save().

A legacy pre-1.0 archive whose __meta__ is a pickled object array is refused by default (object-array unpickling can execute arbitrary code); pass allow_legacy_pickle=True to opt in for a file you trust.

Parameters:

path (str | Path)
allow_legacy_pickle (bool)

Return type:

CampbellResult

mac_to_previous: ndarray

n_blade_modes: int

n_tower_modes: int

omega_rpm: ndarray

participation: ndarray

save(path, *, source_file=None)[source]

Write the sweep result to a .npz archive.

Arrays go in as named keys; labels and the two integer scalars ride in via the embedded JSON __meta__ blob alongside the standard pyBmodes-version / timestamp / source-file / git-hash metadata captured by pybmodes.io._serialize._capture_metadata().

Parameters:

path (str | Path)
source_file (str | Path | None)

Return type:

None

to_csv(path)[source]

Write a spreadsheet-friendly CSV with one row per rotor-speed step.

Columns: rpm, then one frequency column per mode (named by the mode’s label), then one MAC-confidence column per mode suffixed with _mac. Tower-mode MAC columns are NaN throughout because tower modes don’t change with rotor speed — kept as columns for shape-stability across blade-only / tower- only / mixed sweeps.

Parameters:: path (str | Path)
Return type:: None

campbell_sweep(input_path, omega_rpm, n_blade_modes=4, n_tower_modes=4, *, tower_input=None, track_by_mac=True)[source]

Build a Campbell-diagram dataset for the given turbine.

Parameters:

input_path (str | pathlib.Path | RotatingBlade | Tower) –
Either a path or an already-loaded model (issue #51):
- an OpenFAST ElastoDyn main .dat file — the function loads the blade and the tower from the deck and runs both;
- a blade .bmi (beam_type = 1) — blade-only sweep unless tower_input is also supplied;
- a tower .bmi (beam_type = 2) — tower-only result (frequencies are constant across omega_rpm; the result is mostly useful for overlay against the per-rev family);
- an already-constructed RotatingBlade or Tower (from any constructor — __init__, from_elastodyn, from_windio, from_windio_floating, …). The model is used verbatim, with no disk re-read, so a single load point feeds both .run() and the sweep, and a from_windio / from_elastodyn model (whose section properties no path can re-read) can finally be swept. Routed to blade/tower by its beam_type so either may be passed here.
omega_rpm (np.ndarray) – 1-D array of rotor speeds in rpm. Ω = 0 is fine and produces the parked-rotor frequencies.
n_blade_modes (int) – Number of blade modes to extract per speed and report in frequencies[:, :n_blade_modes]. Default 4 covers 1st/2nd flap and 1st/2nd edge — the modes that actually drive resonance design. Pushing this much higher just adds high-order flap modes that no realistic per-rev family crosses inside the operating envelope; raise it deliberately when you need them.
n_tower_modes (int) – Number of tower modes (default 4 — 1st/2nd FA + 1st/2nd SS). Drop to 2 to overlay only the 1st FA + 1st SS pair, or push higher for offshore decks where 3rd-mode crossings matter. Ignored when no tower model is available.
tower_input (str | pathlib.Path | Tower | None) – Optional explicit tower — a tower .bmi path or a loaded :class:`~pybmodes.models.Tower` (keyword-only). Useful when input_path is a blade-only deck or a loaded RotatingBlade. Overrides the deck-supplied tower if input_path was an ElastoDyn .dat. So the single-load-point form is campbell_sweep(blade, omega, tower_input=tower).
track_by_mac (bool) – Whether to use MAC across consecutive rotor speeds to keep each blade output column corresponding to the same physical mode. False returns the eigensolver’s native order (useful for debugging mode re-ordering issues). Tower modes don’t change with rotor speed and are unaffected by this flag.

Return type:

CampbellResult.

plot_campbell(result, excitation_orders=None, rated_rpm=None, ax=None, platform_modes=None, log_freq=False, *, operating_rpm=None, freq_max=None)[source]

Render a Campbell diagram from a CampbellResult.

Engineering-report style (issue #54): structural modes are coloured by family, the legend carries only those four family keys, mode names are written inline next to their lines, and the per-rev family is the only thing the legend enumerates as a group:

Blades — green; per-blade curves, name written at the line.
Tower — black; horizontal lines, name at the line.
Platform — red; floating-platform rigid-body modes (surge / sway / heave / roll / pitch / yaw), near-degenerate symmetric pairs merged (surge/sway) to keep the figure clean.
Blade Passing — blue; the per-rev rays (default 1P / 3P / 6P / 9P), each tagged ↑ nP inline (no legend clutter).

operating_rpm=(lo, hi) shades the operating rotor-speed window grey (outside it stays white) and draws a ↔ Operating Speed Range marker.

Parameters:

result (CampbellResult) – Output of campbell_sweep().
excitation_orders (list[int] | None) – Per-rev orders. Default [1, 3, 6, 9].
rated_rpm (float | None) – If given, a thin reference line at the operating rotor speed.
ax (matplotlib.axes.Axes | None) – Existing Axes to draw into; a fresh figure if None.
platform_modes (list[tuple[str, float]] | None) – Optional [(dof, freq_hz), ...] floating rigid-body modes for the screening path (when the result has no platform columns). Drawn in the red Platform family. The coupled-tower path classifies these natively — see campbell_sweep().
log_freq (bool) – Log-scaled frequency axis (the per-rev rays are densely sampled so they render correctly). Default False.
operating_rpm (tuple[float, float] | None) – (lo, hi) rotor-speed operating window (rpm) — shaded grey with an Operating Speed Range marker. None (default) draws no band — backward compatible.
freq_max (float | None) – Upper frequency-axis limit (Hz). None (default) auto-caps the axis just above the highest structural mode so the modes of interest fill the figure (the steep per-rev rays run off the top, as in a standard Campbell report) instead of the axis stretching to the highest ray. Ignored when log_freq=True.
jitter (Note on blade-line)
-------------------------
shows (For ElastoDyn-derived blade FEMs the 1st-flap line typically)
The ()
axial (BMI adapter floors rotary inertia and forces near-rigid)
1e6) (behaviour (EA / EI ≈)
matrices (leaving the dense FEM)
1e11) (ill-conditioned (κ(M) ≈)
subset (which makes LAPACK's)
the (eigenvalue routines wobble on the lowest mode even when)
tracker (underlying eigenvector is identical step to step. The MAC)
flap-dominant (catches this — the participation array stays > 98 %)
correct (in the 1st-flap slot — so the mode identity is)
only
is (the eigenvalue precision suffers. Centrifugal stiffening)
endpoint-to-endpoint (monotonic in physics (Wright 1982);)
reliable (comparisons (parked vs rated) are)
individual-step
not. (monotonicity is)

Return type:

matplotlib.figure.Figure for the rendered axes.

Modal Assurance Criterion

Modal Assurance Criterion (MAC) utilities for comparing mode shapes.

Two public entry points:

mac_matrix() — pairwise MAC between two lists of mode shapes.
compare_modes() — full mode-by-mode comparison between two ModalResult records: MAC matrix, Hungarian-optimal mode pairing, per-pair frequency shift in percent, and the source labels for display.

The MAC formula:

MAC_ij = |φ_i · φ_j|² / ((φ_i · φ_i) (φ_j · φ_j))

where each φ is the concatenated FEM displacement vector across flap / lag / twist axes (3 × n_nodes entries). The metric is sign- and amplitude-invariant; the only quantity that survives is the directional alignment of the two shapes.

Typical use cases:

compare_modes(baseline_result, patched_result) — confirm the polynomial-coefficient patch didn’t actually change the underlying mode shapes (MAC diagonal stays at ~ 1.0; frequency shifts are the only visible delta).
compare_modes(land_result, monopile_result) — quantify the boundary-condition effect on a tower’s mode shapes; the MAC diagonal drops as the lower modes pick up rigid-body contributions from the soft monopile base.

plot_mac() is a lightweight matplotlib helper that renders the MAC matrix as a heatmap with the Hungarian-paired cells highlighted. matplotlib is imported lazily so the rest of this module works without the [plots] extra installed.

class ModeComparison(mac, frequency_shift, paired_modes, label_A='A', label_B='B', freqs_A=<factory>, freqs_B=<factory>)[source]

Bases: object

Result of a compare_modes() call.

Variables:

mac ((n, m) MAC matrix between result_A.shapes and) – result_B.shapes.
frequency_shift ((n_pairs,) percent change in frequency from) – result_A to result_B for each Hungarian-paired mode. Positive values mean result_B is higher in frequency. NaN for pairs where either side’s frequency is non-positive (e.g. rigid-body modes).
paired_modes (list of (i, j) tuples mapping) – result_A.shapes[i] → result_B.shapes[j] under the Hungarian-optimal MAC assignment. Length is min(n, m).
label_B (label_A /) – by plot_mac as axis titles.
freqs_B (freqs_A /) – the comparison object is self-contained for downstream reporting.

Parameters:

mac (ndarray)
frequency_shift (ndarray)
paired_modes (list[tuple[int, int]])
label_A (str)
label_B (str)
freqs_A (ndarray)
freqs_B (ndarray)

freqs_A: ndarray

freqs_B: ndarray

frequency_shift: ndarray

label_A: str = 'A'

label_B: str = 'B'

mac: ndarray

paired_modes: list[tuple[int, int]]

compare_modes(result_A, result_B, *, label_A='baseline', label_B='modified')[source]

Compare two ModalResult records mode-by-mode.

Builds the full MAC matrix, finds the Hungarian-optimal pairing (each result_A mode mapped to one result_B mode), and computes the per-pair frequency shift (f_B - f_A) / f_A in percent.

Pairs are returned in result_A mode-index order (i.e. paired_modes[0] is (0, j₀) for whatever j₀ matched the first result_A mode).

Parameters:

result_A (ModalResult)
result_B (ModalResult)
label_A (str)
label_B (str)

Return type:

ModeComparison

mac_matrix(shapes_A, shapes_B)[source]

Compute the pairwise MAC matrix between two mode-shape lists.

Returns an (n, m) ndarray where out[i, j] is the MAC value between shapes_A[i] and shapes_B[j]. Values are in [0, 1]; 0 means perfectly orthogonal, 1 means perfectly aligned (or anti-aligned — MAC squares the inner product).

Empty inputs are accepted and return a correctly-shaped zero- dimensional ndarray. Zero-norm shapes (i.e. all-zero displacements) get MAC = 0 against every counterpart, since the denominator is undefined and the closest reasonable answer is “no correlation”.

Parameters:

shapes_A (list[NodeModeShape])
shapes_B (list[NodeModeShape])

Return type:

plot_mac(comparison, ax=None, *, annotate=True, cmap='viridis')[source]

Render the MAC matrix as an annotated heatmap.

Hungarian-paired cells get a red outline; other cells are plain. Cell colour ramps from 0 (dark) to 1 (light) via the supplied cmap. Set annotate=False to drop the numerical overlay (useful for large matrices where labels collide).

Returns the parent Figure so the caller can save / show. Raises ImportError if matplotlib isn’t installed (the optional [plots] extra is required).

Parameters:

comparison (ModeComparison)
ax (Axes | None)
annotate (bool)
cmap (str)

Return type:

Figure

shape_to_vector(shape)[source]

Flatten a NodeModeShape into one (3·n_nodes,) vector.

The concatenation order is [flap_disp, lag_disp, twist] — matches the convention used internally by the Campbell tracker. Slope arrays are not included; the MAC contract operates on displacement directions only.

Parameters:: shape (NodeModeShape)
Return type:: ndarray

Pre-solve sanity checks

Pre-solve sanity checks for Tower and RotatingBlade models.

The check_model() entry point runs a small suite of cheap, deterministic checks on a parsed model and returns a list of ModelWarning records describing any anomalies it found. The list is empty when the model is clean.

The checks run automatically inside Tower.run() and RotatingBlade.run() with check_model=True (the default). WARN- and ERROR-severity findings are routed through Python’s warnings module so they surface at the call site without changing the function’s return type. INFO-severity findings are surfaced only when check_model() is called directly — they’re useful context but not actionable noise on every solve.

Suppress the auto-run with Tower(...).run(n_modes, check_model=False) (symmetric for RotatingBlade.run). Suppression is meant for scripted callers that have already validated their input and want the solver path to stay quiet.

Checks performed (see check_model() for the details):

Every section-property field is finite (no NaN, no ±Inf). Runs first so the per-field checks below don’t have to be NaN-aware (NaN silently passes every <= / > / ratio comparison).
Span stations are strictly increasing.
Mass density is strictly positive at every station.
Bending stiffness (FA + SS) does not jump by more than 5× between adjacent stations.
EI_FA / EI_SS ratio stays within [0.1, 10] at every station.
Tower-top RNA mass is not larger than the integrated tower mass.
PlatformSupport 6×6 inertia / hydro / mooring matrices are well- formed: shape (6, 6), all entries finite, and symmetric (within 1e-6 · max|A|). Rank deficiency is not flagged — surge / sway / yaw hydrostatic restoring is legitimately zero on most floaters and mooring layouts can be low-rank by design.
The horizontal platform CM offset (cm_pform_x / cm_pform_y) does not exceed the platform’s yaw radius of gyration √(I_yaw / m) — a larger value is almost always a coordinate- origin offset leaking into a field that means “CM offset from the tower axis”, which mislabels the rigid-body modes (issue #95).
Floating platform inertia is physical: positive mass_pform and strictly-positive i_matrix diagonal (ERROR otherwise).
Floating model carries hydrodynamic added mass (hydro_M not all zero) — omitting it biases every rigid-body frequency high (WARN).
Floating model has some restoring (hydro_K or mooring_K non-zero); with neither, the rigid-body modes collapse to ~0 Hz (WARN).
The requested n_modes does not exceed the model’s DOF count.
The polynomial-fit design matrix on the mesh stations is not ill-conditioned (cond > 1e4 ⇒ WARN, > 1e6 ⇒ ERROR).

Checks 7–10 are the “floating-model readiness” gates (hub_conn = 2 only): they catch the seakeeping omissions a non-specialist makes when a WindIO .yaml (geometry + material only) is treated as sufficient for a floating system — it is not, the way it is for a land tower (issue #95).

exception ModelValidationError(findings)[source]

Bases: ValueError

Raised by .run() when the pre-solve checks find ERROR-severity, non-physical input and on_error="raise" (the default).

Inherits ValueError so existing except ValueError callers keep catching it. The offending findings are on findings.

Variables:: findings (list of ModelWarning) – The ERROR-severity findings that triggered the fail-closed path. WARN / INFO findings are not collected here (those never raise).
Parameters:: findings (list[ModelWarning])

class ModelWarning(severity, message, location)[source]

Bases: object

One finding from check_model().

Variables:

severity ("INFO", "WARN", or "ERROR". INFO is) – contextual (e.g. “RNA mass dominates the structure”); WARN indicates the solve will probably complete but with degraded accuracy; ERROR indicates a non-physical input that will produce undefined results.
message (human-readable description of the finding.)
location (dotted path to the offending data, e.g.) – "section_properties.mass_den" or "bmi.support.hydro_K".

Parameters:

severity (Literal['INFO', 'WARN', 'ERROR'])
message (str)
location (str)

location: str

message: str

severity: Literal['INFO', 'WARN', 'ERROR']

apply_findings(model, *, n_modes, on_error='raise', stacklevel=2)[source]

Run check_model() and route the findings for a .run() call.

Shared by Tower.run() and RotatingBlade.run(). INFO findings are dropped (contextual, not actionable on the solve path). WARN findings always go through UserWarning. ERROR findings are non-physical input, so they fail closed by default (on_error="raise" raises ModelValidationError); pass on_error="warn" to downgrade them to warnings and continue, the pre-1.14.0 behaviour.

Parameters:

model (the Tower / RotatingBlade being solved.)
n_modes (forwarded to check_model() for the DOF-count gate.)
on_error ("raise" (default, fail closed on ERROR) or "warn".)
stacklevel (forwarded to warnings.warn() so the warning points) – at the user’s .run() call site.

Return type:

None

check_model(model, *, n_modes=None)[source]

Run the full pre-solve check suite on model.

Parameters:

model (a Tower or RotatingBlade instance.)
n_modes (optional. When supplied, additionally checks that) – n_modes doesn’t exceed the model’s DOF count. Skipped when None (e.g. when callers want to validate the static model definition before deciding how many modes to request).

Returns:

Empty when every check passes. The list is ordered roughly by check number (see module docstring), which keeps the output diffable across runs.

Return type:

list of ModelWarning

Numerical options

Centralised numerical-options dataclasses. Defaults preserve every previously-embedded magic number from fem/solver.py, elastodyn/params.py, and checks.py — importing this module changes no behaviour; the value is that callers (and reviewers) can now find every numerical threshold in one place.

Centralised numerical options for pyBmodes.

Three frozen dataclasses gather the magic numbers that used to live as module-level constants scattered across the FEM solver, the polynomial fitter, and the pre-solve sanity checker. Each carries the same default the original constant did, so importing this module changes no behaviour; the value is that callers (and reviewers) can now find every numerical threshold in one place.

The dataclasses are dataclasses.dataclass(frozen=True)() so they hash, compare structurally, and can’t be mutated after construction — suitable as defaults on public APIs and as keys in caches.

Examples

Read the default thresholds:

from pybmodes.options import SolverOptions, FitOptions, CheckOptions

SolverOptions().sparse_ndof_threshold     # 500
FitOptions().polynomial_rms_threshold     # 0.09
CheckOptions().stiffness_jump_factor      # 5.0

Override one threshold while keeping the rest at default:

custom = SolverOptions(sparse_ndof_threshold=2000)

The dataclasses are intentionally additive: future fields gain sensible defaults so existing call sites stay working. Removing or renaming a field is a semver-major change (see API contract).

class CheckOptions(stiffness_jump_factor=5.0, ei_ratio_min=0.1, ei_ratio_max=10.0, support_asymmetry_rtol=1e-06, fit_cond_warn=10000.0, fit_cond_fail=1000000.0, platform_cm_offset_gyradius_factor=1.0)[source]

Bases: object

Pre-solve sanity-check thresholds (consumed by pybmodes.checks.check_model()).

Variables:

stiffness_jump_factor (float, default 5.0) – Maximum allowed ratio of bending stiffness between adjacent section nodes (forward and backward). A jump above this flags a likely transition-piece discontinuity that the polynomial fit will struggle to represent — typically wants extra mesh refinement around it.
ei_ratio_min (float, default 0.1) – Lower bound on EI_FA / EI_SS at every section node. Below this is unphysical for a real tower section; the structural inputs have probably swapped FA / SS.
ei_ratio_max (float, default 10.0) – Upper bound on EI_FA / EI_SS. Symmetric to ei_ratio_min.
support_asymmetry_rtol (float, default 1e-6) – Tolerance (relative to max|support|) for treating a PlatformSupport 6 x 6 matrix as symmetric. Above this, the check warns that the support matrix is asymmetric, which is load-bearing for OC3-Hywind-style cross-coupled platforms but commonly indicates a deck-assembly bug for axisymmetric ones.
fit_cond_warn (float, default 1e4) – Polynomial-fit design-matrix condition above which check_model emits a WARN. Mirrors FitOptions.fit_cond_warn so the check can run independently of the fit module.
fit_cond_fail (float, default 1e6) – Polynomial-fit design-matrix condition above which check_model emits a FAIL. Mirrors FitOptions.fit_cond_fail.
platform_cm_offset_gyradius_factor (float, default 1.0) – The horizontal platform CM offset (cm_pform_x / cm_pform_y on a PlatformSupport) is flagged when its magnitude exceeds this factor times the platform’s yaw radius of gyration √(I_yaw / m). cm_pform_x / cm_pform_y are the CM offset from the tower axis; a value comparable to the platform’s own size is almost always a coordinate-origin error leaking into the field, which injects spurious surge/sway↔yaw coupling and mislabels the rigid-body modes (issue #95).

Parameters:

stiffness_jump_factor (float)
ei_ratio_min (float)
ei_ratio_max (float)
support_asymmetry_rtol (float)
fit_cond_warn (float)
fit_cond_fail (float)
platform_cm_offset_gyradius_factor (float)

ei_ratio_max: float = 10.0

ei_ratio_min: float = 0.1

fit_cond_fail: float = 1000000.0

fit_cond_warn: float = 10000.0

platform_cm_offset_gyradius_factor: float = 1.0

stiffness_jump_factor: float = 5.0

support_asymmetry_rtol: float = 1e-06

class FitOptions(polynomial_rms_threshold=0.09, torsion_contamination_threshold=0.1, fit_cond_warn=10000.0, fit_cond_fail=1000000.0)[source]

Bases: object

Polynomial-fit + family-selection thresholds.

Variables:

polynomial_rms_threshold (float, default 0.09) – Maximum RMS residual a clamped-base polynomial fit is allowed to have for its mode to count as a viable FA / SS family candidate (after the tangent-line root-rigid subtraction). The _select_tower_family algorithm drops candidates above this.
torsion_contamination_threshold (float, default 0.10) – Drops candidates whose modal-kinetic-energy torsion fraction T_tor = sum(phi_tor**2) / sum(phi_total**2) exceeds this. Hybrid bending + twist modes (T_tor in the 1-3 % range are “near-pure-bending”) cannot be expressed by the polynomial ansatz; the filter keeps them out of the FA / SS selection.
fit_cond_warn (float, default 1e4) – Condition-number ceiling above which the polynomial-fit design-matrix conditioning emits a RuntimeWarning. Suggests numerical sensitivity of the polynomial-coefficient solve to perturbations in the input mode shape.
fit_cond_fail (float, default 1e6) – Condition-number ceiling above which the polynomial fit is flagged as FAIL by pybmodes.checks.check_model(). The reconstructed shape may still be a good visual fit, but the coefficient values are not trustworthy beyond a couple of significant figures.

Parameters:

polynomial_rms_threshold (float)
torsion_contamination_threshold (float)
fit_cond_warn (float)
fit_cond_fail (float)

fit_cond_fail: float = 1000000.0

fit_cond_warn: float = 10000.0

polynomial_rms_threshold: float = 0.09

torsion_contamination_threshold: float = 0.1

class SolverOptions(sparse_ndof_threshold=500, symmetry_rtol=1e-12)[source]

Bases: object

FEM solver dispatch + matrix-symmetry tolerances.

Variables:

sparse_ndof_threshold (int, default 500) – Above this many reduced-system DOFs and when a small subset of modes is requested, the solver routes through scipy.sparse.linalg.eigsh() shift-invert instead of the dense scipy.linalg.eigh(). 5-18 x speedup on real-tower meshes; below the threshold the dense path’s LAPACK back-end is faster.
symmetry_rtol (float, default 1e-12) – Relative tolerance (multiplied by max|matrix|) for treating an assembled stiffness or mass matrix as symmetric. Asymmetry beyond this routes the eigenproblem through the general dense scipy.linalg.eig() instead of the symmetric scipy.linalg.eigh(). The OC3 Hywind cross-coupled hydro_K + mooring_K exercises this branch.

Parameters:

sparse_ndof_threshold (int)
symmetry_rtol (float)

sparse_ndof_threshold: int = 500

symmetry_rtol: float = 1e-12

Report generation

Generate a human-readable report on a pyBmodes modal analysis.

The entry point is generate_report(). It builds a structured representation of the report (sections + tables) once, then renders that representation in one of three output formats:

Markdown (format="md") — default. Human-readable plain text with section headers, tables, and inline code blocks. Useful for paste-into-issue or check-in-to-repo workflows.
HTML (format="html") — a self-contained HTML5 document with inline CSS. Same content as the markdown version; emitted directly rather than via an md→html converter so the report module has no runtime dependency on the third-party markdown package.
CSV (format="csv") — frequencies + polynomial-coefficient rows only, suitable for spreadsheet ingestion. The narrative sections (assumptions, validation verdict, warnings) are dropped because CSV can’t represent them faithfully.

The report’s content is driven by what the caller supplies. The minimum is a ModalResult; everything else (model, campbell, validation, check_warnings, tower_params, blade_params) is optional and unlocks the corresponding section when supplied.

generate_report(result, output_path, *, format='md', model=None, campbell=None, validation=None, check_warnings=None, tower_params=None, blade_params=None, elastodyn_compatible=None, source_file=None, status=None)[source]

Render a modal-analysis report to output_path.

source_file is recorded in the Model summary section and the metadata block. When omitted, the function falls back to result.metadata["source_file"] and then to model._bmi.source_file if those are populated.

status is an optional completeness stamp shown at the top of the Model summary section ("complete" / "partial" / "screening" / "invalid"). A workflow that may skip part of its output (e.g. pybmodes.workflows.windio.run_windio()) passes it so a reader can tell at a glance whether the report is the full analysis or a reduced one. None omits the row.

Returns the resolved output path so the caller can chain.

Parameters:

result (ModalResult)
output_path (str | pathlib.Path)
format (ReportFormat)
model (Tower | RotatingBlade | None)
campbell (CampbellResult | None)
validation (ValidationResult | None)
check_warnings (list[ModelWarning] | None)
tower_params (TowerElastoDynParams | None)
blade_params (BladeElastoDynParams | None)
elastodyn_compatible (bool | None)
source_file (str | pathlib.Path | None)
status (str | None)

Return type:

pathlib.Path

Mooring

Quasi-static mooring linearisation for pyBmodes.

Solves the extensible elastic catenary per line, sums fairlead tensions into a platform 6-DOF restoring force, and returns a finite-difference 6×6 mooring stiffness matrix ready for the PlatformSupport.mooring_K block.

References

Jonkman, J. M. (2007). Dynamics Modeling and Loads Analysis of an Offshore Floating Wind Turbine, NREL/TP-500-41958. Appendix B, equations B-1 / B-2 are the extensible-catenary boundary-condition equations implemented in pybmodes.mooring._catenary._catenary_residual(). B-7 / B-8 are the seabed-contact variants (no friction; CB = 0).
Irvine, H. M. (1981). Cable Structures, MIT Press, §2.4 the extensible elastic catenary. Equivalent derivation; the EA-correction terms H · L / EA and L · (V − ½WL) / EA come from §2.4 eqn (2.49).

Module layout

Phase 3 PR C2 of the v1.x architecture refactor split this from a single 1202-line module into a sub-package. The public API is unchanged; internal helpers live in private sub-modules:

pybmodes.mooring.types — LineType, Point, Line dataclasses. Line.solve_static() is here.
pybmodes.mooring.system — MooringSystem: multi-line force assembly, equilibrium Newton, 6×6 stiffness, plus the MooringSystem.from_moordyn() and MooringSystem.from_windio_mooring() classmethod parsers.
_catenary — extensible-elastic-catenary residual + analytical 2×2 Jacobian (Jonkman 2007 B-1 / B-2 + B-7 / B-8).
_rotation — 3-2-1 intrinsic Euler rotation primitive (ElastoDyn attitude convention).
_moordyn_parser — MoorDyn .dat section + row tokenisers.

Scope

Implemented:

Extensible elastic catenary per line (Newton on (H, V); analytical 2 × 2 Jacobian; tol = 1e-6 m, MaxIter = 100).
Fully-suspended (V_F ≥ W · L) and anchor-on-seabed (V_F < W · L, zero friction) profiles, branched inside the residual function.
Multi-line platform restoring force from a 6-DOF body offset.
Central-difference linearisation around an arbitrary or zero offset producing K_mooring (6, 6).
MoorDyn v1 (CONNECTION) and v2 (POINT) .dat parsing.

Known limitations:

Seabed friction (CB > 0) is parsed from LineType but not consumed by the catenary solver.
Sloped seabed, U-shape lines (one line touching the seabed mid- span), and the vertical-line degenerate case (H → 0) are not handled.
Time-domain dynamics, hydrodynamic drag, and added mass on the lines themselves are out of scope — this is a quasi-static linearised model only.
MooringSystem.solve_equilibrium() defaults to the input offset; pure mooring has no z equilibrium without buoyancy / weight, so the Newton iteration is only meaningful for the in-plane DOFs (surge, sway, yaw) of a 3-fold-symmetric layout. Callers wanting platform equilibrium under a full force model should call MooringSystem.restoring_force() and assemble the rest of the forces themselves.

Coordinate / unit conventions

SI throughout: m, N, kg, kg/m, N/m; radians for rotations. Origin at MSL; z positive upward; anchors at negative z (below MSL). Matches OpenFAST / HydroDyn / ElastoDyn without any coordinate transform.

Body rotation uses the 3-2-1 (z-y-x intrinsic) Euler angle convention R = R_z(yaw) · R_y(pitch) · R_x(roll) — the same convention as ElastoDyn’s platform 6-DOF state.

class Line(line_type, point_a, point_b, unstretched_length, seabed_contact=True)[source]

Bases: object

A single mooring line connecting two Point endpoints.

The catenary solve is owned by this class; the multi-line force assembly is delegated to MooringSystem.

Variables:

seabed_contact (bool, default True) – Whether the anchor sits on a seabed. When True, a solver iterate with V_F < W·L triggers the seabed-contact branch of the catenary equations (Jonkman 2007 B-7 / B-8 with CB = 0); the anchor-side portion of the line is treated as resting on the seabed. When False the fully-suspended equations (B-1 / B-2) are used unconditionally — appropriate for analytical tests where both endpoints are in free air and the line just sags between them. FOWT use cases default to True.

Parameters:

line_type (LineType)
point_a (Point)
point_b (Point)
unstretched_length (float)
seabed_contact (bool)

line_type: LineType

point_a: Point

point_b: Point

seabed_contact: bool = True

solve_static(r_a, r_b, tol=1e-06, max_iter=100)[source]

Solve the elastic catenary between r_a (anchor) and r_b (fairlead).

Implements Jonkman 2007 Appendix B equations B-1 and B-2 for the fully-suspended branch; B-7 / B-8 with CB = 0 for the seabed- contact branch. The two unknowns are H (horizontal tension, constant along the line) and V_F (vertical tension at the fairlead, positive when the line pulls the fairlead downward).

Returns:

H (float) – Horizontal tension (N).
V_F (float) – Vertical fairlead tension (N).
f_on_fairlead (ndarray, shape (3,)) – World-frame force the line exerts on r_b (the fairlead): horizontal component pulls toward r_a, vertical component is -V_F (line pulls fairlead down).

Parameters:

r_a (ndarray)
r_b (ndarray)
tol (float)
max_iter (int)

Return type:

tuple[float, float, ndarray]

unstretched_length: float

class LineType(name, diam, mass_per_length_air, EA, w, CB=0.0)[source]

Bases: object

Material spec for a mooring line.

Variables:

name (str) – Identifier referenced by Line.line_type and by MoorDyn LINES rows.
diam (float) – Outer diameter (m).
mass_per_length_air (float) – Mass density in air (kg/m).
EA (float) – Axial stiffness (N). Inextensible limit is EA → ∞.
w (float) – Wet weight per unit length (N/m). For a homogeneous line of diameter d in water of density ρ_w: w = (m - ρ_w · π/4 · d²) · g.
CB (float, default 0) – Seabed friction coefficient (sliding friction along the resting segment of a partly-grounded line). Parsed for round-trip identity with MoorDyn .dat inputs; the current catenary solver only handles the CB = 0 (frictionless seabed) case.

Parameters:

name (str)
diam (float)
mass_per_length_air (float)
EA (float)
w (float)
CB (float)

CB: float = 0.0

EA: float

diam: float

mass_per_length_air: float

name: str

w: float

class MooringSystem(depth, rho=1025.0, g=9.80665, line_types=None, points=None, lines=None)[source]

Bases: object

A collection of catenary lines connecting a platform to anchors.

The system is fully assembled by the constructor or from_moordyn(). The downstream API is:

fairlead_positions() — world-frame fairlead positions at a given body offset.
restoring_force() — 6-vector force / moment from all lines on the body, in world frame, about the body origin.
solve_equilibrium() — Newton iteration over body DOFs to drive restoring_force to zero (may not converge for pure mooring without buoyancy — see module docstring).
stiffness_matrix() — finite-difference 6 × 6 about a chosen offset (or zero by default; see note in the docstring below).

Parameters:

depth (float)
rho (float)
g (float)
line_types (dict[str, LineType] | None)
points (dict[int, Point] | None)
lines (list[Line] | None)

fairlead_positions(body_r6)[source]

World-frame positions of every Vessel-attached point.

Parameters:: body_r6 (ndarray)
Return type:: list[ndarray]

classmethod from_moordyn(dat_path, rho=1025.0, g=9.80665)[source]

Parse a MoorDyn v1 / v2 .dat and return a populated system.

Sections recognised:

LINE TYPES (or LINE DICTIONARY): rows Name Diam MassDenInAir EA …. Wet weight is derived as w = (MassDenInAir − ρ · π/4 · d²) · g.
POINTS or CONNECTION PROPERTIES: rows ID Attachment X Y Z …. Attachment accepted case-insensitively as Fixed, Vessel, or Free.
LINES or LINE PROPERTIES: rows ID LineType AttachA AttachB UnstrLen ….
OPTIONS (or SOLVER OPTIONS): WtrDpth / depth and rhoW / rho if present override the constructor defaults.

Each section header is detected by startswith('---') plus a keyword (case-insensitive); the immediately-following 1-2 rows are skipped as column headers / units rows.

Parameters:

dat_path (Path | str)
rho (float)
g (float)

Return type:

MooringSystem

classmethod from_windio_mooring(floating, *, depth, moordyn_fallback=None, rho=1025.0, g=9.80665)[source]

Build a system from a WindIO mooring block (issue #35).

floating is a parsed pybmodes.io.windio_floating.WindIOFloating — its joints table supplies every anchor / fairlead world position (fairleads are the axial joints resolved during parsing). Line topology (nodes / lines) comes from the WindIO mooring block; line properties (mass/length, EA, wet weight) are resolved in order of preference:

explicit WindIO line_types fields — mass_density / linear_density and stiffness / EA / axial_stiffness — when present;
a companion MoorDyn deck (moordyn_fallback): the accurate path, equivalent to how WISDEM/RAFT delegate chain sizing to MoorPy MoorProps (line types matched by name, or the sole entry);
a documented studless-chain diameter regression (MoorPy MoorProps default coefficients m ≈ 19.9e3·d², EA ≈ 0.854e11·d²) — a rough last resort that emits a UserWarning; supply a deck or explicit props for quantitative work.

Catenary engine + stiffness_matrix are unchanged, so the WindIO-topology system and the companion-MoorDyn system are consistent by construction (cross-path consistency anchor).

Parameters:

depth (float)
moordyn_fallback (Path | str | None)
rho (float)
g (float)

Return type:

MooringSystem

restoring_force(body_r6)[source]

6-vector force/moment from all lines on the platform body.

F[:3] = sum of world-frame forces at every Vessel-attached endpoint; F[3:6] = sum of moments (r_endpoint_world − r_body_origin) × F_endpoint, about the body origin.

For each line, the endpoint forces are derived from the catenary solve in this order:

F_on_B (B = the “fairlead” passed as r_b to solve_static) = H · ê_B→A + (-V_F) ẑ.
F_on_A = H · ê_A→B + V_A ẑ where V_A = max(0, V_F − W·L). Fully suspended: V_A = V_F − W·L > 0 (line pulls anchor up). Seabed contact (V_F < W·L, CB = 0): the anchor is on the seabed; horizontal tension at the anchor is still H (no friction decay), and V_A = 0.

Lines with neither endpoint attached to the body contribute nothing. Lines with both endpoints attached to the body contribute both endpoint reactions.

Parameters:: body_r6 (ndarray)
Return type:: ndarray

solve_equilibrium(body_r6_init=None, tol=0.0001, max_iter=50, dx=0.1, dtheta=0.1)[source]

Newton iteration over body 6-DOF to drive restoring_force to zero.

Warning: pure mooring without buoyancy / weight has no z equilibrium (the lines always pull the platform down). For a 3-fold-symmetric mooring at zero offset the in-plane DOFs (surge, sway, yaw) are already balanced; the heave DOF will not converge. Pass a body_r6_init close to your expected operating point and accept the result as a “best effort” local-minimum.

Parameters:

body_r6_init (ndarray | None)
tol (float)
max_iter (int)
dx (float)
dtheta (float)

Return type:

stiffness_matrix(body_r6=None, dx=0.1, dtheta=0.1)[source]

Linearised 6 × 6 mooring stiffness about body_r6.

Central finite differences with perturbation dx (m) on translational DOFs and dtheta (rad) on rotational DOFs. The trans-rot off-diagonal blocks are symmetrised after differencing — mooring linearised at static equilibrium is the Hessian of a potential and must therefore be symmetric; finite-difference noise gets averaged out.

body_r6 = None is treated as np.zeros(6) (the typical FOWT linearisation point). Pure mooring has no z-direction equilibrium without buoyancy, so a solve-for-equilibrium default would diverge; pass an explicit body_r6 if you want a different linearisation point.

Returns:

K – Stiffness in N/m / N / N·m/rad block structure (trans-trans: N/m, rot-trans / trans-rot: N (= N·m/m), rot-rot: N·m/rad).

Return type:

ndarray, shape (6, 6)

Parameters:

body_r6 (ndarray | None)
dx (float)
dtheta (float)

class Point(id, attachment, r_body)[source]

Bases: object

Endpoint of a mooring line.

Variables:

id (int) – MoorDyn point ID (preserved for round-trip identification).
attachment (str) – One of Fixed / Vessel / Free (case-insensitive on construction; stored title-cased).
r_body (ndarray, shape (3,)) – Position in body frame for Vessel points; world frame for Fixed and Free points. Free points are essentially Fixed placeholders for this quasi-static solver — they don’t participate in the body-equilibrium DOFs.

Parameters:

id (int)
attachment (str)
r_body (ndarray)

attachment: str

id: int

r_body: ndarray

r_world(body_r6)[source]

World-frame position for this point at platform state body_r6.

For Fixed and Free points the platform state is ignored; for Vessel points the rotation R(roll, pitch, yaw) · r_body + r_body_origin is applied with the 3-2-1 intrinsic Euler convention.

Parameters:: body_r6 (ndarray)
Return type:: ndarray

Soil-pile interaction

Mudline coupled-spring foundation for soft monopile soil-pile interaction.

The coupled-spring (CS) model represents the soil-pile reaction at the mudline as three springs: a lateral spring K_hh, a rotational spring K_rr, and a cross-coupling K_hr. The model is well established for monopile-supported offshore wind turbines and is endorsed by Yu and Amdahl (2023) for first three modes when the spring stiffnesses are calculated properly.

This module wires closed-form formulas for K_hh, K_hr, K_rr from pile geometry and soil properties into a 6x6 mooring_K block that drops straight into pybmodes.io.bmi.PlatformSupport of a hub_conn = 3 soft-monopile BMI. The dispatch covers two formula families and three soil profiles, classifying pile behaviour via Randolph (1981).

Scope. MudlineFoundation produces the linearised mudline stiffness used for coupled-frequency prediction. For ElastoDyn polynomial coefficient generation the cantilever path (pybmodes.models.Tower.from_elastodyn() or pybmodes.models.Tower.from_geometry()) is still required regardless of soil flexibility. The mudline stiffness affects the coupled-system frequency but NOT the polynomial basis. ElastoDyn’s SHP ansatz requires clamped-base mode shapes (src/pybmodes/_examples/reference_decks/FLOATING_CASES.md records the source-code citations and cases/ECOSYSTEM_FINDING.md the audit trail).

References

Yu, Z. and Amdahl, J. (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. https://doi.org/10.1016/j.marstruc.2023.103482
Shadlou, M. and Bhattacharya, S. (2016). Dynamic stiffness of monopiles supporting offshore wind turbine generators. Soil Dynamics and Earthquake Engineering 88, 15-32.
Psaroudakis, E. G., Mylonakis, G. and Antonopoulos, A. (2021). Analytical formulas for the lateral response of monopiles in homogeneous Winkler-type foundations. (As cited in Yu and Amdahl 2023, Eq 25.)
Randolph, M. F. (1981). The response of flexible piles to lateral loading. Geotechnique 31, 247-259.

class MudlineFoundation(K_hh, K_hr, K_rr, pile_behaviour, soil_profile, formula)[source]

Bases: object

Three coupled springs at the mudline for a monopile foundation.

Spring convention matches Eq (3) of Yu and Amdahl (2023): the mudline force-moment vector is [F, M] = [[K_hh, K_hr], [K_hr, K_rr]] @ [rho, theta] where rho is the mudline lateral displacement and theta the mudline rotation in the same 2-D plane. By the right-hand convention used in OpenFAST (Jonkman 2010 NREL/TP-500-47535 Table 5-1), the K_hr term is negative for typical sands and clays.

The as_mooring_K() accessor maps the 2-D coupled-spring matrix to the OpenFAST 6-DOF order [surge, sway, heave, roll, pitch, yaw] and is symmetric. Heave and yaw are not modelled by this CS surrogate and stay at zero in the returned matrix. Wire it into the soft monopile via pybmodes.io.bmi.PlatformSupport.mooring_K of a BMI built for hub_conn = 3.

Parameters:

K_hh (float)
K_hr (float)
K_rr (float)
pile_behaviour (Literal['flexible', 'rigid'])
soil_profile (Literal['homogeneous', 'parabolic', 'linear'])
formula (Literal['shadlou', 'psaroudakis'])

K_hh: float

K_hr: float

K_rr: float

as_mooring_K()[source]

Return the 6x6 mudline stiffness in OpenFAST DOF order.

Mapping (Eq 3 of Yu and Amdahl 2023 lifted into 6 DOFs):

K[0, 0] = K[1, 1] = K_hh (lateral on surge and sway, the three-fold axisymmetric assumption that pins both diagonals to the same value).
K[3, 3] = K[4, 4] = K_rr (rotational on roll and pitch).
K[0, 4] = K[4, 0] = K_hr (surge-pitch coupling).
K[1, 3] = K[3, 1] = -K_hr (sway-roll coupling, opposite sign by right-hand rule).
K[2, 2] = 0 (heave not modelled by the CS surrogate).
K[5, 5] = 0 (yaw not modelled by the CS surrogate).

The cross-coupling sign convention is pinned by tests/test_mooring.py::test_oc3hywind_bmi_dof_order_matches_jonkman_2010 against Jonkman (2010) OC3 Table 5-1, which reports K_15 = -2.821e6 N for surge-pitch and K_24 = +2.816e6 N for sway-roll.

Return type:: ndarray

formula: Literal['shadlou', 'psaroudakis']

classmethod from_soil_properties(pile_diameter, pile_length_embedded, pile_EI, soil_E, soil_nu=0.3, soil_profile='homogeneous', pile_behaviour='auto', formula='shadlou')[source]

Compute K_hh, K_hr, K_rr from pile geometry and soil properties.

Parameters:

pile_diameter (float) – Outer diameter of the monopile, D_P in m.
pile_length_embedded (float) – Embedded pile length, L_P in m.
pile_EI (float) – Pile bending stiffness E_P * I_P in N m^2. For a tubular pile with diameter D and wall thickness t, I_P = pi / 64 * (D^4 - (D - 2 t)^4).
soil_E (float) – Soil Young’s modulus E_SO in Pa. For an inhomogeneous profile this is the reference modulus at the depth used by the chosen formula family.
soil_nu (float) – Soil Poisson’s ratio. Default 0.3.
soil_profile (Literal['homogeneous', 'parabolic', 'linear']) – One of "homogeneous" (constant E_SO with depth), "parabolic" (E_SO proportional to sqrt of depth), "linear" (E_SO proportional to depth).
pile_behaviour (str) – One of "flexible", "rigid", "auto". When set to "auto", the pile is classified per Randolph (1981). An intermediate L/D ratio falls back to the flexible formulas with a UserWarning.
formula (Literal['shadlou', 'psaroudakis']) – "shadlou" (default) uses Shadlou and Bhattacharya (2016) per Yu Table 1 and covers all three soil profiles. "psaroudakis" uses Yu Eq 25 and is restricted to the homogeneous profile.

Returns:

Coupled-spring stiffness in SI units.

Return type:

MudlineFoundation

Raises:

ValueError – On non-positive geometry or soil parameters, an unknown soil_profile or pile_behaviour token, or formula="psaroudakis" paired with a non-homogeneous soil profile.

property pile_behavior: Literal['flexible', 'rigid']: US-spelling alias preserved for prompt-style external code.

pile_behaviour: Literal['flexible', 'rigid']

soil_profile: Literal['homogeneous', 'parabolic', 'linear']

Plot helpers (`[plots]` extra)

Professional plotting utilities for pybmodes results.

Requires matplotlib >= 3.7. Install with:

pip install "pybmodes[plots]"

All functions return a matplotlib.figure.Figure object; call fig.show() or fig.savefig(path) as needed.

apply_style()[source]

Apply the MATLAB-style plot defaults.

Mutates matplotlib.rcParams in place. Safe to call multiple times — the second call simply re-applies the same values. Raises ImportError if matplotlib isn’t installed (the optional [plots] extra).

Return type:: None

bir_mode_shape_plot(result, mode_specs, *, title=None, height_label='Tower section height / H', x_label='Modal displacement', annotations=None, coupling_overlay=None, figsize=(5.5, 6.5), xlim=None, colors=None)[source]

Plot mode shapes in the Bir 2010 figure convention.

Parameters:

result (ModalResult) – ModalResult from Tower.run() or RotatingBlade.run().
mode_specs (list[ModeSpec]) – List of (mode_index_1based, component, label) tuples. component is one of "flap" (fore-aft / F-A), "lag" (side-side / S-S), "twist", or "axial". label appears in the legend.
title (str | None) – Optional figure title.
height_label (str) – Y-axis label. Default matches Bir’s notation; pass "Span fraction" for blade plots.
x_label (str) – X-axis label. Default "Modal displacement" matches the paper.
annotations (dict[str, float] | None) – Optional {label: y_fraction} dict drawing horizontal markers at the given normalised heights (e.g. {"Mean Sea Level": 0.40, "Mud Line": 0.25} for Bir Fig 8).
coupling_overlay (list[ModeSpec] | None) – Optional list of (mode_index, component, label) plotted as dashed lines on the same axes; used to show e.g. the twist component of an S-S mode (Bir Fig 5b, 6b).
figsize (tuple[float, float]) – Matplotlib figure size in inches.
xlim (tuple[float, float] | None) – Optional (xmin, xmax); auto-fits with a small pad if omitted.
colors (list | str | None) – Optional colour control (issue #47): a matplotlib colormap name or an explicit list of colours. None (default) uses the styled palette, switching to a continuous colormap when there are more curves than distinct palette colours so no two modes share a hue.

Return type:

bir_mode_shape_subplot(result, panels, *, suptitle=None, height_label='Tower section height / H', x_label='Modal displacement', annotations=None, figsize=None, colors=None)[source]

Multi-panel Bir-convention plot (matches Bir Fig 8 layout).

Parameters:

panels (list[tuple[str, list[ModeSpec]]]) – List of (panel_title, mode_specs) tuples; one subplot per entry.
annotations (dict[str, float] | None) – Drawn on every panel (e.g. MSL + Mud Line).
colors (list | str | None) – Optional colour control (issue #47): a matplotlib colormap name or an explicit list of colours, applied per panel. None (default) uses the styled palette, switching to a continuous colormap when a panel has more curves than distinct palette colours.
result (ModalResult)
suptitle (str | None)
height_label (str)
x_label (str)
figsize (tuple[float, float] | None)

Return type:

blade_fit_pairs(result, params)[source]

Build (label, span_loc, fem_disp, fit) entries for blade modes.

Returns entries for the 1st flap, 2nd flap, and 1st edge modes, matching the order in params.

Parameters:

result (ModalResult) – ModalResult from RotatingBlade.run().
params (BladeElastoDynParams) – BladeElastoDynParams from compute_blade_params(result).

Return type:

list[FitEntry]

jonswap_spectrum(f, *, hs, tp, gamma=3.3)[source]

JONSWAP wave elevation spectrum S(f) (frequency in Hz).

Standard Pierson-Moskowitz core times the peak-enhancement factor gamma. The shape is scaled so the zeroth spectral moment is exactly the significant-wave-height identity m0 = Hs**2 / 16 (Hs = 4 sqrt(m0)). The peak sits at f_p = 1/Tp. Returns 0 for f <= 0.

Parameters:

f (ndarray)
hs (float)
tp (float)
gamma (float)

Return type:

kaimal_spectrum(f, *, mean_speed, length_scale, sigma=None, turbulence_intensity=0.14)[source]

One-sided longitudinal Kaimal turbulence spectrum S_u(f).

S_u(f) = 4 sigma_u^2 (L/U) / (1 + 6 f L/U)^(5/3) (IEC 61400-1 form), monotonically decreasing in f with the finite low-frequency plateau S_u(0) = 4 sigma_u^2 L / U. Units are m^2/s (PSD of wind speed) when mean_speed is in m/s and length_scale in m; the plot normalises it, so only the shape matters there.

sigma (the longitudinal standard deviation) defaults to turbulence_intensity * mean_speed when not given.

Parameters:

f (ndarray)
mean_speed (float)
length_scale (float)
sigma (float | None)
turbulence_intensity (float)

Return type:

plot_environmental_spectra(*, tower_fa_hz=None, tower_ss_hz=None, rpm_design=None, rpm_constraint=None, harmonics=(1, 3), wind=None, wave=None, freq_max=0.6, n_points=2000, ax=None, title=None)[source]

Draw the environmental-loading frequency-placement diagram.

Parameters:

tower_fa_hz (float | None) – Tower 1st fore-aft / side-side natural frequencies (Hz), drawn as a solid / dashed black vertical line. Either may be None to omit.
tower_ss_hz (float | None) – Tower 1st fore-aft / side-side natural frequencies (Hz), drawn as a solid / dashed black vertical line. Either may be None to omit.
rpm_design (tuple[float, float] | None) – (rpm_min, rpm_max) of the actual rotor operating range — the darker hatched design band for each requested harmonic.
rpm_constraint (tuple[float, float] | bool | None) – (rpm_min, rpm_max) of the allowable placement window — the lighter solid constraint band drawn behind the design band. Defaults to the design range widened by +/-15 % when omitted (None). Pass False to suppress the constraint band (and its legend entries) entirely and draw only the operating design band — for callers who have a fixed operating range but no separate placement envelope.
harmonics (Sequence[int]) – Per-rev orders to draw (default (1, 3) -> 1P and 3P).
wind (dict | None) – dict of kaimal_spectrum() keyword arguments (mean_speed, length_scale, optionally sigma / turbulence_intensity). None skips the wind curve.
wave (dict | None) – dict of jonswap_spectrum() keyword arguments (hs, tp, optionally gamma). None skips the wave curve.
freq_max (float) – Frequency-axis upper bound (Hz) and sample count.
n_points (int) – Frequency-axis upper bound (Hz) and sample count.
ax (Axes | None) – Existing matplotlib Axes to draw into; a new figure is created when None.
title (str | None) – Optional figure title.

Return type:

plot_fit_quality(fits, title=None, figsize=None)[source]

Plot polynomial fit vs FEM data with residuals for each mode.

Each subplot shows: * FEM nodal values (circles) * Fitted polynomial (solid line over a fine grid) * Residual band (shaded region between FEM and fit) * RMS residual and tip-slope annotation

Parameters:

fits (list[FitEntry]) – List of (label, span_loc, fem_disp, fit) tuples as returned by blade_fit_pairs() or tower_fit_pairs().
title (str | None) – Overall figure title.
figsize (tuple[float, float] | None) – Matplotlib figure size. Defaults to (4.5 * n_cols, 4.0 * n_rows).

Return type:

plot_mode_shapes(result, n_modes=6, component='both', title=None, figsize=None, *, normalize='mode', colors=None)[source]

Plot normalised mode shape displacements vs normalised span.

Parameters:

result (ModalResult) – Output from RotatingBlade.run() or Tower.run().
n_modes (int) – Number of modes to overlay (at most len(result.shapes)).
component (str) – "flap" — fore-aft (w) panel only. "lag" — side-side (v) panel only. "both" — two side-by-side panels (default).
title (str | None) – Overall figure title. Uses a sensible default when None.
figsize (tuple[float, float] | None) – Matplotlib figure size (width_in, height_in). Defaults to (9, 4) for one panel and (14, 4) for two panels.
normalize (str) –
How each mode’s two curves are scaled (issue #47):
- "mode" (default) — both panels share one scale per mode, the peak |displacement| across flap and lag. The dominant direction then reaches ±1 and the minor direction stays proportionally small, so a fore-aft mode reads as a full-amplitude curve in the flap panel and a flat near-zero line in the lag panel. This is what makes FA vs SS (and a rigid platform DOF vs a real bending mode) visually distinguishable.
- "component" — the legacy behaviour: each panel is normalised independently to its own peak, so even a 1 % cross-coupling component is blown up to full height and a predominantly-flap mode looks identical in both panels. Kept for reproducing pre-1.5 figures.
colors (list | str | None) – Optional colour control passed to the internal palette picker: a matplotlib colormap name (sampled into n_modes distinct colours) or an explicit list of colours. None (default) uses the styled prop_cycle, automatically switching to a continuous colormap when there are more modes than distinct palette colours so no two modes collide on the same hue.

Return type: