"""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 :mod:`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]
"""
from __future__ import annotations
from dataclasses import dataclass
from pathlib import Path
from typing import TYPE_CHECKING
import numpy as np
import scipy.sparse as sp
from .._threads import limit as _thread_limit
from .linear import get_linear_solver
if TYPE_CHECKING:
from femorph_solver.model import Model
[docs]
@dataclass(frozen=True, slots=True)
class StaticResult:
"""Result of a linear static solve.
Attributes
----------
displacement : numpy.ndarray
``(N,)`` DOF-indexed displacement produced by the solve.
reaction : numpy.ndarray
``(N,)`` reaction force; nonzero only at constrained DOFs.
free_mask : numpy.ndarray
``(N,)`` bool — ``True`` at solved-for DOFs.
"""
displacement: np.ndarray
reaction: np.ndarray
free_mask: np.ndarray
def __repr__(self) -> str:
# The default dataclass repr dumps full ndarrays — useless on
# any realistic mesh (10-100 k DOFs). Compact summary covers
# the same shape signal without spamming the terminal.
n = int(self.displacement.shape[0])
n_free = int(self.free_mask.sum())
max_u = float(np.max(np.abs(self.displacement))) if n else 0.0
return f"StaticResult(N={n}, free={n_free}, max|u|={max_u:.3e})"
[docs]
def elastic_strain(self, *, model: Model) -> np.ndarray:
"""Recover ``(n_points, 6)`` Voigt elastic strain at every grid node.
Strain is evaluated at each element's own nodes by
:math:`\\varepsilon = B(\\xi_\\text{node})\\,u_e` and averaged
across every element that touches a shared node.
Parameters
----------
model : femorph_solver.model.Model
Model the static solve ran on — supplies the element-type
registry needed to recover element-wise strain.
Returns
-------
numpy.ndarray
``(n_points, 6)`` engineering-shear Voigt strain
``[ε_xx, ε_yy, ε_zz, γ_xy, γ_yz, γ_xz]``.
"""
from femorph_solver.result._stress_recovery import compute_nodal_strain
return compute_nodal_strain(model, self.displacement)
[docs]
def elastic_strain_per_element(self, *, model: Model) -> dict[int, np.ndarray]:
"""Per-element elastic strain at every element node — no averaging.
Returns the raw ``{elem_num: (n_nodes_in_elem, 6)}`` dict that
:meth:`elastic_strain` averages. Useful when a verification
check needs to compare element-by-element values rather than
the smoothed nodal-averaged field.
"""
return model.strain(self.displacement)
[docs]
def plot(
self,
model: Model,
*,
comp: str = "norm",
displacement_factor: float = 1.0,
overlay_wireframe: bool = True,
**plot_kwargs: object,
) -> object:
"""Plot the deformed shape coloured by displacement.
Convenience wrapper that projects the DOF-indexed
:attr:`displacement` onto per-node UX / UY / UZ and hands off
to the existing
:func:`femorph_solver.plotting.plot_displacement_snapshot`
helper.
Parameters
----------
model : femorph_solver.Model
Model whose K matrix this result was computed on; supplies
the mesh + DOF layout.
comp : str, default ``"norm"``
Scalar to colour by — ``"norm"`` (magnitude), ``"ux"``,
``"uy"``, or ``"uz"``.
displacement_factor : float, default 1.0
Scale applied to the mesh deformation. Pass a larger
value (e.g. 100.0) to make small elastic deformations
visible at the mesh's geometric scale.
overlay_wireframe : bool, default True
Overlay the undeformed mesh as a faint wireframe.
**plot_kwargs
Forwarded to
:func:`femorph_solver.plotting.plot_displacement_snapshot`
(and through to ``pyvista.Plotter`` — accepts ``off_screen``,
``screenshot``, ``cmap``, ``window_size``, ``cpos``, …).
Returns
-------
Any
Whatever ``pyvista.Plotter.show()`` returns (typically
``None`` in interactive mode, or the camera position
tuple when ``off_screen=True``).
See Also
--------
femorph_solver.plotting.plot_displacement_snapshot : the
underlying helper with full control over the warp + scene.
save : persist as a disk-backed StaticResultDisk that has its
own ``.plot_nodal_displacement()`` method (no model arg
required since the disk handle embeds the grid).
Examples
--------
Requires a working OpenGL / OSMesa render context — every
line skipped under doctest gating because headless CI runners
typically lack one::
import femorph_solver as fs
model = fs.templates.cantilever_hex8()
result = model.solve_static()
result.plot(model, off_screen=True, displacement_factor=1e6)
"""
from femorph_solver.plotting._displacement import ( # noqa: PLC0415
plot_displacement_snapshot,
)
# Project DOF-indexed displacement onto per-node (UX, UY, UZ).
dof_map = model.dof_map()
node_nums = np.asarray(model.grid.point_data["ansys_node_num"], dtype=np.int64)
node_to_row = {int(n): i for i, n in enumerate(node_nums)}
node_row = np.fromiter(
(node_to_row[int(n)] for n in dof_map[:, 0]),
dtype=np.int64,
count=dof_map.shape[0],
)
comp_idx = dof_map[:, 1].astype(np.int64)
# comp_idx 0/1/2 = UX/UY/UZ; ignore rotational DOFs (3-5).
translational_mask = comp_idx < 3
displacement_xyz = np.zeros((model.grid.n_points, 3), dtype=np.float64)
rows = node_row[translational_mask]
cols = comp_idx[translational_mask]
displacement_xyz[rows, cols] = self.displacement[translational_mask]
return plot_displacement_snapshot(
model.grid,
displacement_xyz,
comp=comp,
displacement_factor=displacement_factor,
overlay_wireframe=overlay_wireframe,
**plot_kwargs,
)
[docs]
def save(
self,
path: str | Path,
model: Model,
*,
unit_system: str = "SI",
extra_field_data: dict[str, np.ndarray] | None = None,
) -> Path:
"""Serialise this static result to a disk-backed ``.pv`` file.
Re-projects the DOF-indexed :attr:`displacement` (and
:attr:`reaction` as the force field) onto per-node
``(n_points, n_dof_per_node)`` arrays using the model's
:meth:`~femorph_solver.model.Model.dof_map`, then hands off
to :func:`~femorph_solver.result.static.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
-------
pathlib.Path
Canonicalised absolute path to the written file.
"""
from femorph_solver.result._material import MaterialTable
from femorph_solver.result.static import write_static_result
grid = model.grid
dof_map = model.dof_map()
node_nums = np.asarray(grid.point_data["ansys_node_num"], dtype=np.int64)
node_to_row = {int(n): i for i, n in enumerate(node_nums)}
node_row = np.fromiter(
(node_to_row[int(n)] for n in dof_map[:, 0]),
dtype=np.int64,
count=dof_map.shape[0],
)
comp = dof_map[:, 1].astype(np.int64)
disp = np.zeros((grid.n_points, 6), dtype=np.float64)
disp[node_row, comp] = self.displacement
force = np.zeros((grid.n_points, 6), dtype=np.float64)
force[node_row, comp] = self.reaction
material_table = MaterialTable.from_model(model, unit_system=str(unit_system))
return write_static_result(
path,
grid,
disp,
force=force,
material_table=material_table,
unit_system=unit_system,
extra_field_data=extra_field_data,
)
def _build_coupling_transform(
n_dof: int, coupling: list[tuple[int, int]]
) -> tuple[sp.csr_array, np.ndarray]:
"""Build the ``(N, N_eff)`` master-slave transformation ``T``.
Each ``(slave, master)`` pair forces ``u[slave] = u[master]`` after
the reduction by mapping the slave row in ``T`` to the master's
column. Returns ``(T, eff_idx)`` where ``eff_idx[i]`` is the
reduced-space column for full DOF ``i`` (slaves share their
master's column).
A coupling pair containing a previously-introduced slave (chained
coupling) is normalised to point at the chain's root master so
the resulting ``T`` is rank-correct and there are no transitive
surprises during ``T.T @ K @ T``.
"""
# Resolve transitive masters first: walk each slave's chain until
# we hit a non-slave so the eff_idx assignment below is canonical.
slave_to_master: dict[int, int] = {}
for s, m in coupling:
slave_to_master[int(s)] = int(m)
root_master: dict[int, int] = {}
for s in slave_to_master:
cur = s
while cur in slave_to_master and slave_to_master[cur] != cur:
cur = slave_to_master[cur]
root_master[s] = cur
eff_idx = np.full(n_dof, -1, dtype=np.int64)
counter = 0
for i in range(n_dof):
if i in root_master and root_master[i] != i:
continue
eff_idx[i] = counter
counter += 1
for s, m in root_master.items():
eff_idx[s] = eff_idx[m]
n_eff = int(counter)
rows = np.arange(n_dof, dtype=np.int64)
cols = eff_idx
data = np.ones(n_dof, dtype=np.float64)
transform = sp.csr_array((data, (rows, cols)), shape=(n_dof, n_eff))
return transform, eff_idx
[docs]
def solve_static(
K: sp.csr_array,
F: np.ndarray,
prescribed: dict[int, float],
*,
coupling: list[tuple[int, int]] | None = None,
linear_solver: str = "auto",
thread_limit: int | None = None,
mixed_precision: bool | None = False,
) -> StaticResult:
"""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}``
coupling : list of ``(slave, master)`` pairs, optional
Multi-point constraints (MAPDL ``CP``). Each pair forces
``u[slave] = u[master]`` after solve. Implemented via the
canonical master-slave transformation
``K_eff = Tᵀ K T``, ``F_eff = Tᵀ F`` where ``T`` maps the
reduced-DOF vector back to the full one (slaves get the
master's column).
linear_solver : str
Name of a registered linear backend — see
:func:`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
:func:`femorph_solver.set_thread_limit`).
"""
n_dof = K.shape[0]
if coupling:
transform, _eff_idx = _build_coupling_transform(n_dof, coupling)
# Prescribed DOFs flow through the transform: a prescribed
# slave constrains the master to the same value. Translate
# each prescribed full-DOF entry into a reduced-DOF entry.
slaves = {int(s) for s, _ in coupling}
eff_prescribed: dict[int, float] = {}
for idx, val in prescribed.items():
eff_idx_i = int(_eff_idx[int(idx)])
existing = eff_prescribed.get(eff_idx_i)
if existing is not None and not np.isclose(existing, val):
raise ValueError(
f"coupling: prescribed DOFs {idx} (= {val}) and its master "
f"are inconsistent (already pinned to {existing})"
)
eff_prescribed[eff_idx_i] = float(val)
_ = slaves # silence linter; the slaves set is used implicitly via _eff_idx
k_red = (transform.T @ K @ transform).tocsr()
f_red = transform.T @ F
result_eff = solve_static(
k_red,
f_red,
eff_prescribed,
linear_solver=linear_solver,
thread_limit=thread_limit,
)
# Scatter the reduced solution back to the full DOF space.
u_full = transform @ result_eff.displacement
# Reactions: in the reduced space, the master holds the
# combined reaction of the master + slave DOFs; replicate it
# onto every coupled slave so a downstream consumer sees the
# right physical force at each node. Free mask is rebuilt
# against the full-DOF system so callers see the full layout.
r_full = transform @ result_eff.reaction
free_full = transform @ result_eff.free_mask.astype(np.float64)
return StaticResult(
displacement=u_full,
reaction=r_full,
free_mask=(free_full > 0.5),
)
n_dof_full = n_dof
u = np.zeros(n_dof_full, dtype=np.float64)
mask_fixed = np.zeros(n_dof_full, dtype=bool)
for idx, val in prescribed.items():
mask_fixed[int(idx)] = True
u[int(idx)] = float(val)
diag = K.diagonal()
scale = float(np.max(np.abs(diag))) if diag.size else 1.0
tol = 1e-12 * scale if scale > 0 else 1e-30
mask_fixed |= np.abs(diag) <= tol
free = ~mask_fixed
K_csc = K.tocsc()
K_ff = K_csc[free][:, free]
K_fc = K_csc[free][:, mask_fixed]
rhs = F[free] - K_fc @ u[mask_fixed]
with _thread_limit(thread_limit):
if free.any():
Solver = get_linear_solver(linear_solver, n_dof=K_ff.shape[0])
try:
lu = Solver(K_ff, assume_spd=True, mixed_precision=mixed_precision)
u[free] = lu.solve(rhs)
except Exception:
# Factor failed — usually a deck-side gap where
# implicit BCs were assumed but never set (vm21 / vm71
# 2D-in-3D BEAM models leave UZ+ROTY unanchored). Try
# the null-space pin recovery on small problems.
if K_ff.shape[0] > 2000:
raise
# Build an "avoid" mask for the auto-pin: DOFs the
# user has loaded (force ≠ 0) or prescribed
# (Dirichlet) shouldn't be pinned away.
free_idx_local = np.where(free)[0]
avoid = np.zeros(K_ff.shape[0], dtype=bool)
F_free = F[free_idx_local]
avoid |= np.abs(F_free) > 0.0
rb_pins = _detect_rigid_body_pins(K_ff, avoid=avoid)
if not rb_pins:
raise
free_idx = np.where(free)[0]
for local_i in rb_pins:
mask_fixed[free_idx[local_i]] = True
free = ~mask_fixed
K_ff = K_csc[free][:, free]
K_fc = K_csc[free][:, mask_fixed]
rhs = F[free] - K_fc @ u[mask_fixed]
Solver = get_linear_solver(linear_solver, n_dof=K_ff.shape[0])
lu = Solver(K_ff, assume_spd=True, mixed_precision=mixed_precision)
u[free] = lu.solve(rhs)
R = np.zeros(n_dof_full, dtype=np.float64)
if mask_fixed.any():
R[mask_fixed] = K_csc[mask_fixed] @ u - F[mask_fixed]
return StaticResult(displacement=u, reaction=R, free_mask=free)
def _detect_rigid_body_pins(K_ff, avoid: np.ndarray | None = None) -> list[int]: # type: ignore[no-untyped-def]
"""Identify free-DOF rows that lie in K_ff's null space.
Returns a list of local indices (within ``K_ff``) to pin so the
remaining sub-block is invertible. Walks the dense null-space
via SVD — fine for the small VM-corpus problems where K_ff is
typically well below 1k DOFs. Larger meshes should rely on the
deck's explicit BCs and skip this path entirely.
Accepts either a sparse ``csr_array`` or a dense ndarray.
``avoid``: optional bool ``(n,)`` mask over local indices. ``True``
entries signal "do not pin this DOF if any other null-vector
component is even ε-non-zero" — used by the static path to keep
the auto-pin away from DOFs the user has loaded (forces or
prescribed displacements). Without this guard the pin can land
on a load-carrying DOF and silently zero the answer there.
"""
n = int(K_ff.shape[0])
if n == 0:
return []
if hasattr(K_ff, "todense"):
dense = np.asarray(K_ff.todense(), dtype=np.float64)
else:
dense = np.asarray(K_ff, dtype=np.float64)
sym = 0.5 * (dense + dense.T)
w, v = np.linalg.eigh(sym)
scale = float(np.abs(w).max()) if w.size else 1.0
tol_eig = 1e-10 * scale
null_cols = np.where(np.abs(w) < tol_eig)[0]
if null_cols.size == 0:
return []
pins: list[int] = []
pinned_set: set[int] = set()
avoid_arr = avoid if avoid is not None else np.zeros(n, dtype=bool)
for col in null_cols:
vec = v[:, col]
order = np.argsort(-np.abs(vec))
# Two-pass: first try to pin a non-loaded DOF; if all
# significant null entries are on loaded DOFs, fall back to
# the absolute-dominant entry (still better than no pin →
# singular solve).
chosen: int | None = None
for idx in order:
if int(idx) in pinned_set:
continue
if abs(float(vec[idx])) < 1e-12:
break
if not bool(avoid_arr[int(idx)]):
chosen = int(idx)
break
if chosen is None:
for idx in order:
if int(idx) in pinned_set:
continue
if abs(float(vec[idx])) < 1e-12:
break
chosen = int(idx)
break
if chosen is None:
continue
pins.append(chosen)
pinned_set.add(chosen)
return pins