"""Per-node stress and strain recovery for disk-backed results (T1-P5).
Bridges three pieces that live on different modules:
* :meth:`femorph_solver.model.Model.strain` — per-element, per-node
engineering-shear Voigt strain from a DOF-ordered displacement.
* :meth:`femorph_solver.result._material.MaterialTable.eval_c_isotropic`
— per-material ``(6, 6)`` elasticity matrix in Voigt order.
* :func:`femorph_solver.elements._stress.stress_from_strain` — the
``C · ε`` contraction.
The public entry points :func:`compute_nodal_strain` and
:func:`compute_nodal_stress` turn a displacement vector into an
``(n_points, 6)`` Voigt strain or stress array at the grid's nodes,
using the Model's element → material → node graph and averaging
where multiple elements share a node.
References
----------
* Nodal stress recovery by per-element evaluation at element nodes
followed by multi-element averaging (the canonical nodal-recovery
parlance): Cook, Malkus, Plesha, Witt, *Concepts and Applications
of Finite Element Analysis*, 4th ed., Wiley, 2002, §6.14–§6.15
(stress recovery from displacement; nodal averaging discussion).
* Alternative Zienkiewicz–Zhu superconvergent patch recovery (SPR)
— roadmap upgrade from arithmetic averaging to SPR on the same
patch data: Zienkiewicz, O.C. and Zhu, J.Z., "The superconvergent
patch recovery and a posteriori error estimates. Part 1: The
recovery technique," *IJNME* 33 (1992), pp. 1331–1364.
* ``σ = C·ε`` contraction mechanics (Voigt, engineering shear):
Cook et al. §3.7 (see also :mod:`femorph_solver.elements._stress`).
"""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
from femorph_solver.elements._stress import stress_from_strain
from femorph_solver.result._material import MaterialTable
if TYPE_CHECKING:
from femorph_solver.model import Model
__all__ = [
"compute_nodal_strain",
"compute_nodal_stress",
"nodal_to_dof_vector",
]
[docs]
def nodal_to_dof_vector(
nodal: np.ndarray,
model: Model,
) -> np.ndarray:
"""Gather a ``(n_points, 6)`` per-node array into a DOF-ordered vector.
Inverse of the scatter that each solver-result ``.save()`` applies
when projecting a DOF vector onto per-node per-component rows.
Used by the stress-recovery path to feed a DOF vector into
:meth:`Model.strain`.
Parameters
----------
nodal : numpy.ndarray
``(n_points, ≥ 3)`` per-node displacement. Extra columns
(e.g. rotations) are indexed by :meth:`Model.dof_map`.
model : femorph_solver.model.Model
Provides the DOF layout and the mapping from node
numbers back to grid point indices.
Returns
-------
numpy.ndarray
``(n_dof,)`` DOF-indexed float64 vector.
"""
nodal = np.asarray(nodal, dtype=np.float64)
if nodal.ndim != 2 or nodal.shape[1] < 3:
raise ValueError(f"expected (n_points, >= 3) nodal array, got {nodal.shape}")
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_rows = 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)
return np.ascontiguousarray(nodal[node_rows, comp])
def _average_per_node_field(
model: Model,
per_elem_field: dict[int, np.ndarray],
) -> np.ndarray:
"""Average a ``{elem_num: (n_nodes_in_elem, 6)}`` field onto grid nodes.
Walks every element in ``per_elem_field``, scatters each row onto
the matching grid point (looked up by ANSYS node number), and
returns the unweighted arithmetic mean across every element that
touches a node. Nodes with no incident elements are zero.
Parameters
----------
model : femorph_solver.model.Model
Source of the element → node graph (``model.grid``).
per_elem_field : dict[int, numpy.ndarray]
``{elem_num: (n_nodes_in_elem, 6)}`` Voigt field, indexed by
ANSYS element number. The row order along axis 0 must match
the element's local node order (the ``eel_batch`` convention).
Returns
-------
numpy.ndarray
``(n_points, 6)`` Voigt field on the grid points, in the
``ansys_node_num`` order.
Notes
-----
Averaging is unweighted by element volume — same convention as
MAPDL's PRNSOL. Both the stress and the strain recovery paths
share this primitive so they produce identical averaging.
"""
grid = model.grid
if not per_elem_field:
return np.zeros((grid.n_points, 6), dtype=np.float64)
# Lift the cell-data / connectivity arrays once so the per-element
# loop is a pure indexing path (no per-element ``Model.element_info``
# call, no per-element ``np.searchsorted``).
elem_nums_all = np.asarray(grid.cell_data["ansys_elem_num"], dtype=np.int64)
cell_conn = np.asarray(grid.cell_connectivity, dtype=np.int64)
cell_offset = np.asarray(grid.offset, dtype=np.int64)
pt_nums_by_vtk = np.asarray(grid.point_data["ansys_node_num"], dtype=np.int64)
elem_to_cell = {int(en): i for i, en in enumerate(elem_nums_all)}
node_to_row = {int(n): i for i, n in enumerate(pt_nums_by_vtk)}
accum = np.zeros((grid.n_points, 6), dtype=np.float64)
count = np.zeros(grid.n_points, dtype=np.int64)
for elem_num, field in per_elem_field.items():
cell_idx = elem_to_cell[int(elem_num)]
vtk_ids = cell_conn[cell_offset[cell_idx] : cell_offset[cell_idx + 1]]
nns = [int(pt_nums_by_vtk[j]) for j in vtk_ids]
rows = np.asarray(field, dtype=np.float64)
for local_i, n in enumerate(nns):
r = node_to_row[int(n)]
accum[r] += rows[local_i]
count[r] += 1
out = np.zeros_like(accum)
mask = count > 0
out[mask] = accum[mask] / count[mask, None]
return out
[docs]
def compute_nodal_strain(
model: Model,
displacement: np.ndarray,
) -> np.ndarray:
"""Compute ``(n_points, 6)`` Voigt strain at every mesh node.
Per-element strain comes from :meth:`Model.strain`, which evaluates
:math:`\\varepsilon = B(\\xi_\\text{node})\\,u_e` directly at each
element's own nodes. Where multiple elements touch a node the
strains are averaged with equal weight (the same convention as the
stress recovery in :func:`compute_nodal_stress` and as MAPDL's
PRNSOL).
Parameters
----------
model : femorph_solver.model.Model
Source of the strain kernel (``model.strain(u)``). Element
cells / nodes are read off ``model.grid``.
displacement : numpy.ndarray
``(N,)`` DOF-indexed displacement in :meth:`Model.dof_map`
ordering. Typically a static-solve displacement or a modal
mode shape.
Returns
-------
numpy.ndarray
``(n_points, 6)`` engineering-shear Voigt strain
``[ε_xx, ε_yy, ε_zz, γ_xy, γ_yz, γ_xz]``, averaged across
every element that touches a node. Nodes with no incident
elements are zero.
See Also
--------
compute_nodal_stress :
Same averaging pattern, but contracted with the per-material
elasticity matrix to produce Voigt stress.
femorph_solver.model.Model.strain :
Per-element source of the strain rows averaged here.
"""
eps_by_elem = model.strain(np.asarray(displacement, dtype=np.float64))
return _average_per_node_field(model, eps_by_elem)
[docs]
def compute_nodal_stress(
model: Model,
displacement: np.ndarray,
) -> np.ndarray:
"""Compute ``(n_points, 6)`` Voigt stress at every mesh node.
The caller is responsible for passing a Model with complete
materials (at least ``EX`` + ``PRXY`` on every referenced MAT id)
and element-type registrations matching the grid's
``ansys_elem_type_num`` values.
Parameters
----------
model : femorph_solver.model.Model
Source of both the strain kernel (``model.strain(u)``) and the
materials (``model.materials``). Element cells / nodes are
read off ``model.grid``.
displacement : numpy.ndarray
``(N,)`` DOF-indexed displacement in :meth:`Model.dof_map`
ordering. Typically a static-solve displacement or a modal
mode shape.
Returns
-------
numpy.ndarray
``(n_points, 6)`` Voigt stress, averaged across every element
that touches a node. Nodes with no incident elements are
zero.
See Also
--------
compute_nodal_strain :
Same averaging pattern on the underlying strain field.
"""
grid = model.grid
strain_by_elem = model.strain(np.asarray(displacement, dtype=np.float64))
if not strain_by_elem:
return np.zeros((grid.n_points, 6), dtype=np.float64)
# One C matrix per MAT id — cheap to precompute and reuse across
# every element of that material.
table = MaterialTable.from_model(model)
c_by_mat: dict[int, np.ndarray] = {}
for mat_id in table.ids:
c_by_mat[int(mat_id)] = table.eval_c_isotropic(int(mat_id))
mat_ids_all = np.asarray(grid.cell_data["ansys_material_type"], dtype=np.int64)
elem_nums_all = np.asarray(grid.cell_data["ansys_elem_num"], dtype=np.int64)
elem_to_cell = {int(en): i for i, en in enumerate(elem_nums_all)}
# Apply ``σ = C·ε`` per element, then hand the per-element stress
# dict to the same multi-element averaging primitive that
# :func:`compute_nodal_strain` uses — so both fields are guaranteed
# to use identical averaging weights.
stress_by_elem: dict[int, np.ndarray] = {}
for elem_num, eps in strain_by_elem.items():
cell_idx = elem_to_cell[int(elem_num)]
mat_id = int(mat_ids_all[cell_idx])
try:
c = c_by_mat[int(mat_id)]
except KeyError as exc: # pragma: no cover - defensive
raise KeyError(
f"element {elem_num} references MAT id {mat_id} which has no "
"C matrix; set EX/PRXY on that material and retry"
) from exc
stress_by_elem[int(elem_num)] = stress_from_strain(np.asarray(eps, dtype=np.float64), c)
return _average_per_node_field(model, stress_by_elem)