"""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
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 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,
)
[docs]
def solve_static(
K: sp.csr_array,
F: np.ndarray,
prescribed: dict[int, float],
*,
linear_solver: str = "auto",
thread_limit: int | None = None,
) -> 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}``
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 = K.shape[0]
u = np.zeros(N, dtype=np.float64)
mask_fixed = np.zeros(N, 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])
lu = Solver(K_ff, assume_spd=True)
u[free] = lu.solve(rhs)
R = np.zeros(N, 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)