{
  "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# TET10 \u2014 prescribed uniform strain on a 10-node tet\n\nSingle TET10 (10-node quadratic tet) driven by a fully-prescribed\ndisplacement field that matches a uniform strain state. Under\nall-Dirichlet BCs the solve has a unique solution equal to the\nprescribed field, so the recovered strain tensor matches the analytical\none to machine precision. Useful as a self-contained smoke test that\n``TET10`` + :meth:`femorph_solver.result.StaticResult.elastic_strain` are in lockstep.\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_QUADRATIC_TETRA\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reference 10-node unit tet\nCorners at (0,0,0), (1,0,0), (0,1,0), (0,0,1); mid-edge nodes in\nVTK_QUADRATIC_TETRA order halfway along each edge.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.1e11  # Pa\nNU = 0.30\nEPS_XX = 5.0e-4\nEPS_YY = -NU * EPS_XX  # Poisson contraction\n\nI = np.array([0.0, 0.0, 0.0])\nJ = np.array([1.0, 0.0, 0.0])\nK = np.array([0.0, 1.0, 0.0])\nL = np.array([0.0, 0.0, 1.0])\ncoords = np.array(\n    [\n        I,\n        J,\n        K,\n        L,\n        0.5 * (I + J),  # 5\n        0.5 * (J + K),  # 6\n        0.5 * (K + I),  # 7\n        0.5 * (I + L),  # 8\n        0.5 * (J + L),  # 9\n        0.5 * (K + L),  # 10\n    ],\n    dtype=np.float64,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the model + prescribed displacement\nApply ``u = (\u03b5xx\u00b7x, \u03b5yy\u00b7y, \u03b5yy\u00b7z)`` at every node.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cells = np.concatenate([[10], np.arange(10, dtype=np.int64)])\ncell_types = np.array([VTK_QUADRATIC_TETRA], dtype=np.uint8)\ngrid = pv.UnstructuredGrid(cells, cell_types, coords)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.TET10, material={\"EX\": E, \"PRXY\": NU})\n\nfor i, (x, y, z) in enumerate(coords, start=1):\n    m.fix(nodes=[i], dof=\"UX\", value=EPS_XX * x)\n    m.fix(nodes=[i], dof=\"UY\", value=EPS_YY * y)\n    m.fix(nodes=[i], dof=\"UZ\", value=EPS_YY * z)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve + verify strain recovery\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\neps = res.elastic_strain(model=m)\n\nprint(f\"\u03b5xx prescribed = {EPS_XX:.3e}  |  recovered (mean) = {eps[:, 0].mean():.3e}\")\nprint(f\"\u03b5yy prescribed = {EPS_YY:.3e}  |  recovered (mean) = {eps[:, 1].mean():.3e}\")\nnp.testing.assert_allclose(eps[:, 0], EPS_XX, rtol=1e-10)\nnp.testing.assert_allclose(eps[:, 1], EPS_YY, rtol=1e-10)\nnp.testing.assert_allclose(eps[:, 2], EPS_YY, rtol=1e-10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deformed tet coloured by \u03b5xx\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid = femorph_solver.io.static_result_to_grid(m, res)\ngrid.point_data[\"eps_xx\"] = eps[:, 0]\n\nplotter = pv.Plotter(off_screen=True)\nplotter.add_mesh(grid, style=\"wireframe\", color=\"gray\", line_width=1, opacity=0.4)\nplotter.add_mesh(\n    grid.warp_by_vector(\"displacement\", factor=200.0),\n    scalars=\"eps_xx\",\n    show_edges=True,\n    cmap=\"viridis\",\n    scalar_bar_args={\"title\": \"\u03b5xx\"},\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
}