{
  "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# Simply-supported plate under uniform pressure (Navier series)\n\nClassical thin-plate verification: a simply-supported square\nplate of side $a$ and thickness $h$ carries a\nuniform transverse pressure $q$.  The Navier double-sine\nexpansion (Timoshenko & Woinowsky-Krieger 1959 \u00a730 eqs. 115\u2013117)\ngives the centre deflection\n\n\\begin{align}w(x, y) = \\frac{16\\, q}{\\pi^{6} D}\n    \\sum_{m \\,\\text{odd}} \\sum_{n \\,\\text{odd}}\n    \\frac{\\sin(m \\pi x / a)\\, \\sin(n \\pi y / a)}\n         {m\\, n\\, \\bigl(m^{2}/a^{2} + n^{2}/a^{2}\\bigr)^{2}},\n    \\qquad\n    D = \\frac{E\\, h^{3}}{12\\, (1 - \\nu^{2})}.\\end{align}\n\nEvaluated at $(a/2, a/2)$ with summation truncated at the\nfirst ~25 odd-modes per direction the series converges to\n\n\\begin{align}w_\\mathrm{max} \\;\\approx\\; 0.00406\\, \\frac{q\\, a^{4}}{D}\n    \\qquad \\text{(Timoshenko \u00a731 Table 8).}\\end{align}\n\n## Implementation\n\nA $30 \\times 30 \\times 2$-element HEX8 mesh on the unit\nsquare plate, with $L/h = 50$ (default geometry).\n:class:`~femorph_solver.elements.hex8.Hex8` is invoked in\n**enhanced-strain** mode (Simo\u2013Rifai 1990) so the bending\nresponse of the thin-plate slab is recovered without\nshear-locking \u2014 the locked plain-Gauss form lands at ~50 %\nerror on this mesh; EAS closes that gap to ~5 %.\n\nSymmetry / SS BCs:\n\n* $u_z = 0$ along all four edges, full thickness \u2014 the\n  standard solid-element analogue of a Kirchhoff simple\n  support.\n* In-plane rigid-body translation pinned at one corner;\n  in-plane $u_y$ pinned at the opposite-x corner to kill\n  the residual rotation mode.\n\nPressure is lumped onto the top-face nodes by uniform area\nweighting (the lumping error is $O(1/n_x^{2})$, well below\nthe plate-theory error on this mesh).\n\n## References\n\n* Timoshenko, S. P. and Woinowsky-Krieger, S. (1959) *Theory\n  of Plates and Shells*, 2nd ed., McGraw-Hill, \u00a730 (Navier\n  series), \u00a731 (centre-deflection tables).\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n  *Concepts and Applications of Finite Element Analysis*, 4th\n  ed., Wiley, \u00a712.5.\n* Simo, J. C. and Rifai, M. S. (1990) \"A class of mixed\n  assumed-strain methods and the method of incompatible\n  modes,\" *IJNME* 29 (8), 1595\u20131638 (HEX8 enhanced strain).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport math\n\nimport numpy as np\nimport pyvista as pv\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problem data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "E = 2.0e11  # Pa (steel)\nNU = 0.3\nRHO = 7850.0  # kg/m^3\na = 1.0  # plate edge length [m]\nh = 0.02  # thickness [m]; L/h = 50\nq = 1.0e5  # uniform pressure [Pa]\n\nNX = NY = 30\nNZ = 2\n\nD = E * h**3 / (12.0 * (1.0 - NU**2))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Closed-form Navier series for the centre deflection\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def navier_w_max(q: float, a: float, D: float, n_terms: int = 25) -> float:\n    \"\"\"Sum the Navier double-sine series for the centre deflection.\"\"\"\n    total = 0.0\n    for mm in range(1, 2 * n_terms + 1, 2):\n        for nn in range(1, 2 * n_terms + 1, 2):\n            sign = ((-1) ** ((mm - 1) // 2)) * ((-1) ** ((nn - 1) // 2))\n            denom = mm * nn * (mm**2 / a**2 + nn**2 / a**2) ** 2\n            total += sign / denom\n    return (16.0 * q / (math.pi**6 * D)) * total\n\n\nw_max_published = navier_w_max(q, a, D)\nprint(f\"a = {a} m, h = {h} m, L/h = {a / h:.0f}\")\nprint(f\"E = {E / 1e9:.0f} GPa, \u03bd = {NU}, D = {D:.3e} N m\")\nprint(f\"q = {q / 1e3:.1f} kPa\")\nprint(f\"w_max  (Navier sum, 25 odd modes per axis) = {w_max_published * 1e6:.3f} \u00b5m\")\nprint(f\"w_max  (Timoshenko Table 8 \u2248 0.00406 q a^4/D) = {0.00406 * q * a**4 / D * 1e6:.3f} \u00b5m\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 30 \u00d7 30 \u00d7 2 HEX8 mesh\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "xs = np.linspace(0.0, a, NX + 1)\nys = np.linspace(0.0, a, NY + 1)\nzs = np.linspace(0.0, h, NZ + 1)\ngrid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\nprint(f\"mesh: {grid.n_points} nodes, {grid.n_cells} HEX8 cells\")\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n    ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary conditions\n\nUZ pinned along all four edges (full thickness) \u2192 SS support;\nin-plane rigid-body translation pinned at one corner;\nresidual UY rotation pinned at the diagonally-opposite corner.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pts = np.asarray(m.grid.points)\ntol = 1e-9\nedge = (\n    (np.abs(pts[:, 0]) < tol)\n    | (np.abs(pts[:, 0] - a) < tol)\n    | (np.abs(pts[:, 1]) < tol)\n    | (np.abs(pts[:, 1] - a) < tol)\n)\nm.fix(nodes=(np.where(edge)[0] + 1).tolist(), dof=\"UZ\")\n\ncorner = int(np.where(np.linalg.norm(pts, axis=1) < tol)[0][0])\nm.fix(nodes=[corner + 1], dof=\"UX\")\nm.fix(nodes=[corner + 1], dof=\"UY\")\n\ncorner_x = int(\n    np.where((np.abs(pts[:, 0] - a) < tol) & (np.abs(pts[:, 1]) < tol) & (np.abs(pts[:, 2]) < tol))[\n        0\n    ][0]\n)\nm.fix(nodes=[corner_x + 1], dof=\"UY\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Lumped pressure on the top face\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "top = np.where(np.abs(pts[:, 2] - h) < tol)[0]\nf_per_node = -q * (a * a) / len(top)\nfor n_idx in top:\n    m.apply_force(int(n_idx + 1), fz=f_per_node)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve + centre-deflection extraction\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\nu = np.asarray(res.displacement, dtype=np.float64).reshape(-1, 3)\ncentre = np.array([a / 2, a / 2, h / 2])\nidx_centre = int(np.argmin(np.linalg.norm(pts - centre, axis=1)))\nw_computed = -float(u[idx_centre, 2])\n\nerr_pct = (w_computed - w_max_published) / w_max_published * 100.0\nprint()\nprint(f\"w_max computed (HEX8 EAS, 30\u00d730\u00d72) = {w_computed * 1e6:+.3f} \u00b5m\")\nprint(f\"w_max published (Navier)            = {w_max_published * 1e6:+.3f} \u00b5m\")\nprint(f\"relative error                      = {err_pct:+.2f} %\")\n\n# 10% tolerance \u2014 HEX8 EAS at L/h = 50 lands within ~5\u20136% on a 30\u00d730\u00d72\n# mesh.  Refining to 60\u00d760\u00d72 closes this further; the validation suite\n# studies the convergence ladder explicitly.\nassert abs(err_pct) < 10.0, f\"w_max deviation {err_pct:.2f}% exceeds 10%\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deflected plate, coloured by UZ\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid_render = m.grid.copy()\ngrid_render.point_data[\"displacement\"] = u\ngrid_render.point_data[\"UZ\"] = u[:, 2]\n\nwarp = grid_render.warp_by_vector(\"displacement\", factor=2.0e2)\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 480))\nplotter.add_mesh(\n    grid_render,\n    style=\"wireframe\",\n    color=\"grey\",\n    opacity=0.35,\n    line_width=1,\n    label=\"undeformed\",\n)\nplotter.add_mesh(\n    warp,\n    scalars=\"UZ\",\n    cmap=\"viridis\",\n    show_edges=False,\n    scalar_bar_args={\"title\": \"UZ [m]  (deformed \u00d7200)\"},\n)\nplotter.view_isometric()\nplotter.camera.zoom(1.05)\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
}