SOLID187 — prescribed uniform strain on a 10-node tet

Single SOLID187 (10-node quadratic tet) driven by a fully-prescribed displacement field that matches a uniform strain state. Under all-Dirichlet BCs the solve has a unique solution equal to the prescribed field, so the recovered strain tensor matches the analytical one to machine precision. Useful as a self-contained smoke test that SOLID187 + femorph_solver.Model.eel() are in lockstep.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Reference 10-node unit tet

Corners at (0,0,0), (1,0,0), (0,1,0), (0,0,1); mid-edge nodes in MAPDL / VTK order halfway along each edge.

E = 2.1e11  # Pa
NU = 0.30
EPS_XX = 5.0e-4
EPS_YY = -NU * EPS_XX  # Poisson contraction

I = np.array([0.0, 0.0, 0.0])
J = np.array([1.0, 0.0, 0.0])
K = np.array([0.0, 1.0, 0.0])
L = np.array([0.0, 0.0, 1.0])
coords = np.array(
    [
        I,
        J,
        K,
        L,
        0.5 * (I + J),  # 5
        0.5 * (J + K),  # 6
        0.5 * (K + I),  # 7
        0.5 * (I + L),  # 8
        0.5 * (J + L),  # 9
        0.5 * (K + L),  # 10
    ],
    dtype=np.float64,
)

Build the model + prescribed displacement

Apply u = (εxx·x, εyy·y, εyy·z) at every node.

m = femorph_solver.Model()
m.et(1, "SOLID187")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
for i, (x, y, z) in enumerate(coords, start=1):
    m.n(i, x, y, z)
m.e(*range(1, 11))

for i, (x, y, z) in enumerate(coords, start=1):
    m.d(i, "UX", EPS_XX * x)
    m.d(i, "UY", EPS_YY * y)
    m.d(i, "UZ", EPS_YY * z)

/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:57: DeprecationWarning: Model.et(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).et(et_id, name)` for line-by-line APDL deck porting, or the native `Model.assign("HEX8", material)` for new code.
  m.et(1, "SOLID187")
/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:58: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("EX", 1, E)
/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:59: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("PRXY", 1, NU)
/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:61: DeprecationWarning: Model.n(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).n(num, x, y, z)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.n(i, x, y, z)
/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:62: DeprecationWarning: Model.e(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).e(*node_nums)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.e(*range(1, 11))
/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:65: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(i, "UX", EPS_XX * x)
/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:66: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(i, "UY", EPS_YY * y)
/home/runner/_work/solver/solver/examples/elements/solid187/example_solid187.py:67: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(i, "UZ", EPS_YY * z)

Solve + verify strain recovery

res = m.solve()
eps = m.eel(res.displacement)

print(f"εxx prescribed = {EPS_XX:.3e}  |  recovered (mean) = {eps[:, 0].mean():.3e}")
print(f"εyy prescribed = {EPS_YY:.3e}  |  recovered (mean) = {eps[:, 1].mean():.3e}")
np.testing.assert_allclose(eps[:, 0], EPS_XX, rtol=1e-10)
np.testing.assert_allclose(eps[:, 1], EPS_YY, rtol=1e-10)
np.testing.assert_allclose(eps[:, 2], EPS_YY, rtol=1e-10)

εxx prescribed = 5.000e-04  |  recovered (mean) = 5.000e-04
εyy prescribed = -1.500e-04  |  recovered (mean) = -1.500e-04

Render the deformed tet coloured by εxx

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)}
grid.point_data["eps_xx"] = np.array(
    [eps[node_to_idx[int(nn)], 0] for nn in grid.point_data["ansys_node_num"]]
)

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(grid, style="wireframe", color="gray", line_width=1, opacity=0.4)
plotter.add_mesh(
    grid.warp_by_vector("displacement", factor=200.0),
    scalars="eps_xx",
    show_edges=True,
    cmap="viridis",
    scalar_bar_args={"title": "εxx"},
)
plotter.add_axes()
plotter.show()