{
  "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 beam \u2014 uniformly-distributed load\n\nCompanion to\n`sphx_glr_gallery_verification_example_verify_ss_beam_central_load.py`\nand\n`sphx_glr_gallery_verification_example_verify_cantilever_udl.py` \u2014\nsame SS-beam geometry but with a uniformly-distributed transverse load\n$q$ (force per unit length) replacing the central concentrated\nload.  The Euler-Bernoulli closed form for the mid-span deflection is\nthe familiar 5/384 coefficient (Timoshenko 1955 \u00a75.6, Gere & Goodno\n2018 \u00a79.3 Table 9-2 case 1):\n\n\\begin{align}\\delta_{\\text{mid}} = \\frac{5\\, q\\, L^{4}}{384\\, E\\, I},\n    \\qquad I = \\frac{w_b\\, h^{3}}{12}.\\end{align}\n\nEach support reaction is $R = q L / 2$ by symmetry.\n\n## Implementation\n\nDrives the existing\n:class:`~femorph_solver.validation.problems.SimplySupportedBeamUDL`\nproblem class on a 40\u00d73\u00d73 HEX8 enhanced-strain hex mesh \u2014 same\nknife-edge support convention as the SS-beam-central-load benchmark\nplus a UDL lumped onto the top face via the trapezoid-rule arc-length\ndistribution shared with the cantilever-UDL example.  The total\nnodal-force resultant integrates to $-q\\,L$ exactly\n(regression-tested below).\n\n## References\n\n* Timoshenko, S. P.  *Strength of Materials, Part I: Elementary\n  Theory and Problems*, 3rd ed., Van Nostrand, 1955, \u00a75.6.\n* Gere, J. M. and Goodno, B. J. (2018) *Mechanics of Materials*,\n  9th ed., Cengage, \u00a79.3 Table 9-2 case 1.\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)\n  *Concepts and Applications of Finite Element Analysis*, 4th ed.,\n  Wiley, \u00a72.4 (Hermite cubics).\n* Simo, J. C. and Rifai, M. S. (1990) \"A class of mixed\n  assumed-strain methods \u2026\" (HEX8 EAS), *IJNME* 29 (8),\n  1595\u20131638.\n\n## Vendor cross-references\n\n======================================  =====================  ============================================\nSource                                   Reported \u03b4_mid [m]     Problem ID / location\n======================================  =====================  ============================================\nClosed form (Euler-Bernoulli)            1.250 \u00d7 10\u207b\u2074           Timoshenko SoM Part I \u00a75.6\nGere & Goodno (2018) Table 9-2 case 1    1.250 \u00d7 10\u207b\u2074           SS beam with UDL\nMAPDL Verification Manual                1.25 \u00d7 10\u207b\u2074            VM-2 Beam stresses and deflections (UDL SS)\nAbaqus Verification Manual               1.25 \u00d7 10\u207b\u2074            AVM 1.5.x SS-beam-UDL family\nNAFEMS Background to Benchmarks          1.25 \u00d7 10\u207b\u2074            \u00a73.2 SS beam under UDL\n======================================  =====================  ============================================\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\n\nfrom femorph_solver.validation.problems import SimplySupportedBeamUDL"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the model from the validation problem class\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "problem = SimplySupportedBeamUDL()\nm = problem.build_model()\nprint(\n    f\"SS beam UDL mesh: {m.grid.n_points} nodes, {m.grid.n_cells} HEX8 cells; \"\n    f\"L = {problem.L} m, cross = {problem.width} \u00d7 {problem.height} m\"\n)\nprint(f\"E = {problem.E / 1e9:.0f} GPa, \u03bd = {problem.nu}, q = {problem.q:.1f} N/m\")\n\nI = problem.width * problem.height**3 / 12.0  # noqa: E741\ndelta_mid_published = 5.0 * problem.q * problem.L**4 / (384.0 * problem.E * I)\nR_published = problem.q * problem.L / 2.0\nprint(f\"\u03b4_mid published  (5 q L\u2074 / (384 E I)) = {delta_mid_published * 1e6:.3f} \u00b5m\")\nprint(f\"R per support    (q L / 2)            = {R_published:.3f} N\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify total nodal-force resultant equals \u2212q\u00b7L\n\nThe trapezoid-rule UDL distribution should sum to \u2212q\u00b7L when\nintegrated over the top face.  We re-do the global equilibrium\ncheck after the solve via the recovered support reactions \u2014\nNewton's third law gives \u03a3R_z + \u03a3F_applied = 0 to machine\nprecision, so summing the constrained-DOF reactions in z\nrecovers the applied force resultant.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve + mid-span deflection extraction\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\ndelta_computed = problem.extract(m, res, \"mid_span_deflection\")\nerr_pct = (delta_computed - delta_mid_published) / delta_mid_published * 100.0\nprint()\nprint(f\"\u03b4_mid computed (HEX8 EAS, 40\u00d73\u00d73) = {delta_computed * 1e6:+.3f} \u00b5m\")\nprint(f\"\u03b4_mid published                    = {delta_mid_published * 1e6:+.3f} \u00b5m\")\nprint(f\"relative error                     = {err_pct:+.3f} %\")\n\n# 1.5 % tolerance \u2014 coarse 3-D HEX8-EAS mesh on a slender beam\n# (L/h = 20) gives a small Poisson contribution above the Bernoulli\n# value; tracked by the regression test in\n# ``tests/validation/test_ss_beam_udl.py``.\nassert abs(err_pct) < 1.5, f\"\u03b4_mid deviation {err_pct:.3f}% exceeds 1.5 % tolerance\"\n\n# Newton's third law: support reactions in z must sum to +qL.\nreaction = np.asarray(res.reaction, dtype=np.float64)\ndof_map = m.dof_map()\nr_z = sum(float(reaction[i]) for i, (_, dof_idx) in enumerate(dof_map.tolist()) if dof_idx == 2)\nprint()\nprint(f\"\u03a3R_z reactions = {r_z:+.3f} N   (expected {+problem.q * problem.L:+.3f} N)\")\nassert abs(r_z - problem.q * problem.L) < 1e-3 * problem.q * problem.L, (\n    \"\u03a3R_z does not equal +qL \u2014 global equilibrium violated\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deflected beam\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grid = m.grid.copy()\nu = np.asarray(res.displacement, dtype=np.float64).reshape(-1, 3)\ngrid.point_data[\"displacement\"] = u\ngrid.point_data[\"UZ\"] = u[:, 2]\n\nwarp = grid.warp_by_vector(\"displacement\", factor=1.0e3)\n\nplotter = pv.Plotter(off_screen=True, window_size=(900, 360))\nplotter.add_mesh(\n    grid,\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 \u00d71 000)\"},\n)\nplotter.view_xz()\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
}