{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# HEX8 \u2014 uniaxial tension of a unit cube\n\nSingle HEX8 element (a 1 m unit cube) with one face clamped in the\npulling direction and the opposite face loaded. The displacement field\nis uniform \u2014 ``u_x = \u03c3 L / E`` on the loaded face \u2014 which gives a clean\ncheck on the full 24 \u00d7 24 stiffness plus the boundary/load handling.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport numpy as np\nimport pyvista as pv\nfrom vtkmodules.util.vtkConstants import VTK_HEXAHEDRON\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problem data\nSteel. Unit cube \u21d2 cross-section area ``A = 1`` and length ``L = 1``,\nso applying a total force ``F`` produces nominal stress ``\u03c3 = F``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.1e11  # Pa\nNU = 0.3\nRHO = 7850.0\nF_TOTAL = 1.0e6  # N distributed over the +x face (4 corner nodes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the model\nEight corners in VTK_HEXAHEDRON order. Clamp the -x face (nodes 1,\n4, 5, 8) in x and pin one node's y/z to prevent rigid-body rotation.\nSplit the total +x force onto the four +x-face nodes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "coords = np.array(\n    [\n        [0.0, 0.0, 0.0],  # 1\n        [1.0, 0.0, 0.0],  # 2\n        [1.0, 1.0, 0.0],  # 3\n        [0.0, 1.0, 0.0],  # 4\n        [0.0, 0.0, 1.0],  # 5\n        [1.0, 0.0, 1.0],  # 6\n        [1.0, 1.0, 1.0],  # 7\n        [0.0, 1.0, 1.0],  # 8\n    ],\n    dtype=np.float64,\n)\n\ncells = np.array([8, 0, 1, 2, 3, 4, 5, 6, 7], dtype=np.int64)\ncell_types = np.array([VTK_HEXAHEDRON], dtype=np.uint8)\ngrid = pv.UnstructuredGrid(cells, cell_types, coords)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.HEX8, material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})\n\n# Clamp x=0 face (nodes 1, 4, 5, 8) in UX\nfor nn in (1, 4, 5, 8):\n    m.fix(nodes=[nn], dof=\"UX\")\n# Roller supports to kill rigid-body y/z motion without restraining Poisson\nfor nn in (1, 2, 5, 6):\n    m.fix(nodes=[nn], dof=\"UY\")\nfor nn in (1, 2, 3, 4):\n    m.fix(nodes=[nn], dof=\"UZ\")\n\n# Apply F_TOTAL on x=1 face, split evenly over its 4 corners\nF_each = F_TOTAL / 4.0\nfor nn in (2, 3, 6, 7):\n    m.apply_force(nn, fx=F_each)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve + analytical comparison\nExpected: ``u_x = \u03c3 L / E`` with ``\u03c3 = F_TOTAL / A = F_TOTAL`` and\n``L = 1`` \u2192 ``u_x = F_TOTAL / E`` on the +x face. Poisson contraction\nis ``u_y = u_z = \u2212\u03bd u_x`` but only the free nodes on those faces\nactually display it (here only UY at the +y face, UZ at the +z face).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\ndof = m.dof_map()\n\nux_expected = F_TOTAL / E\nuy_expected = -NU * ux_expected\n\nprint(f\"Expected u_x (+x face) = {ux_expected:.6e} m\")\nfor nn in (2, 3, 6, 7):\n    row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 0))[0][0]\n    print(f\"  node {nn} UX           = {res.displacement[row]:.6e}\")\n    assert np.isclose(res.displacement[row], ux_expected, rtol=1e-8)\n\nprint(f\"Expected u_y (+y face) = {uy_expected:.6e} m\")\nfor nn in (3, 4, 7, 8):\n    row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 1))[0][0]\n    val = res.displacement[row]\n    print(f\"  node {nn} UY           = {val:.6e}\")\n    assert np.isclose(val, uy_expected, rtol=1e-6)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the deformed cube, coloured by displacement magnitude\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid = m.grid.copy()\ndisp = np.zeros((grid.n_points, 3), dtype=np.float64)\nfor i, nn in enumerate(grid.point_data[\"ansys_node_num\"]):\n    rows = np.where(dof[:, 0] == int(nn))[0]\n    for r in rows:\n        disp[i, int(dof[r, 1])] = res.displacement[r]\ngrid.point_data[\"displacement\"] = disp\n\nplotter = pv.Plotter(off_screen=True)\nplotter.add_mesh(grid, style=\"wireframe\", color=\"gray\", line_width=2)\nplotter.add_mesh(\n    grid.warp_by_vector(\"displacement\", factor=5.0e4),\n    scalars=np.linalg.norm(disp, axis=1),\n    show_edges=True,\n    scalar_bar_args={\"title\": \"disp [m]\"},\n)\nplotter.add_axes()\nplotter.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}