"""Linear transient solver (Newmark-β time integration).
Integrates ``M ü + C u̇ + K u = F(t)`` over a uniform timestep grid.
Dirichlet constraints fix both the displacement and (implicitly) the
velocity / acceleration at those DOFs to zero.
Newmark-β recurrence (standard form with parameters β, γ)::
u_{n+1} = u_n + dt · u̇_n + dt² · [(1/2 − β) ü_n + β ü_{n+1}]
u̇_{n+1} = u̇_n + dt · [(1 − γ) ü_n + γ ü_{n+1}]
Rearranging ``M ü_{n+1} + C u̇_{n+1} + K u_{n+1} = F_{n+1}`` in terms of
``u_{n+1}`` gives the effective system::
A u_{n+1} = b_{n+1}
A = a0·M + a1·C + K
b = F_{n+1} + M(a0 u_n + a2 u̇_n + a3 ü_n) + C(a1 u_n + a4 u̇_n + a5 ü_n)
with the usual coefficient shorthands ``a0 = 1/(β dt²)``,
``a1 = γ/(β dt)``, ``a2 = 1/(β dt)``, ``a3 = 1/(2β) − 1``,
``a4 = γ/β − 1``, ``a5 = dt/2 · (γ/β − 2)``.
Defaults ``(β, γ) = (1/4, 1/2)`` are the unconditionally-stable
average-acceleration method. This implementation is implicit and
requires ``β > 0`` — explicit (β=0) central-difference would need a
different solver path that factors the mass matrix directly.
References
----------
- Newmark, N. M. "A method of computation for structural dynamics."
J. Eng. Mech. Div. (ASCE) 85 (EM3), 67-94 (1959). [original
Newmark-β integration scheme]
- Hughes, T. J. R. *The Finite Element Method: Linear Static and
Dynamic Finite Element Analysis*, Dover (2000), §9.1-9.2.
[stability + accuracy analysis of the (β, γ) family; conditions
for second-order accuracy and unconditional stability]
- Bathe, K.-J. *Finite Element Procedures*, 2nd ed., §9.2.
[implementation form of the effective-stiffness recurrence
used above and the reduction to a per-step linear solve]
"""
from __future__ import annotations
from collections.abc import Callable
from dataclasses import dataclass
from pathlib import Path
from typing import TYPE_CHECKING
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
from .._threads import limit as _thread_limit
if TYPE_CHECKING:
from femorph_solver.model import Model
[docs]
@dataclass
class TransientResult:
"""Result of a linear transient solve, sampled on the time grid.
Attributes
----------
time : numpy.ndarray
``(n_steps + 1,)`` sample times [s].
displacement : numpy.ndarray
``(n_steps + 1, N)`` DOF-indexed displacement; row per step.
velocity : numpy.ndarray
``(n_steps + 1, N)`` DOF-indexed velocity.
acceleration : numpy.ndarray
``(n_steps + 1, N)`` DOF-indexed acceleration.
"""
time: np.ndarray
displacement: np.ndarray
velocity: np.ndarray
acceleration: np.ndarray
def __repr__(self) -> str:
n_steps = int(self.time.shape[0])
n_dof = int(self.displacement.shape[1]) if self.displacement.ndim == 2 else 0
if n_steps:
t0 = float(self.time[0])
t1 = float(self.time[-1])
return f"TransientResult(n_steps={n_steps}, N={n_dof}, t∈[{t0:.4g}, {t1:.4g}] s)"
return f"TransientResult(n_steps=0, N={n_dof})"
[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 transient result to a disk-backed ``.pv`` file.
Re-projects the DOF-indexed :attr:`displacement` /
:attr:`velocity` / :attr:`acceleration` onto per-step, per-node
``(n_steps, n_points, n_dof_per_node)`` arrays via
:meth:`~femorph_solver.model.Model.dof_map`, then hands off to
:func:`~femorph_solver.result.transient.write_transient_result`.
Parameters
----------
path : str or pathlib.Path
Destination file.
model : femorph_solver.model.Model
Model whose K / M matrices this result was integrated on;
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.transient import write_transient_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)
n_steps = int(self.displacement.shape[0])
def _project(dof_matrix: np.ndarray) -> np.ndarray:
"""``(n_steps, N)`` → ``(n_steps, n_points, 6)`` scatter."""
out = np.zeros((n_steps, grid.n_points, 6), dtype=np.float64)
for k in range(n_steps):
out[k, node_row, comp] = dof_matrix[k]
return out
material_table = MaterialTable.from_model(model, unit_system=str(unit_system))
return write_transient_result(
path,
grid,
self.time,
_project(self.displacement),
velocity=_project(self.velocity),
acceleration=_project(self.acceleration),
material_table=material_table,
unit_system=unit_system,
extra_field_data=extra_field_data,
)
[docs]
def solve_transient(
K: sp.csr_array,
M: sp.csr_array,
F: np.ndarray | Callable[[float], np.ndarray],
*,
dt: float,
n_steps: int,
prescribed: dict[int, float] | None = None,
C: sp.csr_array | None = None,
u0: np.ndarray | None = None,
v0: np.ndarray | None = None,
beta: float = 0.25,
gamma: float = 0.5,
thread_limit: int | None = None,
) -> TransientResult:
"""Integrate ``M ü + C u̇ + K u = F(t)`` with Newmark-β.
Parameters
----------
K, M
Global stiffness and mass matrices (``(N, N)`` sparse).
F
Load vector. Either a fixed ``(N,)`` array (step-constant load) or
a callable ``t -> (N,)``.
dt
Uniform timestep.
n_steps
Number of steps to take. Returns ``n_steps + 1`` samples including
the initial state at ``t = 0``.
prescribed
Dirichlet constraints ``{global_dof_index: prescribed_value}``.
Constrained DOFs are held at ``u_c`` with zero velocity and
acceleration; their rows/cols are dropped from ``A``.
C
Damping matrix. ``None`` means no damping.
u0, v0
Initial displacement and velocity. Default to zero.
beta, gamma
Newmark parameters. Defaults ``(1/4, 1/2)`` give the
unconditionally-stable average-acceleration method. ``(0, 1/2)``
reduces to central-difference (explicit, conditionally stable).
thread_limit : int or None, optional
Cap on BLAS / OpenMP threads used inside the factorisation and
per-step back-solve. ``None`` defers to the global limit (see
:func:`femorph_solver.set_thread_limit`).
"""
if dt <= 0:
raise ValueError(f"dt must be positive, got {dt}")
if n_steps <= 0:
raise ValueError(f"n_steps must be positive, got {n_steps}")
if beta <= 0:
raise ValueError(f"beta must be positive for this implicit integrator, got {beta}")
N = K.shape[0]
prescribed = prescribed or {}
mask_fixed = np.zeros(N, dtype=bool)
u_fixed = np.zeros(N, dtype=np.float64)
for idx, val in prescribed.items():
mask_fixed[int(idx)] = True
u_fixed[int(idx)] = float(val)
free = ~mask_fixed
K_csc = K.tocsc()
M_csc = M.tocsc()
C_csc = C.tocsc() if C is not None else None
K_ff = K_csc[free][:, free]
M_ff = M_csc[free][:, free]
C_ff = C_csc[free][:, free] if C_csc is not None else None
# Newmark coefficients.
a0 = 1.0 / (beta * dt * dt)
a1 = gamma / (beta * dt)
a2 = 1.0 / (beta * dt)
a3 = 1.0 / (2.0 * beta) - 1.0
a4 = gamma / beta - 1.0
a5 = dt / 2.0 * (gamma / beta - 2.0)
A_ff = a0 * M_ff + K_ff
if C_ff is not None:
A_ff = A_ff + a1 * C_ff
A_ff = A_ff.tocsc()
time = np.arange(n_steps + 1, dtype=np.float64) * dt
U = np.zeros((n_steps + 1, N), dtype=np.float64)
V = np.zeros((n_steps + 1, N), dtype=np.float64)
A_acc = np.zeros((n_steps + 1, N), dtype=np.float64)
u = np.zeros(N, dtype=np.float64) if u0 is None else np.asarray(u0, dtype=np.float64).copy()
v = np.zeros(N, dtype=np.float64) if v0 is None else np.asarray(v0, dtype=np.float64).copy()
u[mask_fixed] = u_fixed[mask_fixed]
v[mask_fixed] = 0.0
F_arr = F if callable(F) else np.asarray(F, dtype=np.float64)
F0 = F_arr(0.0) if callable(F) else F_arr
with _thread_limit(thread_limit):
A_factor = spla.splu(A_ff)
# Initial acceleration: M a0 = F0 − K u0 − C v0, on free DOFs only.
rhs0 = F0[free] - K_ff @ u[free]
if C_ff is not None:
rhs0 = rhs0 - C_ff @ v[free]
a = np.zeros(N, dtype=np.float64)
a[free] = spla.spsolve(M_ff.tocsc(), rhs0)
U[0, :] = u
V[0, :] = v
A_acc[0, :] = a
for step in range(n_steps):
t_next = time[step + 1]
F_next = F_arr(t_next) if callable(F) else F_arr
# Effective right-hand side on free DOFs.
rhs = F_next[free].copy()
rhs = rhs + M_ff @ (a0 * u[free] + a2 * v[free] + a3 * a[free])
if C_ff is not None:
rhs = rhs + C_ff @ (a1 * u[free] + a4 * v[free] + a5 * a[free])
u_next_free = A_factor.solve(rhs)
a_next = np.zeros(N, dtype=np.float64)
v_next = np.zeros(N, dtype=np.float64)
a_next[free] = a0 * (u_next_free - u[free]) - a2 * v[free] - a3 * a[free]
v_next[free] = v[free] + dt * ((1.0 - gamma) * a[free] + gamma * a_next[free])
u = np.zeros(N, dtype=np.float64)
u[free] = u_next_free
u[mask_fixed] = u_fixed[mask_fixed]
v = v_next
a = a_next
U[step + 1, :] = u
V[step + 1, :] = v
A_acc[step + 1, :] = a
return TransientResult(time=time, displacement=U, velocity=V, acceleration=A_acc)