"""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
:data:`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).
"""
from __future__ import annotations
from pathlib import Path
import numpy as np
from femorph_solver.mapdl_api.io.binary import _open_binary, read_record
from femorph_solver.mapdl_api.io.rst_schema import parse_header
emat_header_keys: list[str] = [
"fun02",
"nume",
"numdof",
"lenu",
"lenbac",
"maxn",
"nlgeEMA",
"sstEMAT",
"nodref",
"lumpm",
"kygst",
"kygm",
"kycd",
"kygss",
"kygaf",
"kygrf",
"0",
"Glblenbac",
"ptrGBkl",
"ptrGBkh",
"ptrElmh",
"ptrFSTh",
"ptrLSTh",
"ptrBITh",
"ptrEHDh",
"ptrIDXh",
"numCE",
"maxLeng",
"ptrCEl",
"ptrCEh",
"ptrDOF",
"ptrBAC",
"ptrElml",
"ptrFSTl",
"ptrLSTl",
"ptrBITl",
"ptrEHDl",
"ptrIDXl",
"ptrendH",
"ptrendL",
]
element_matrix_header_keys: list[str] = [
"stkey",
"mkey",
"dkey",
"sskey",
"akey",
"nrkey",
"ikey",
"_",
"_",
"nmrow",
]
def _unpack_symmetric_lower(values: np.ndarray, n: int) -> np.ndarray:
"""Expand a packed lower-triangular array into an ``(n, n)`` dense matrix.
MAPDL stores symmetric element matrices in row-major lower-triangular
order: row 0 → col 0, row 1 → cols 0..1, row 2 → cols 0..2, ….
"""
expected = n * (n + 1) // 2
if values.size != expected:
raise ValueError(f"lower-tri packed length {values.size} != expected {expected} for n={n}")
out = np.zeros((n, n), dtype=values.dtype)
idx = 0
for row in range(n):
count = row + 1
out[row, :count] = values[idx : idx + count]
idx += count
# mirror lower triangle to upper to form a symmetric full matrix
out = out + np.tril(out, -1).T
return out
[docs]
def read_emat_element(
path: str | Path,
index: int,
) -> tuple[np.ndarray, dict[str, np.ndarray]]:
"""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, matrices
``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.
"""
header = read_emat_header(path)
nume = header["nume"]
if not (0 <= index < nume):
raise IndexError(f"element index {index} out of range (nume={nume})")
fp = _open_binary(path)
try:
idx_table, _ = read_record(fp, offset=header["ptrIDX"] * 4, dtype="<i4")
# Index table is (lo, hi) pairs; nume entries.
lo = int(idx_table[index])
hi = int(idx_table[index + nume]) if idx_table.size >= 2 * nume else 0
elem_word = lo | (hi << 32)
em_header_arr, _ = read_record(fp, offset=elem_word * 4, dtype="<i4")
em = parse_header(em_header_arr, element_matrix_header_keys)
nmrow = abs(em["nmrow"])
# nmrow<0 ⇒ lower-triangular packing; nmrow>0 is treated identically
# because real EMAT files from modern MAPDL always pack symmetric
# lower-tri for LINK180/BEAM188/etc.
unpack = _unpack_symmetric_lower
# DOF index table — 1-based global DOF numbers.
dof_arr, _ = read_record(fp, dtype="<i4")
dof_idx = dof_arr.astype(np.int64) - 1
matrices: dict[str, np.ndarray] = {}
def _read_sym():
vals, info = read_record(fp, dtype="<f8")
if vals.size == nmrow * nmrow:
return vals.reshape(nmrow, nmrow)
return unpack(vals, nmrow)
if em["stkey"]:
matrices["k"] = _read_sym()
if em["mkey"]:
matrices["m"] = _read_sym()
if em["dkey"]:
matrices["d"] = _read_sym()
if em["sskey"]:
matrices["ss"] = _read_sym()
if em["akey"]:
vec, _ = read_record(fp, dtype="<f8")
matrices["applied_force"] = vec
if em["nrkey"]:
vec, _ = read_record(fp, dtype="<f8")
matrices["newton_raphson"] = vec
if em["ikey"]:
vec, _ = read_record(fp, dtype="<f8")
matrices["imaginary_load"] = vec
finally:
fp.close()
return dof_idx, matrices
[docs]
def read_emat_element_matrices(
path: str | Path,
) -> list[tuple[np.ndarray, dict[str, np.ndarray]]]:
"""Read every element's (dof_idx, matrices) from the EMAT file."""
header = read_emat_header(path)
return [read_emat_element(path, i) for i in range(header["nume"])]