SOLID185 — elastic-strain post-processing

Solve a SOLID185 flat plate under uniaxial tension and recover the full 6-component elastic-strain tensor on the mesh with femorph_solver.Model.eel() — the MAPDL S,EPEL equivalent.

Model.eel(u) returns the nodal-averaged Voigt strain (n_nodes, 6) (PLNSOL-style, default) or the per-element dict {elem_num: (n_nodes_in_elem, 6)} (PLESOL-style) when called with nodal_avg=False. Strain is computed at each element’s own nodes as \(\varepsilon(\xi_\text{node}) = B(\xi_\text{node})\cdot u_e\) — no RST round-trip, no disk write.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Problem setup

A 1 m × 0.4 m × 0.05 m steel plate meshed as a 20 × 8 × 1 SOLID185 brick (160 elements). The x = 0 face is held in UX (symmetry), a single pin at the origin kills the UY / UZ rigid-body modes, and the x = LX face is pulled by a total force F split over its corner nodes.

E = 2.1e11  # Pa
NU = 0.30
RHO = 7850.0
LX, LY, LZ = 1.0, 0.4, 0.05
NX, NY, NZ = 20, 8, 1
F_TOTAL = 1.0e5  # N

xs = np.linspace(0.0, LX, NX + 1)
ys = np.linspace(0.0, LY, NY + 1)
zs = np.linspace(0.0, LZ, NZ + 1)
xx, yy, zz = np.meshgrid(xs, ys, zs, indexing="ij")
points = np.stack([xx.ravel(), yy.ravel(), zz.ravel()], axis=1)


# Hex connectivity in MAPDL / VTK order.
def _node_id(i: int, j: int, k: int) -> int:
    return (i * (NY + 1) + j) * (NZ + 1) + k + 1  # 1-based


cells = []
for i in range(NX):
    for j in range(NY):
        for k in range(NZ):
            cells.append(
                [
                    _node_id(i, j, k),
                    _node_id(i + 1, j, k),
                    _node_id(i + 1, j + 1, k),
                    _node_id(i, j + 1, k),
                    _node_id(i, j, k + 1),
                    _node_id(i + 1, j, k + 1),
                    _node_id(i + 1, j + 1, k + 1),
                    _node_id(i, j + 1, k + 1),
                ]
            )

Build the femorph-solver model

m = femorph_solver.Model()
m.et(1, "SOLID185")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
m.mp("DENS", 1, RHO)
for nn, (x, y, z) in enumerate(points, start=1):
    m.n(nn, x, y, z)
for conn in cells:
    m.e(*conn)

# Symmetry BC: x=0 face clamped in UX; single pin at the origin in UY/UZ.
x0_nodes = [nn for nn, p in enumerate(points, start=1) if p[0] < 1e-9]
for nn in x0_nodes:
    m.d(nn, "UX")
m.d(1, "UY")
m.d(1, "UZ")

# Traction on x=LX face: split F_TOTAL over its nodes.
x_end_nodes = [nn for nn, p in enumerate(points, start=1) if p[0] > LX - 1e-9]
for nn in x_end_nodes:
    m.f(nn, "FX", F_TOTAL / len(x_end_nodes))

Static solve

res = m.solve()

Recover elastic strain

Default call returns nodal-averaged strain of shape (n_nodes, 6): columns are [εxx, εyy, εzz, γxy, γyz, γxz] with engineering shears (matching MAPDL’s S,EPEL output).

eps = m.eel(res.displacement)
print(f"eps shape: {eps.shape}")

# Analytical: uniform σxx = F_TOTAL / (LY · LZ), εxx = σ / E,
# εyy = εzz = -ν · εxx.
sigma_xx = F_TOTAL / (LY * LZ)
eps_xx_expected = sigma_xx / E
eps_yy_expected = -NU * eps_xx_expected
print(f"εxx expected = {eps_xx_expected:.3e}")
print(f"εxx recovered (mean over nodes) = {eps[:, 0].mean():.3e}")
print(f"εyy recovered (mean)            = {eps[:, 1].mean():.3e}")
print(f"εyy analytical                  = {eps_yy_expected:.3e}")

eps shape: (378, 6)
εxx expected = 2.381e-05
εxx recovered (mean over nodes) = 2.416e-05
εyy recovered (mean)            = -7.365e-06
εyy analytical                  = -7.143e-06

nodal_avg=False returns per-element arrays keyed by element number — the PLESOL equivalent. Useful when you want to see jumps at element boundaries or compute element-wise strain norms.

per_elem = m.eel(res.displacement, nodal_avg=False)
first_elem = next(iter(per_elem))
print(
    f"per-element dict has {len(per_elem)} elements; "
    f"first key = {first_elem}, "
    f"strain block shape = {per_elem[first_elem].shape}"
)

per-element dict has 160 elements; first key = 1, strain block shape = (8, 6)

Visualise εxx on the deformed mesh

femorph_solver.io.static_result_to_grid() scatters the DOF-indexed displacement vector onto (n_points, 3) UX/UY/UZ point data in one call — no hand-rolled dof-map loop required. We then paint εxx onto the same grid by mapping the eel output (indexed by femorph_solver.Model.node_numbers) onto ansys_node_num.

grid = femorph_solver.io.static_result_to_grid(m, res)
node_nums = m.node_numbers
node_to_idx = {int(nn): i for i, nn in enumerate(node_nums)}
point_eps_xx = np.array([eps[node_to_idx[int(nn)], 0] for nn in grid.point_data["ansys_node_num"]])
grid.point_data["eps_xx"] = point_eps_xx

warped = grid.warp_by_vector("displacement", factor=200.0)

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(
    grid,
    style="wireframe",
    color="gray",
    line_width=1,
    opacity=0.4,
    label="undeformed",
)
plotter.add_mesh(
    warped,
    scalars="eps_xx",
    show_edges=True,
    cmap="viridis",
    scalar_bar_args={"title": "εxx"},
    label="εxx (deformed ×200)",
)
plotter.add_legend()
plotter.add_axes()
plotter.view_xy()
plotter.show()