API reference#

Model#

class femorph_solver.model.Model(grid: UnstructuredGrid | None = None, *, unit_system: UnitSystem | str = UnitSystem.UNSPECIFIED)[source]#

Bases: object

apply_force(node: int, *, fx: float = 0.0, fy: float = 0.0, fz: float = 0.0, mx: float = 0.0, my: float = 0.0, mz: float = 0.0) → None[source]#

Apply a keyword-driven nodal force / moment vector.

Sparse by construction — zero-valued components are skipped. Any of the six ForceLabel entries can be supplied as keyword arguments; unused ones stay zero (and unregistered on the Model).

Parameters:

node – MAPDL node number the load is applied at.
fx – Translational force components (units inherited from the Model’s unit_system).
fy – Translational force components (units inherited from the Model’s unit_system).
fz – Translational force components (units inherited from the Model’s unit_system).
mx – Moment components. Relevant for elements with rotational DOFs (BEAM2 / QUAD4_SHELL / etc.).
my – Moment components. Relevant for elements with rotational DOFs (BEAM2 / QUAD4_SHELL / etc.).
mz – Moment components. Relevant for elements with rotational DOFs (BEAM2 / QUAD4_SHELL / etc.).

Examples

>>> model.apply_force(tip_node, fz=-100.0)      # 100 N down
>>> model.apply_force(pin_node, fx=200.0, mz=5.0)

assert_unit_system(expected: UnitSystem | str) → None[source]#

Raise ValueError if this Model’s unit_system is not expected.

Use this in scripts that require a specific unit convention — a mismatch between the caller’s assumption and the Model’s stamp usually means numbers will be silently wrong by a factor of \(10^n\). UnitSystem.UNSPECIFIED always passes (no assertion possible on an unstamped Model).

Declare an element kernel and a material in one call.

A single call replaces the common three-step MAPDL sequence et(id, name) / mp(prop, mat_id, val) × n / r(real_id, ...).

Parameters:

element – Neutral name ("HEX8", "HEX20", "TET10", …) or MAPDL alias ("SOLID185" etc.) — either works. Register the kernel under et_id (default 1).
material – Mapping of material properties ({"EX": ..., "PRXY": ..., "DENS": ...}). Values are stamped under mat_id (default 1). Accepts the output of femorph_solver.materials.fetch() directly.
et_id – Numeric ids under which to register the kernel / material / real constants. Only relevant when a model mixes multiple element types or materials; most single-physics demos use the defaults.
mat_id – Numeric ids under which to register the kernel / material / real constants. Only relevant when a model mixes multiple element types or materials; most single-physics demos use the defaults.
real_id – Numeric ids under which to register the kernel / material / real constants. Only relevant when a model mixes multiple element types or materials; most single-physics demos use the defaults.
real – Real-constant values (array-like of floats) to stamp under real_id. None skips the real-constant declaration.

d(node: int, label: str | DOFLabel, value: float = 0.0) → None[source]#

MAPDL D — constrain a nodal DOF.

label accepts a raw string, a DOFLabel member, or the sentinel string "ALL", which expands to every DOF active at that node (as determined by the elements that touch it).

dof_map() → ndarray[source]#

Return an (N, 2) array of (mapdl_node_num, local_dof_idx).

Row i of K/M corresponds to dof_map()[i]. local_dof_idx is the 0-based position in MAPDL’s canonical order (UX, UY, UZ, ROTX, ROTY, ROTZ). Nodes only get rows for the DOFs the elements at them actually need, so a truss model stays 3-DOF while a beam model gets 6.

e(*node_nums: int, num: int | None = None) → int[source]#

MAPDL E — define an element from node numbers.

Uses the current TYPE / MAT / REAL stamps. Returns the element number (auto-assigned to max(existing) + 1 when num is omitted).

eel(u: ndarray, *, nodal_avg: bool = True) → ndarray | dict[int, ndarray][source]#

Elastic strain from a displacement vector (MAPDL S,EPEL equivalent).

Strain is evaluated at each element’s own nodes by ε(ξ_node) = B(ξ_node)·u_e — direct nodal evaluation rather than Gauss extrapolation, but the two agree exactly for uniform-strain fields and converge on mesh refinement. Output is 6-component Voigt with engineering shears [εxx, εyy, εzz, γxy, γyz, γxz].

Parameters:

u ((N,) float64) – Global displacement vector indexed by dof_map() — i.e. the StaticResult.displacement from solve() or a mode shape from modal_solve().
nodal_avg (bool, default True) – If True (MAPDL PLNSOL default), values at shared nodes are averaged across every element touching them — returns a (n_nodes_global, 6) array indexed by node_numbers. If False, returns the per-element dict {elem_num: (n_nodes_in_elem, 6)} — unaveraged, like PLESOL.

element_info(elem_num: int) → tuple[list[int], int, int, int][source]#

Return (node_nums, et_id, mat_id, real_id) for an element.

Convenience accessor for the stress-recovery / strain-export paths that previously reached into model._elements[en]. O(log n) lookup; lift grid.cell_data / grid.cell_connectivity yourself for bulk iteration.

estimate_solve(n_modes: int = 10, *, linear_solver: str = 'auto', eigen_solver: str = 'arpack', ooc: bool = False, host_spec=None)[source]#

Predicted (wall_s, peak_rss_mb) for a modal solve on this model and machine, without running it.

Wraps Estimator — reads training rows from any femorph-benchmark-*.json files in the current working directory and fits per-(host, backend) log-log power laws on the fly. No training data → falls back to a shape-of-universe prior tuned from the repo’s own 14900KF benchmarks; the returned p95 band then sits at 2× p50 so callers can tell they’re looking at a prior, not a calibrated prediction.

Parameters:

n_modes (int) – Modes requested (matches modal_solve()).
linear_solver (str) – Backend names — picks which per-bucket fit applies.
eigen_solver (str) – Backend names — picks which per-bucket fit applies.
ooc (bool, default False) – True applies the measured OOC multipliers (+3.5× wall, −22 % peak) on top of the in-core prediction; useful when modal_solve would otherwise exceed RAM.
host_spec (HostSpec or None) – Override host detection (useful for “what would this take on a different box?” queries). None → auto- detect.

Returns:

p50 / p95 bands + the training-bucket diagnostic the Estimator picked.

Return type:

femorph_solver.estimators.Estimate

et(et_id: int, name: str | ElementType) → None[source]#

MAPDL ET — declare an element type.

name accepts either the neutral topology-first name ("HEX8", "HEX20", "TET10", "WEDGE15", "PYR13", "BEAM2", "QUAD4_SHELL", "QUAD4_PLANE", "TRUSS2", "SPRING", "POINT_MASS"), the MAPDL catalogue name ("SOLID185" / "SOLID186" / "BEAM188" / …), or an ElementType member (ElementType.HEX8). The two forms resolve to the same kernel; neutral is canonical and MAPDL is a compatibility alias.

f(node: int, label: str | ForceLabel, value: float) → None[source]#

MAPDL F — apply a nodal force or moment.

label accepts a raw string ("FZ") or a ForceLabel member.

Pin nodal DOFs to zero (Dirichlet boundary condition).

Parameters:

nodes – Single MAPDL node number, or an array / iterable of node numbers. Mutually exclusive with where.
where – Boolean mask of length Model.n_nodes selecting grid points to pin. Indexed by VTK point order — the same order used by grid.points. Mutually exclusive with nodes.
dof (str or DOFLabel, default "ALL") – DOF to pin. "ALL" clamps every DOF active at the selected node(s); any other value passes through to d() verbatim.

Examples

Clamp everything at the x = 0 face of a bar:

model.fix(where=grid.points[:, 0] < 1e-9)

Pin a single node’s UZ:

model.fix(nodes=42, dof="UZ")

classmethod from_grid(grid: UnstructuredGrid) → Model[source]#

Wrap an existing pyvista.UnstructuredGrid as a Model.

The grid is expected to carry the MAPDL-style cell/point data arrays (ansys_node_num, ansys_elem_num, ansys_elem_type_num, ansys_material_type, ansys_real_constant). Missing arrays are auto-filled with sequential 1-based ids so meshes built directly in pyvista (e.g. a StructuredGrid cast to unstructured) are usable without the caller stamping every field themselves.

The grid becomes the Model’s sole source of truth for geometry and connectivity — the staging dicts stay empty for the lifetime of the Model unless the caller later mutates via n() / e(). Two incidental passes happen here:

cell_connectivity / offset are coerced to int32 if they came in as int64 (VTK’s default on 64-bit builds) — halves the cell-structure footprint and removes a copy-cast the hot path would otherwise do on every assembly call.
Points and cells are reordered so ansys_node_num and ansys_elem_num are strictly increasing. Every downstream hot path assumes MAPDL-sorted order to produce a deterministic DOF layout; canonicalising here keeps the invariant local to Model construction rather than leaking into every consumer.

mass_matrix(lumped: bool = False) → csr_array[source]#

Assemble the global mass matrix as a scipy.sparse.csr_array.

Same caching semantics as stiffness_matrix() — consistent and lumped variants are cached independently.

mat(mat_id: int) → None[source]#: MAPDL MAT — set the current material stamp for E.

modal_solve(n_modes: int = 10, lumped: bool = False, *, eigen_solver: str = 'arpack', linear_solver: str = 'auto', tol: float = 1e-12, release_inputs: bool = True, mixed_precision: bool | None = False, progress: bool = False) → ModalResult[source]#

Run a linear modal analysis (MAPDL ANTYPE,MODAL equivalent).

Solves K φ = ω² M φ with D-staged Dirichlet BCs applied. Returns the lowest n_modes modes, ordered by frequency.

eigen_solver picks the sparse eigenbackend ("arpack", "lobpcg", "primme", "dense"); linear_solver picks the inner shift-invert factorizer. The default "auto" picks the fastest SPD direct solver available (pardiso → cholmod → umfpack → superlu); pass an explicit name ("superlu", "cholmod", …) to override. See femorph_solver.solvers.eigen.list_eigen_solvers() and femorph_solver.solvers.linear.list_linear_solvers().

release_inputs=True (the default) drops the Model’s cached full K / M the moment the BC-reduced K_ff / M_ff are built, and before the Pardiso factor allocates. On a 1.2 M-DOF SOLID186 problem that saves ~2 GB of peak RSS (the full K and M each run ~1 GB + scaling with density).

Time cost: negligible on the first call — CPython’s refcount frees the full K / M immediately when the BC-reduce helper returns (no gc.collect involved; the sparse-matrix graph has no cycles so a full GC would be pure overhead). On repeat calls on the same Model however the full K / M must be re-assembled from elements (~1-3 s at large sizes) because the cache was released; parametric sweeps should pass release_inputs=False to keep the warm-assembly cache.

mixed_precision controls the inner Pardiso factor precision (ignored by non-Pardiso backends):

False (default) — factor and solve in double precision. Identical frequencies to MAPDL at machine-epsilon.
True — factor in single precision and refine with double-precision residuals. Halves the factor’s permanent footprint when the linked MKL honours the request; the refinement typically recovers ≥13-digit eigenvalue accuracy on well-conditioned SPD stiffness matrices.
None — let femorph-solver decide: opts in for Pardiso SPD factorisation above ~50 k free DOFs, where the memory saving matters most. Not recommended for production models that are known to be ill-conditioned (e.g. highly anisotropic composites, degenerate meshes with near-zero-Jacobian cells).

Note

Mixed-precision Pardiso depends on MKL exposing iparm(28) through pypardiso; some MKL builds silently no-op the request. Check lu._solver.get_iparm(28) after a factorize if you need to confirm the request was honoured.

mp(prop: str | MaterialProperty, mat_id: int, value: float) → None[source]#

MAPDL MP — set a scalar material property.

prop accepts either a raw string ("EX") or a MaterialProperty member (MaterialProperty.EX). Unknown names are rejected to catch typos like "EXX" before they silently mask a zero default in an element kernel.

n(num: int, x: float = 0.0, y: float = 0.0, z: float = 0.0) → None[source]#: MAPDL N — define or overwrite a node.

node_coord(node_num: int) → tuple[float, float, float][source]#

Return (x, y, z) for the given MAPDL node number.

The grid is the only permanent storage, so single-node lookups aren’t O(1) the way the old _nodes[nn] dict lookup was — but callers that do this in a hot loop should lift the full grid.points array via Model.grid() instead.

r(real_id: int, *values: float) → None[source]#: MAPDL R — set the real-constant set.

real(real_id: int) → None[source]#: MAPDL REAL — set the current real-constant stamp for E.

release_matrix_cache() → None[source]#

Drop the cached assembled K / M matrices.

Each cached matrix is roughly nnz × 12 bytes — on a 1.2 M-DOF plate with ~100 M nnz that’s ~1.2 GB each, so releasing both saves ~2.4 GB. Call after a modal_solve / static_solve finishes and the caller is done with further analyses on this Model. Subsequent stiffness_matrix() or mass_matrix() calls will re-assemble.

set_unit_system(unit_system: UnitSystem | str) → None[source]#

Declare the UnitSystem this Model’s numeric inputs follow.

Bookkeeping only — the solver’s math is dimensionless. The stamp travels with the Model through assembly, solve, and the result serialiser so downstream readers can interpret the numbers.

solve(*, linear_solver: str = 'auto', progress: bool = False) → StaticResult[source]#

Run a linear static analysis (MAPDL ANTYPE,STATIC equivalent).

Uses the K assembled from current elements, the D-staged Dirichlet BCs, and the F-staged nodal forces. Returns a StaticResult whose displacement/reaction arrays are indexed by dof_map(). linear_solver picks the sparse backend (see femorph_solver.solvers.linear.list_linear_solvers()).

Parameters:

linear_solver (str, default "auto") – Sparse linear backend name.
progress (bool, default False) – Show a per-stage progress display (rich → tqdm → logger fallback). Off by default for scripted use; the CLI front-end flips it on.

stiffness_matrix() → csr_array[source]#

Assemble the global stiffness matrix as a scipy.sparse.csr_array.

Row/column i corresponds to dof_map() row i.

Cached on the Model: repeat calls on an unchanged Model skip the ~2 s assembly phase entirely. The cache is invalidated by any structural (n / e / et) or material (mp / r) mutation.

strain(displacement: ndarray) → dict[int, ndarray][source]#

Compute per-element elastic strain from a global displacement vector.

Returns a dict mapping MAPDL element number → (n_nodes, 6) engineering-shear Voigt strain sampled at each element node. Shares the kernel path with write_static_result() / write_modal_result() but stays in memory — nothing hits disk. Use this to recover strain from any displacement field (static result, one modal shape, a morphed guess) without round-tripping through .pv.

displacement must be indexed by dof_map() (the same layout the solvers produce). Elements whose kernel returns None from eel_batch are skipped.

transient_solve(*, dt: float, n_steps: int, load: Callable[[float], np.ndarray] | np.ndarray | None = None, lumped: bool = False, damping: _sp.csr_array | None = None, u0: np.ndarray | None = None, v0: np.ndarray | None = None, beta: float = 0.25, gamma: float = 0.5) → TransientResult[source]#

Run a linear transient analysis (MAPDL ANTYPE,TRANS equivalent).

Integrates M ü + C u̇ + K u = F(t) with Newmark-β. The load F(t) defaults to the F-staged nodal forces held constant over the interval; pass a callable t -> (N,) for time-varying loads or a fixed array to override the staged values.

type(et_id: int) → None[source]#: MAPDL TYPE — set the current element-type stamp for E.

property unit_system: UnitSystem#

Current UnitSystem stamp.

Defaults to UnitSystem.UNSPECIFIED for legacy Models that pre-date this attribute. Set via set_unit_system(), pass via the constructor, or let Model.from_grid() / femorph_solver.mapdl_api.from_cdb() detect it from the input.

Solvers#

Linear static solve.

Partitions the global system K u = F into free / fixed DOF blocks:

[ K_ff  K_fc ] [ u_f ]   [ F_f ]
[ K_cf  K_cc ] [ u_c ] = [ R_c ]

with u_c prescribed. u_f is recovered with a pluggable linear backend from femorph_solver.solvers.linear. Reaction forces at constrained DOFs are recovered as R_c = K_cf · u_f + K_cc · u_c - F_c (sign convention: R_c is the external reaction needed to hold u_c).

References

Zienkiewicz, O. C., Taylor, R. L., & Zhu, J. Z. The Finite Element Method: Its Basis and Fundamentals, 7th ed., §1 and §9. [the DOF partitioning / master-slave reduction formulated above]
Cook, R. D., Malkus, D. S., Plesha, M. E., & Witt, R. J. Concepts and Applications of Finite Element Analysis, 4th ed., §2. [recovery of reaction forces from the constrained rows]
Bathe, K.-J. Finite Element Procedures, 2nd ed., §3.4. [direct elimination of prescribed DOFs — the form used here, as opposed to penalty or Lagrange-multiplier alternatives]

class femorph_solver.solvers.static.StaticResult(displacement: ndarray, reaction: ndarray, free_mask: ndarray)[source]#

Bases: object

Result of a linear static solve.

displacement#

(N,) DOF-indexed displacement produced by the solve.

Type:: numpy.ndarray

reaction#

(N,) reaction force; nonzero only at constrained DOFs.

Type:: numpy.ndarray

free_mask#

(N,) bool — True at solved-for DOFs.

Type:: numpy.ndarray

save(path: str | Path, model: Model, *, unit_system: str = 'SI', extra_field_data: dict[str, np.ndarray] | None = None) → Path[source]#

Serialise this static result to a disk-backed .pv file.

Re-projects the DOF-indexed displacement (and reaction as the force field) onto per-node (n_points, n_dof_per_node) arrays using the model’s dof_map(), then hands off to write_static_result().

Parameters:

path (str or pathlib.Path) – Destination file (.pv extension conventional).
model (femorph_solver.model.Model) – Model whose K matrix this result was computed from; supplies the mesh + DOF layout.
extra_field_data (mapping, optional) – Extra field_data entries (e.g. solver statistics).

Returns:

Canonicalised absolute path to the written file.

Return type:

pathlib.Path

femorph_solver.solvers.static.solve_static(K: csr_array, F: ndarray, prescribed: dict[int, float], *, linear_solver: str = 'auto', thread_limit: int | None = None) → StaticResult[source]#

Solve K u = F with Dirichlet BCs at the indices in prescribed.

Parameters:

K ((N, N) scipy.sparse.csr_array)
F ((N,) float64 load vector)
prescribed (mapping {global_dof_index: prescribed_value})
linear_solver (str) – Name of a registered linear backend — see femorph_solver.solvers.linear.list_linear_solvers().
thread_limit (int or None, optional) – Cap on BLAS / OpenMP threads used inside the linear solve. None defers to the global limit (see femorph_solver.set_thread_limit()).

Solves K φ = ω² M φ for the lowest n_modes free-vibration modes with Dirichlet BCs partitioned out. DOFs with zero diagonal stiffness (for example the transverse DOFs of a pure axial truss) are dropped as rigid-body modes — MAPDL’s default behaviour.

For n_free <= DENSE_CUTOFF we always use the dense LAPACK path (scipy.linalg.eigh()). For larger systems the eigen backend is configurable via the eigen_solver / linear_solver kwargs — see femorph_solver.solvers.eigen and femorph_solver.solvers.linear for the registered options.

References

Bathe, K.-J. Finite Element Procedures, 2nd ed., Prentice-Hall (2014), §11. [structural eigenproblem K φ = ω² M φ and discussion of subspace / Lanczos methods]
Hughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, §9. [free-vibration modal analysis, M-orthonormalisation, rigid-body mode treatment]
Wilkinson, J. H. The Algebraic Eigenvalue Problem, Oxford (1965). [reference text for the symmetric generalized eigenvalue problem solved by scipy.linalg.eigh on the dense branch (LAPACK’s DSYGVD)]

class femorph_solver.solvers.modal.ModalResult(omega_sq: ndarray, frequency: ndarray, mode_shapes: ndarray, free_mask: ndarray)[source]#

Bases: object

Result of a modal solve.

omega_sq#

(n_modes,) eigenvalues ω² = (2π f)².

Type:: numpy.ndarray

frequency#

(n_modes,) cyclic frequencies f [Hz].

Type:: numpy.ndarray

mode_shapes#

(N, n_modes) DOF-indexed mode shapes — each column is one mode in the global DOF ordering produced by femorph_solver.model.Model.dof_map().

Type:: numpy.ndarray

free_mask#

(N,) bool — DOFs kept in the eigenproblem after Dirichlet reduction.

Type:: numpy.ndarray

save(path: str | Path, model: Model, *, unit_system: str = 'SI', extra_field_data: dict[str, np.ndarray] | None = None) → Path[source]#

Serialise this modal result to a disk-backed .pv file.

Re-projects the DOF-indexed mode_shapes onto per-node (n_points, n_dof_per_node) arrays using the model’s dof_map(), then hands off to write_modal_result().

Parameters:

path (str or pathlib.Path) – Destination file (.pv extension conventional).
model (femorph_solver.model.Model) – The model whose K/M matrices this result was computed from — supplies the mesh geometry and DOF layout needed to reshape mode_shapes.
extra_field_data (mapping, optional) – Extra field_data entries (e.g. solver statistics).

Returns:

Canonicalised absolute path to the written file.

Return type:

pathlib.Path

Notes

After saving, load the file back via femorph_solver.result.read(); the resulting object is a ModalResult (disk-backed) with the same frequencies and per-mode displacements plus the full plotting / animation surface.

femorph_solver.solvers.modal.solve_modal(K: csr_array, M: csr_array, prescribed: dict[int, float], n_modes: int = 10, *, eigen_solver: str = 'arpack', linear_solver: str = 'auto', sigma: float = 0.0, tol: float = 1e-12, thread_limit: int | None = None, mixed_precision: bool | None = False, v0: ndarray | None = None) → ModalResult[source]#

Solve the generalized eigenproblem K φ = ω² M φ.

Parameters:

K (scipy.sparse.csr_array) – Global stiffness and mass matrices, same shape.
M (scipy.sparse.csr_array) – Global stiffness and mass matrices, same shape.
prescribed (mapping {global_dof_index: 0}) – Indices to drop from the eigenproblem (clamped / zero-stiffness DOFs).
n_modes (int) – Number of lowest modes to return.
eigen_solver (str) – Name of a registered eigen backend ("arpack", "lobpcg", "primme", "dense"), or "auto" to let the solver pick based on problem size. The autotune prefers PRIMME for n_modes > 20 when the package is installed (faster on the clustered eigenvalue spectra typical of rotor modal analyses); it falls back to ARPACK otherwise. See femorph_solver.solvers.eigen.list_eigen_solvers().
linear_solver (str) – Name of the registered linear backend used by ARPACK’s shift-invert step. Pass "auto" (the default) to pick the fastest available SPD direct solver (pardiso → cholmod → umfpack → superlu). See femorph_solver.solvers.linear.list_linear_solvers() for the full registry.
tol (float) – Convergence tolerance passed to the sparse eigensolver. Default 1e-12 — tight enough that the returned frequencies agree with a full machine-epsilon (tol=0) solve to at least 13 significant digits on typical modal problems, while saving ~25 % of ARPACK iterations. Pass 0.0 to force ARPACK to converge to full machine precision; pass a looser value (1e-10, 1e-8) when wall-time matters more than the last couple of digits.
thread_limit (int or None, optional) – Cap on BLAS / OpenMP threads used inside the dense and sparse eigenpaths. None (default) defers to the global limit set via femorph_solver.set_thread_limit() or the FEMORPH_SOLVER_NUM_THREADS environment variable; pass an explicit integer to override for a single call.
v0 (numpy.ndarray or None, optional) – Starting vector for ARPACK’s Lanczos iteration. Default is None — ARPACK seeds from a random vector. Pass a previous mode shape (full-system length N or reduced-system length n_free) to warm-start a parametric sweep — ARPACK’s convergence iterations typically drop 1.5–3× when the new eigenvectors are close to the old ones (e.g. material sensitivity studies, mistuning identification). Ignored by non-ARPACK eigen backends.

femorph_solver.solvers.modal.solve_modal_reduced(K_ff: csr_array, M_ff: csr_array, free_mask: ndarray, n_modes: int = 10, *, eigen_solver: str = 'arpack', linear_solver: str = 'auto', sigma: float = 0.0, tol: float = 1e-12, thread_limit: int | None = None, mixed_precision: bool | None = False, v0: ndarray | None = None) → ModalResult[source]#

Solve the eigenproblem directly on pre-reduced K_ff / M_ff.

Skips the full-matrix BC slicing that solve_modal() does, so the caller can build K_ff / M_ff via _bc_reduce(), release the original K / M (Model cache + their own refs), then call this to run the eigensolve with a peak RSS roughly K_ff + M_ff + triu + MKL factor instead of the default path’s K + M + K_ff + M_ff + triu + MKL factor.

Cyclic-symmetry modal analysis for a single base sector.

A rotor with N identical sectors spanning 360° satisfies, at each harmonic index (nodal diameter) k ∈ {0, 1, …, N/2},

u_sector_{j+1}(global) = e^{i k α} · R(α) · u_sector_j(global)

with α = 2π/N and R(α) the spatial rotation about the symmetry axis. Applied at the base-sector’s high-angle face (which is identified with sector 1’s low-angle face), this gives the constraint

u_high(global) = e^{i k α} · R(α) · u_low(global)

relating paired low-/high-face DOFs. That constraint lets a modal solve work on one sector and produce the full-rotor spectrum one k at a time — an N-times speedup relative to meshing the whole rotor.

The reduction here matches MAPDL’s CYCLIC approach in its complex form. Base-sector DOFs are partitioned into interior I and low-face L / high-face H sets, with a one-to-one (low-node, high-node) correspondence given by the user. High DOFs are eliminated via the constraint above; the remaining [u_I; u_L] vector satisfies a Hermitian generalised eigenproblem

K_red φ = ω² M_red φ

with K_red = P^H K P and M_red = P^H M P (K, M real SPD; P the complex constraint matrix that embeds the reduced vector back into the full one). For k=0 and k=N/2 (even N) the phase is ±1 and the reduced problem can be real; for intermediate k it is complex Hermitian.

Per-pair rotation#

Stiffness / mass are usually assembled in a global XYZ frame, which means the cyclic phase relation picks up the same rotation that maps sector j to sector j+1. Pass pair_rotation as the (d, d) matrix R(α) acting on each node pair’s d translational DOFs (3 for a 3D solid sweeping about z). Omit it (None) when face DOFs are already in a cylindrical / per-sector local frame, in which case R(α) reduces to the identity.

Harmonic-index counting#

k=0 — single real eigenmodes (in-phase family)
k=N/2 (even N) — single real eigenmodes (anti-phase family)
0 < k < N/2 — each eigenvalue corresponds to a degenerate pair in the full rotor (travelling-wave directions). The complex reduced solve returns them once per k; the full-rotor multiset is recovered by counting them twice.

References

Thomas, D. L. “Dynamics of rotationally periodic structures.” International Journal for Numerical Methods in Engineering 14 (1), 81-102 (1979). [the foundational complex-valued cyclic reduction applied here — one base sector factored at each harmonic index k recovers the full-rotor spectrum]
Wildheim, S. J. “Vibrations of rotationally periodic structures.” Ph.D. dissertation, Linköping University (1979). [companion treatment of the constraint matrix P and its Hermitian reduction P^H K P]
Bathe, K.-J. Finite Element Procedures, Prentice Hall (1996), §10.3.4 (cyclic symmetry). [textbook derivation matching the partitioned [u_I; u_L] variable convention used below]

class femorph_solver.solvers.cyclic.CyclicModalResult(harmonic_index: int, n_sectors: int, omega_sq: ndarray, frequency: ndarray, mode_shapes: ndarray, axis_point: ndarray | None = None, axis_dir: ndarray | None = None)[source]#

Bases: object

Modes for one harmonic index of a cyclic-symmetry modal solve.

axis_dir: ndarray | None = None#: (3,) — unit vector along the rotor axis. None means “not recorded”; writers default to +Z. Populated by solve_cyclic_modal() from the supplied pair_rotation.

axis_point: ndarray | None = None#: (3,) — point on the rotor rotation axis. None means “not recorded by the caller”; writers default to the origin in that case. Populated by solve_cyclic_modal() when a pair_rotation is supplied.

frequency: ndarray#: (n_modes,) — cyclic frequencies f [Hz].

mode_shapes: ndarray#: (N, n_modes) — complex mode shapes on the base sector, indexed by the original DOF layout (not the reduced layout). For k=0 and k=N/2 the imaginary part is identically zero.

n_sectors: int#: Number of sectors in the full rotor (same for every result in a sweep).

omega_sq: ndarray#: (n_modes,) — ω² = (2πf)²; numerically negative values clipped to 0.

femorph_solver.solvers.cyclic.save_cyclic_modal_results(results: Sequence[CyclicModalResult], path: str | Path, model: Model, *, unit_system: str = 'SI', extra_field_data: dict[str, np.ndarray] | None = None) → Path[source]#

Serialise a cyclic-modal sweep to a single disk-backed .pv file.

Parameters:

results (sequence of CyclicModalResult) – Per-harmonic-index results as returned by solve_cyclic_modal(). All entries must share the same n_sectors.
path (str or pathlib.Path) – Destination file.
model (femorph_solver.model.Model) – Model whose K / M matrices this sweep was computed on; supplies the mesh and DOF layout.
extra_field_data (mapping, optional) – Extra field_data entries (e.g. solver statistics).

Returns:

Canonicalised absolute path to the written file.

Return type:

pathlib.Path

Notes

The DOF-indexed (N, n_modes) complex mode shapes on each per-k result get re-projected onto per-(k, mode), per-node (n_modes, n_points, 6) complex arrays — the layout write_cyclic_modal_result() expects.

femorph_solver.solvers.cyclic.solve_cyclic_modal(K: spmatrix | sparray, M: spmatrix | sparray, low_node_dofs: ndarray, high_node_dofs: ndarray, n_sectors: int, n_modes: int = 10, *, pair_rotation: ndarray | None = None, harmonic_indices: Sequence[int] | None = None, prescribed: dict[int, float] | None = None, thread_limit: int | None = None, linear_solver: str = 'auto') → list[CyclicModalResult][source]#

Solve the cyclic-symmetry modal problem on a single base sector.

Parameters:

K (scipy.sparse) – Global stiffness / mass of the base sector, before cyclic or Dirichlet reduction. Assumed real SPD.
M (scipy.sparse) – Global stiffness / mass of the base sector, before cyclic or Dirichlet reduction. Assumed real SPD.
low_node_dofs ((P, d) int arrays) – Global DOF indices at P paired boundary nodes, each with d DOFs (typically 3 for a 3D solid). low_node_dofs[i, :] and high_node_dofs[i, :] must address the same physical DOFs (same labels, same order) at the node pair that’s identified by the cyclic BC. A 1-D array is accepted for scalar (d=1) problems.
high_node_dofs ((P, d) int arrays) – Global DOF indices at P paired boundary nodes, each with d DOFs (typically 3 for a 3D solid). low_node_dofs[i, :] and high_node_dofs[i, :] must address the same physical DOFs (same labels, same order) at the node pair that’s identified by the cyclic BC. A 1-D array is accepted for scalar (d=1) problems.
n_sectors (int) – Number of repetitions that close the rotor (N).
n_modes (int, default 10) – Number of lowest modes per harmonic index.
pair_rotation ((d, d) array, optional) – Spatial rotation R(α) that maps sector j’s local frame to sector j+1’s, applied per node pair. Pass None when face DOFs are already in a per-sector local frame (then R=I implied).
harmonic_indices (sequence of int, optional) – Harmonic indices to solve. Defaults to 0, 1, …, N//2.
prescribed (mapping, optional) – Dirichlet BCs on the sector (global DOF index → 0). Currently only zero-value Dirichlet is supported; if a cyclic-face DOF is prescribed its partner must also be prescribed (so the cyclic constraint is trivially satisfied at that DOF).
thread_limit (int or None, optional) – Cap on BLAS / OpenMP threads used inside the per-harmonic dense Hermitian eigensolve. None defers to the global limit (see femorph_solver.set_thread_limit()).
linear_solver (str, default "auto") – Name of the registered linear backend used for the inner shift-invert factorisation in the sparse path. "auto" picks the fastest available SPD direct solver in the priority chain pardiso → cholmod → umfpack → superlu and emits a one-shot warning if pardiso is missing on a problem large enough to benefit from it. See femorph_solver.solvers.linear.list_linear_solvers().

Returns:

One CyclicModalResult per requested harmonic index, in the order requested.

Return type:

list[CyclicModalResult]

Elements#

femorph-solver element library.

Each element class implements ElementBase and registers itself via register(). The registry is how the assembler (and anything else walking a mesh) dispatches per-element math from the ansys_elem_type_num cell-data array on a model’s grid.

class femorph_solver.elements.ElementBase[source]#

Bases: ABC

abstractmethod static ke(coords: ndarray, material: dict[str, float], real: ndarray) → ndarray[source]#

Return the element stiffness matrix in global DOF ordering.

Parameters:

coords ((n_nodes, 3) float64)
material (mapping with MAPDL property names as keys (EX, PRXY, DENS, ...))
real (1-D array of real constants)

abstractmethod static me(coords: ndarray, material: dict[str, float], real: ndarray, lumped: bool = False) → ndarray[source]#: Return the element mass matrix in global DOF ordering.

MAPDL compatibility#

MAPDL compatibility layer.

This subpackage isolates everything MAPDL-specific in femorph-solver:

femorph_solver.mapdl_api.io — readers and writers for MAPDL binary formats (RST / FULL / EMAT / HBMAT) and their framing.
femorph_solver.mapdl_api.from_cdb() — load a MAPDL CDB input deck.

The core femorph-solver library (femorph_solver.Model, elements, solvers, generic VTK helpers) has no dependency on anything in this subpackage. Install the optional mapdl extra to pull in mapdl_archive, the only external dependency required by this layer:

pip install femorph_solver[mapdl]

femorph_solver.mapdl_api.from_cdb(path: str, **kwargs: Any) → Model[source]#

Load a MAPDL CDB file via mapdl_archive.

Requires the mapdl extra: pip install femorph_solver[mapdl].

The deck’s /UNITS command (if present) is parsed and stamped onto Model.unit_system. Decks without a /UNITS line get UnitSystem.UNSPECIFIED.

Load MAPDL CDB input decks into an femorph-solver Model.

Thin wrapper around mapdl_archive. Only the file-format reader is involved — MAPDL itself is never invoked.

The CDB’s ET table (MAPDL element-type declarations such as ET, 1, 186) is translated on load so the returned Model already knows each ET id’s kernel — callers don’t have to re-declare them.

References

CDB format definition (MAPDL input deck, produced by the CDWRITE command): Ansys, Inc., Ansys Mechanical APDL Command Reference, Release 2022R2, entry CDWRITE; and Ansys Mechanical APDL Basic Analysis Guide, chapter on the database file.
Open-source parser: mapdl-archive (akaszynski/mapdl-archive, MIT) — the production parser we delegate to. Does not invoke MAPDL, only reads the text grammar.

femorph_solver.mapdl_api.cdb.from_cdb(path: str, **kwargs: Any) → Model[source]#

Load a MAPDL CDB file via mapdl_archive.

Requires the mapdl extra: pip install femorph_solver[mapdl].

The deck’s /UNITS command (if present) is parsed and stamped onto Model.unit_system. Decks without a /UNITS line get UnitSystem.UNSPECIFIED.

Binary IO (MAPDL formats)#

Low-level reader for MAPDL binary files (RST / FULL / EMAT / …).

MAPDL binary files use Fortran-style sequential records. Each record is laid out as:

[ size:int32 ] [ flags:int32 ] [ payload : size * 4 bytes ] [ pad:int32 ]

size counts 4-byte words of payload (so a double array of length n gives size = 2 * n).
flags only uses its least-significant byte (byte 7 from the start of the record header). femorph-solver exposes the raw byte but does not use it to infer payload dtype or compression state. pymapdl-reader’s Python and C++ backends disagree on the bit positions and neither matches files produced by recent MAPDL — e.g. bit 3 (declared “sparse” by both) is routinely set on perfectly dense records emitted under /FCOMP,RST,0. Callers must pass the dtype appropriate to the record they expect (int32 for index tables, float64 for coordinates and results, etc.).
the trailing pad int32 is a Fortran record-length repeater.

femorph-solver’s decoder supports uncompressed, dense records in the four standard dtypes (int32, int16, float64, float32). Compression (bsparse, wsparse, zlib) is not attempted — run MAPDL with /FCOMP,<file>,0 to disable.

class femorph_solver.mapdl_api.io.binary.RecordInfo(size: int, dtype: dtype, flags: int, offset: int)[source]#

Bases: object

Metadata for a single MAPDL record.

femorph_solver.mapdl_api.io.binary.decompress_bsparse(payload: ndarray, dtype: dtype | str = '<f8') → ndarray[source]#

Decompress a bit-sparse (bsparse) record payload.

MAPDL’s bsparse encoding for a vector of length N:

word[0]  = N                     # decompressed length
word[1]  = bitcode               # 32-bit mask; bit i set -> position i nonzero
payload  = K packed values       # K = number of set bits in bitcode

Output is a dense ndarray of length N with zeros in unset slots. Only N <= 32 is supported — MAPDL uses a different encoding for longer vectors.

femorph_solver.mapdl_api.io.binary.decompress_wsparse(payload: ndarray, dtype: dtype | str = '<f8') → ndarray[source]#

Decompress a windowed-sparse (wsparse) record payload.

MAPDL’s wsparse encoding for a vector of length N:

word[0]  = N           # decompressed length
word[1]  = NWin        # number of windows
for each window:
    iLoc = int32
    if iLoc > 0:       # isolated single value at position iLoc
        one value (dtype-sized)
    else:              # window starts at position -iLoc
        iLen = int32
        if iLen > 0:   # iLen explicit values
            iLen values
        else:          # constant run of (-iLen) copies
            one value

All missing positions are zero.

femorph_solver.mapdl_api.io.binary.iter_records(fp: BinaryIO, start: int = 0, dtype: dtype | str = '<i4')[source]#

Yield (offset_bytes, RecordInfo, ndarray) tuples sequentially.

Stops at EOF or when the size field is non-positive (MAPDL marks the end of the sequential stream with a -1 sentinel or zero padding).

femorph_solver.mapdl_api.io.binary.open_for_write(path: str | Path) → BinaryIO[source]#: Open path for writing as a new MAPDL binary file.

femorph_solver.mapdl_api.io.binary.read_record(fp: BinaryIO, offset: int | None = None, dtype: dtype | str = '<i4') → tuple[ndarray, RecordInfo][source]#

Read one MAPDL record starting at offset (bytes) or current position.

dtype selects how to cast the payload. MAPDL doesn’t encode dtype reliably in the flag byte, so callers must pass the dtype appropriate to the record they expect (int32 for index tables, float64 for coordinates and results, etc.). The default is little-endian int32.

Returns (array, info) where array is the dense payload and info carries the header metadata. Advances the file pointer past the trailing Fortran pad word.

femorph_solver.mapdl_api.io.binary.read_record_at(path: str | Path, word_offset: int, dtype: dtype | str = '<i4') → tuple[ndarray, RecordInfo][source]#: Convenience: open, seek to a 4-byte-word offset, read one record.

femorph_solver.mapdl_api.io.binary.read_standard_header(path: str | Path) → dict[source]#

Parse the MAPDL standard file header.

Every MAPDL binary file begins with a 100-word standard header whose fields live at fixed int32 word offsets. Numeric fields occupy one int32 word followed by a single padding word; string fields are packed 4 chars per word with the bytes of each word reversed (MAPDL stores them with the Fortran byte order flipped).

The returned mapping covers the canonical fields: file_number, file_format, time, date, units, verstring, mainver, subver, machine, jobname, product, username, machine_identifier, system_record_size, jobname2, title, subtitle, split_point.

femorph_solver.mapdl_api.io.binary.write_record(fp: BinaryIO, payload: ndarray, flags: int | None = None) → int[source]#

Append one MAPDL record to fp and return its byte offset.

payload is written verbatim as bytes. size is derived from len(payload.tobytes()) // 4 so any numpy dtype whose itemsize is a multiple of 4 works. flags defaults to 0x80 for int32 payloads and 0x00 for float64; any other dtype defaults to 0x00.

femorph_solver.mapdl_api.io.binary.write_standard_header(fp: BinaryIO, *, file_format: int, jobname: str = 'file', time: str = '00:00:00', date: str = '', units: int = 0, verstring: str = '22.2', machine: str = 'LINUX x64', product: str = 'FULL', username: str = '', title: str = '') → None[source]#

Write the fixed-layout standard-header record to fp.

Mirrors read_standard_header(). The record holds a 100-word payload; the payload indices below are payload-relative. Read side sees these words at file-absolute positions [i+2] because the record framer prepends the 2-word (size, flags) header.

Typed record schemas for MAPDL RST / FULL / EMAT files.

The key lists mirror those in fdresu.inc (the Fortran include file that defines the layout of result-file records). They are reproduced verbatim from pymapdl-reader’s _rst_keys.py.

parse_header maps a flat int32 table onto a named dict, handling two conventions that MAPDL uses inside these tables:

Array fields — keys that appear multiple times (e.g. DOFS, title) accumulate into a list.
Split 64-bit pointers — pairs named ptrXXXl / ptrXXXh are reassembled into a single 64-bit ptrXXX value.

References

Ansys MAPDL result-file layout (primary compatibility source): Ansys, Inc., Ansys Mechanical APDL Theory Reference, Release 2022R2 — binary-file chapter describing the block / record / pointer structure of .rst, .full, and .emat.
Open-source cross-reference used to port the Fortran record keys into Python: pymapdl-reader (MIT licence, ansys/pymapdl-reader) — the key tables here (solution_header_keys, result_header_keys, …) are reproduced from its _rst_keys.py. The underlying layout originates in MAPDL’s own fdresu.inc Fortran include file.

femorph_solver.mapdl_api.io.rst_schema.find_record(path: str | Path, index: int, dtype: dtype | str = '<i4') → ndarray[source]#: Return the payload of record index (0-based) as an ndarray.

femorph_solver.mapdl_api.io.rst_schema.parse_header(table: ndarray, keys: list[str]) → dict[source]#

Map a flat int32 table onto a named header dict.

Array-valued keys (those that appear more than once in keys) are gathered into a list. Pairs ptrXXXl / ptrXXXh are reassembled into a single 64-bit ptrXXX entry.

femorph_solver.mapdl_api.io.rst_schema.read_nodal_displacement(path: str | Path, set_index: int = 0) → ndarray[source]#

Return the decompressed nodal displacement vector for result set 0.

Output shape is (nnod, numdof). Entries for DOFs with no stored value (bits unset in the bsparse bitmask) come back as zero.

femorph_solver.mapdl_api.io.rst_schema.read_rst_geometry_header(path: str | Path) → dict[source]#

Read the geometry header from an RST file using ptrGEO.

The result header’s ptrGEO is a byte-offset expressed in 4-byte words from the start of the file. Seek there and read one record.

femorph_solver.mapdl_api.io.rst_schema.read_rst_result_header(path: str | Path) → dict[source]#: Read record #1 of an RST file and parse it as a result header.

femorph_solver.mapdl_api.io.rst_schema.read_solution_data_header(path: str | Path, set_index: int = 0) → tuple[dict, int][source]#

Parse the solution-data header for result set set_index.

Returns (header_dict, solution_word_offset). The word offset is the absolute file position of the solution header’s 8-byte record frame — useful because ptrNSL / ptrESL / ptrBC inside the header are relative to that position.

Schema + reader for MAPDL .full files (assembled K/M/C matrices).

A FULL file stores, for each column of the upper triangle, two framed records: one holding the 1-based row indices of that column’s nonzero entries, and one holding the matching double-precision values. The record pointers ptrSTF / ptrMAS / ptrDMP address those column blocks; ptrDOF addresses a (ndof, const) pair of records that describes which equations are constrained.

Layout sketch (for one column c of an ntermK-nonzero stiffness matrix, upper triangle only):

[ size=nitems | flags | rows(nitems int32, 1-based) | pad4 ]
[ size=2·nitems | flags | vals(nitems doubles) | pad4 ]

The constrained-DOF mask (const[i] < 0) is applied at read time to drop fixed rows and columns, matching pymapdl-reader’s behavior.

References

Ansys MAPDL .full file layout (compatibility source): Ansys, Inc., Ansys Mechanical APDL Theory Reference, Release 2022R2 — binary-file chapter, section on .full (assembled global K / M / C matrices).
Open-source cross-reference: pymapdl-reader (MIT, ansys/pymapdl-reader) full_file.py — used here only as a key-map yardstick; the underlying CSC-like column-indexed layout comes from MAPDL’s own Fortran code.

femorph_solver.mapdl_api.io.full_schema.read_full_header(path: str | Path) → dict[source]#: Read the FULL file header (record index 1) as a named dict.

femorph_solver.mapdl_api.io.full_schema.read_full_km(path: str | Path, reduce: bool = True) → tuple[csc_matrix | None, csc_matrix | None, ndarray][source]#

Read stiffness + mass matrices from a MAPDL FULL file.

Parameters:

path – Path to the .full file.
reduce – If True (the default), drop rows and columns whose const[i] < 0 (i.e. the constrained DOFs MAPDL still includes in the assembled matrix but which the user almost always wants eliminated). This matches the behavior of pymapdl_reader.FullFile.load_km.

Returns:

k and m are scipy CSC matrices (or None if the file doesn’t store that matrix). dof_ref is an (neqn_out, 2) int32 array of (node, dof_id) pairs aligned with the matrix rows; dof_id uses the 1-based MAPDL DOF numbering (see femorph_solver.mapdl_api.io.rst_schema.DOF_REF).

Return type:

k, m, dof_ref

Schema + reader for MAPDL .emat (element matrix) files.

An EMAT file stores the per-element stiffness, mass, damping, stress- stiffening, applied-force, Newton–Raphson, and imaginary-load contributions that MAPDL summed to build the global K/M/C matrices in the matching FULL file. Layout (per pymapdl-reader’s emat.py and fdemat.inc):

Record 0: standard header
Record 1 (at word 103): EMAT header (40 words; see emat_header_keys)
Record 2: time-info record (20 doubles)
Record at ptrDOF: DOF identifier list (numdof ints)
Record at ptrBAC: nodal equivalence (lenbac ints)
Record at ptrElm: element equivalence (nume ints)
Record at ptrBIT: DOF bits (lenu ints)
Record at ptrIDX: element-matrix index table (2 * nume ints, (lo, hi) pairs pointing to each element’s block)

Per-element block (repeated nume times, addressed via the index table):

[ element matrix header (10 ints) ]
[ DOF index table (nmrow ints, 1-based global DOF numbers) ]
[ K matrix     if stkey ]  ← symmetric stored lower triangular
[ M matrix     if mkey  ]
[ D matrix     if dkey  ]
[ SS matrix    if sskey ]
[ applied force (2*nmrow doubles) if akey ]
[ NR load       (nmrow doubles) if nrkey ]
[ imaginary    (nmrow doubles) if ikey  ]

Symmetric element matrices are packed lower-triangular in column-major order (n(n+1)/2 doubles). If nmrow is positive the matrix is stored unsymmetric full (n*n doubles).

References

Ansys MAPDL element-matrix file layout (compatibility source): Ansys, Inc., Ansys Mechanical APDL Theory Reference, Release 2022R2 — binary-file chapter, section on .emat.
Fortran include reference used as the source for key names and record ordering: MAPDL’s internal fdemat.inc.
Open-source cross-reference: pymapdl-reader (MIT, ansys/pymapdl-reader) emat.py. We rely on it as a compatibility yardstick for record offsets, not for any algorithmic content.

femorph_solver.mapdl_api.io.emat_schema.read_emat_element(path: str | Path, index: int) → tuple[ndarray, dict[str, ndarray]][source]#

Read one element’s block from an EMAT file.

Parameters:

path – Path to the .emat file.
index – 0-based element index (into the element equivalence table; not the MAPDL element number unless the file is already sorted).

Returns:

dof_idx is a zero-based array of the global DOF positions this element contributes to ((N-1)*numdof + dof). matrices maps the keys {'k', 'm', 'd', 'ss', 'applied_force', 'newton_raphson', 'imaginary_load'} to dense (nmrow, nmrow) matrices (for K/M/D/SS) or length-nmrow / length-2·nmrow vectors (for the force terms), whichever are present in the file.

Return type:

dof_idx, matrices

femorph_solver.mapdl_api.io.emat_schema.read_emat_element_matrices(path: str | Path) → list[tuple[ndarray, dict[str, ndarray]]][source]#: Read every element’s (dof_idx, matrices) from the EMAT file.

femorph_solver.mapdl_api.io.emat_schema.read_emat_header(path: str | Path) → dict[source]#: Read the EMAT header record (word 103) as a named dict.