SOLID185 — uniaxial tension of a unit cube

Single SOLID185 element (a 1 m unit cube) with one face clamped in the pulling direction and the opposite face loaded. The displacement field is uniform — u_x = σ L / E on the loaded face — which gives a clean check on the full 24 × 24 stiffness plus the boundary/load handling.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Problem data

Steel. Unit cube ⇒ cross-section area A = 1 and length L = 1, so applying a total force F produces nominal stress σ = F.

E = 2.1e11  # Pa
NU = 0.3
RHO = 7850.0
F_TOTAL = 1.0e6  # N distributed over the +x face (4 corner nodes)

Build the model

Eight corners in MAPDL / VTK hex order. Clamp the -x face (nodes 1, 4, 5, 8) in x and pin one node’s y/z to prevent rigid-body rotation. Split the total +x force onto the four +x-face nodes.

coords = np.array(
    [
        [0.0, 0.0, 0.0],  # 1
        [1.0, 0.0, 0.0],  # 2
        [1.0, 1.0, 0.0],  # 3
        [0.0, 1.0, 0.0],  # 4
        [0.0, 0.0, 1.0],  # 5
        [1.0, 0.0, 1.0],  # 6
        [1.0, 1.0, 1.0],  # 7
        [0.0, 1.0, 1.0],  # 8
    ],
    dtype=np.float64,
)

m = femorph_solver.Model()
m.et(1, "SOLID185")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
m.mp("DENS", 1, RHO)
for i, (x, y, z) in enumerate(coords, start=1):
    m.n(i, x, y, z)
m.e(1, 2, 3, 4, 5, 6, 7, 8)

# Clamp x=0 face (nodes 1, 4, 5, 8) in UX
for nn in (1, 4, 5, 8):
    m.d(nn, "UX")
# Roller supports to kill rigid-body y/z motion without restraining Poisson
for nn in (1, 2, 5, 6):
    m.d(nn, "UY")
for nn in (1, 2, 3, 4):
    m.d(nn, "UZ")

# Apply F_TOTAL on x=1 face, split evenly over its 4 corners
F_each = F_TOTAL / 4.0
for nn in (2, 3, 6, 7):
    m.f(nn, "FX", F_each)

/home/runner/_work/solver/solver/examples/elements/solid185/example_solid185.py:51: 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, "SOLID185")
/home/runner/_work/solver/solver/examples/elements/solid185/example_solid185.py:52: 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/solid185/example_solid185.py:53: 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/solid185/example_solid185.py:54: 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("DENS", 1, RHO)
/home/runner/_work/solver/solver/examples/elements/solid185/example_solid185.py:56: 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/solid185/example_solid185.py:57: 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(1, 2, 3, 4, 5, 6, 7, 8)
/home/runner/_work/solver/solver/examples/elements/solid185/example_solid185.py:61: 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(nn, "UX")
/home/runner/_work/solver/solver/examples/elements/solid185/example_solid185.py:64: 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(nn, "UY")
/home/runner/_work/solver/solver/examples/elements/solid185/example_solid185.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(nn, "UZ")
/home/runner/_work/solver/solver/examples/elements/solid185/example_solid185.py:71: DeprecationWarning: Model.f(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).f(node, label, value)` for line-by-line APDL deck porting, or the native `Model.apply_force(node, fx=..., fy=..., fz=...)` for new code.
  m.f(nn, "FX", F_each)

Static solve + analytical comparison

Expected: u_x = σ L / E with σ = F_TOTAL / A = F_TOTAL and L = 1 → u_x = F_TOTAL / E on the +x face. Poisson contraction is u_y = u_z = −ν u_x but only the free nodes on those faces actually display it (here only UY at the +y face, UZ at the +z face).

res = m.solve()
dof = m.dof_map()

ux_expected = F_TOTAL / E
uy_expected = -NU * ux_expected

print(f"Expected u_x (+x face) = {ux_expected:.6e} m")
for nn in (2, 3, 6, 7):
    row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 0))[0][0]
    print(f"  node {nn} UX           = {res.displacement[row]:.6e}")
    assert np.isclose(res.displacement[row], ux_expected, rtol=1e-8)

print(f"Expected u_y (+y face) = {uy_expected:.6e} m")
for nn in (3, 4, 7, 8):
    row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 1))[0][0]
    val = res.displacement[row]
    print(f"  node {nn} UY           = {val:.6e}")
    assert np.isclose(val, uy_expected, rtol=1e-6)

Expected u_x (+x face) = 4.761905e-06 m
  node 2 UX           = 4.761905e-06
  node 3 UX           = 4.761905e-06
  node 6 UX           = 4.761905e-06
  node 7 UX           = 4.761905e-06
Expected u_y (+y face) = -1.428571e-06 m
  node 3 UY           = -1.428571e-06
  node 4 UY           = -1.428571e-06
  node 7 UY           = -1.428571e-06
  node 8 UY           = -1.428571e-06

Plot the deformed cube, coloured by displacement magnitude

grid = m.grid.copy()
disp = np.zeros((grid.n_points, 3), dtype=np.float64)
for i, nn in enumerate(grid.point_data["ansys_node_num"]):
    rows = np.where(dof[:, 0] == int(nn))[0]
    for r in rows:
        disp[i, int(dof[r, 1])] = res.displacement[r]
grid.point_data["displacement"] = disp

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(grid, style="wireframe", color="gray", line_width=2)
plotter.add_mesh(
    grid.warp_by_vector("displacement", factor=5.0e4),
    scalars=np.linalg.norm(disp, axis=1),
    show_edges=True,
    scalar_bar_args={"title": "disp [m]"},
)
plotter.add_axes()
plotter.show()