{
  "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# Continuous beam over three supports \u2014 Clapeyron's three-moment theorem\n\nClassical statically-indeterminate beam benchmark: a prismatic\nstraight beam of total length $2L$ rests on three simple\nsupports \u2014 two at the ends and one at the centre \u2014 and carries a\n**uniformly distributed load** $q$ over both spans.  Because\nthere is one redundant reaction, the moment at the central support\nfollows from the three-moment theorem of Clapeyron (1857):\n\n\\begin{align}M_{\\mathrm{mid}} \\;=\\; -\\,\\frac{q\\, L^{2}}{8}.\\end{align}\n\nClosed-form deflections, reactions, and span moments follow by\nintegrating $EI\\, v'' = M(x)$ with $M(x) = R_e\\, x -\nq\\, x^{2}/2$ on the left span (Timoshenko 1955 \u00a717; Roark Table 8\ncase 9):\n\n\\begin{align}R_e = \\tfrac{3}{8}\\, q\\, L,\\qquad\n    R_{\\mathrm{mid}} = \\tfrac{5}{4}\\, q\\, L,\n    \\qquad\n    M_{\\max}^{(+)} = \\tfrac{9}{128}\\, q\\, L^{2}\\ \\text{at}\\ x = \\tfrac{3}{8}\\, L,\\end{align}\n\nand\n\n\\begin{align}:label: continuous-beam-deflection\n\n    v(x) \\;=\\;\n    -\\,\\frac{q\\, L^{3}}{48\\, E\\, I}\\, x\n    \\;+\\;\n    \\frac{q\\, L\\, x^{3}}{16\\, E\\, I}\n    \\;-\\;\n    \\frac{q\\, x^{4}}{24\\, E\\, I},\\end{align}\n\nso two clean check points are available:\n\n* mid-span ($x = L/2$):\n  $v_{\\mathrm{mid}} \\;=\\; -\\, q\\, L^{4} / (192\\, E\\, I)$.\n* peak-deflection point ($x \\approx 0.4215\\, L$, the positive\n  root of $8 u^{3} - 9 u^{2} + 1 = 0$):\n  $v_{\\max} \\;\\approx\\; -\\, q\\, L^{4} / (185\\, E\\, I)$.\n\nBy symmetry the right span carries identical reactions and an\nidentical mirror-image deflection profile.\n\n## Implementation\n\nA 60-element BEAM2 (Hermite-cubic Bernoulli) line spans the full\n$2L$ beam \u2014 30 elements per span.  The UDL $q$ is\nconverted to consistent nodal forces (interior nodes receive\n$q\\, h$, edge nodes $q\\, h / 2$, where $h = L /\n30$ is the element length).  Three nodal supports pin $u_y =\n0$ at $x = 0$, $x = L$, and $x = 2L$; out-of-\nplane / rotational DOFs are pinned at one end node to remove the\nremaining rigid-body modes without introducing spurious moments.\n\nHermite-cubic shape functions interpolate the deflection\nanalytically between nodes, so the FE response under a piecewise-\nlinear distributed load matches the closed form to machine\nprecision at every node \u2014 well within the 0.5 % tolerance the\nverification framework specifies.\n\n## References\n\n* Clapeyron, B. P. E. (1857) \"Calcul d'une poutre \u00e9lastique\n  reposant librement sur des appuis in\u00e9galement espac\u00e9s,\"\n  *Comptes Rendus de l'Acad\u00e9mie des Sciences* 45, 1076\u20131080 \u2014\n  the original three-moment theorem.\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I:\n  Elementary Theory and Problems*, 3rd ed., Van Nostrand,\n  \u00a717 (continuous beams + three-moment equation).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 9\n  (continuous beam, equal spans, uniform load).\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.\n  (2002) *Concepts and Applications of Finite Element\n  Analysis*, 4th ed., Wiley, \u00a72.5 (Hermite cubic shape\n  functions); \u00a716.2 (statically indeterminate beam FE).\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\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\nNU = 0.30\nRHO = 7850.0\nWIDTH = 0.05\nHEIGHT = 0.05\nA = WIDTH * HEIGHT\nI_z = WIDTH * HEIGHT**3 / 12.0\nI_y = HEIGHT * WIDTH**3 / 12.0\nJ = (1.0 / 3.0) * min(WIDTH, HEIGHT) ** 3 * max(WIDTH, HEIGHT)\n\nL = 1.0  # span length [m]; total beam length = 2 L\nq = 1.0e3  # uniform distributed load magnitude [N/m, downward = -y]\n\n# Closed-form references (Roark Table 8 case 9 / Timoshenko \u00a717).\nEI = E * I_z\nM_mid_pub = -q * L**2 / 8.0  # moment at central support\nR_e_pub = 3.0 * q * L / 8.0  # reaction at end support\nR_mid_pub = 5.0 * q * L / 4.0  # reaction at central support\nv_midspan_pub = -q * L**4 / (192.0 * EI)  # deflection at x = L/2\n# Peak-deflection: factor 8u^3 - 9u^2 + 1 = (u - 1)(8u^2 - u - 1).\n# The (u - 1) root corresponds to the end of the span; the\n# physically meaningful peak sits at the positive root of the\n# quadratic, ``u = (1 + \u221a33) / 16 \u2248 0.4216``.\nu_peak_root = (1.0 + 33.0**0.5) / 16.0\nv_peak_pub = (\n    q * L * (u_peak_root * L) ** 3 / (16.0 * EI)\n    - q * (u_peak_root * L) ** 4 / (24.0 * EI)\n    - q * L**3 * (u_peak_root * L) / (48.0 * EI)\n)\n\nprint(\"Continuous beam \u2014 three equal supports, both spans uniformly loaded\")\nprint(f\"  span L = {L} m, total length = {2 * L} m, q = {q} N/m down\")\nprint(f\"  E = {E:.2e} Pa, I = {I_z:.3e} m^4, EI = {EI:.3e} N m^2\")\nprint()\nprint(\"Closed-form references (Clapeyron / Roark Table 8 case 9):\")\nprint(f\"  M_middle      = {M_mid_pub:+.4e} N m   (= -q L^2 / 8)\")\nprint(f\"  R_end         = {R_e_pub:+.4e} N      (= 3 q L / 8)\")\nprint(f\"  R_middle      = {R_mid_pub:+.4e} N      (= 5 q L / 4)\")\nprint(f\"  v(L/2)        = {v_midspan_pub * 1e3:+.4e} mm    (= -q L^4 / (192 E I))\")\nprint(\n    f\"  v_peak        = {v_peak_pub * 1e3:+.4e} mm    \"\n    f\"(at x = {u_peak_root:.4f} L; \u2248 -q L^4/(185 E I))\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 60-element BEAM2 line\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_PER_SPAN = 30\nn_elem = 2 * N_PER_SPAN\nxs = np.linspace(0.0, 2.0 * L, n_elem + 1)\npts = np.column_stack((xs, np.zeros_like(xs), np.zeros_like(xs)))\n\ncells = np.empty((n_elem, 3), dtype=np.int64)\ncells[:, 0] = 2\ncells[:, 1] = np.arange(n_elem)\ncells[:, 2] = np.arange(1, n_elem + 1)\ngrid = pv.UnstructuredGrid(\n    cells.ravel(),\n    np.full(n_elem, 3, dtype=np.int64),\n    pts,\n)\n\nm = femorph_solver.Model.from_grid(grid)\nm.assign(\n    ELEMENTS.BEAM2,\n    material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n    real=(A, I_z, I_y, J),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary conditions: simple supports at x = 0, L, 2L\n\nEach support pins UY (the load direction) \u2014 that's all a *simple*\nsupport requires.  In addition we pin UZ + ROTX everywhere along\nthe line to keep the problem strictly 2D, and pin UX + ROTY +\nROTZ at one node to remove rigid-body translation along x and\nrotation about y / z.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "end_left = 1\nend_right = n_elem + 1\nmid_node = N_PER_SPAN + 1\n\nm.fix(nodes=int(end_left), dof=\"UY\")\nm.fix(nodes=int(mid_node), dof=\"UY\")\nm.fix(nodes=int(end_right), dof=\"UY\")\n\n# Suppress out-of-plane motion + axial / rotational rigid bodies.\n# Apply UZ = 0, ROTX = 0 to every node \u2014 keeps the problem 2D in\n# the x-y plane.\nfor nn in range(1, n_elem + 2):\n    m.fix(nodes=int(nn), dof=\"UZ\")\n    m.fix(nodes=int(nn), dof=\"ROTX\")\n\n# Single anchor at end_left for UX and the in-plane rotations\n# perpendicular to the bending plane \u2014 kills the remaining\n# rigid-body modes.\nm.fix(nodes=int(end_left), dof=\"UX\")\nm.fix(nodes=int(end_left), dof=\"ROTY\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Apply the UDL as consistent nodal forces\n\nFor Hermite-cubic BEAM2 the work-equivalent nodal load from a\nconstant UDL ``q`` over an element of length ``h`` is ``q h``\nsplit evenly across the two end nodes (``q h / 2`` each).\nInterior nodes therefore receive ``q h`` total (the contributions\nfrom both adjacent elements add); the two beam ends receive\n``q h / 2``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "h = L / N_PER_SPAN  # element length\nfor i in range(n_elem + 1):\n    # Pre-compute weight: q*h for interior nodes, q*h/2 at the\n    # two extreme ends.\n    weight = q * h / 2.0 if i == 0 or i == n_elem else q * h\n    m.apply_force(int(i + 1), fy=-weight)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve + extract\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()\ndof_map = m.dof_map()\ndisp = np.asarray(res.displacement, dtype=np.float64)\n\n\ndef _node_dof(node_id: int, dof_idx: int) -> float:\n    \"\"\"Return the displacement of ``dof_idx`` (0=UX, 1=UY, 2=UZ) at ``node_id``.\"\"\"\n    rows = np.where(dof_map[:, 0] == node_id)[0]\n    for r in rows:\n        if int(dof_map[r, 1]) == dof_idx:\n            return float(disp[r])\n    return 0.0\n\n\n# Mid-span deflection at x = L/2 (node L/2 / h + 1).\ni_mid = N_PER_SPAN // 2\nv_midspan_fe = _node_dof(int(i_mid + 1), 1)\n\n# Peak deflection: scan the left span for the most negative UY.\nleft_span_uys = np.array([_node_dof(int(i + 1), 1) for i in range(N_PER_SPAN + 1)])\ni_peak = int(np.argmin(left_span_uys))\nv_peak_fe = float(left_span_uys[i_peak])\nx_peak_fe = float(xs[i_peak])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify deflection check points\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "err_mid = (v_midspan_fe - v_midspan_pub) / abs(v_midspan_pub) * 100.0\nerr_peak = (v_peak_fe - v_peak_pub) / abs(v_peak_pub) * 100.0\n\nprint()\nprint(f\"{'quantity':<22}  {'FE':>14}  {'published':>14}  {'err':>9}\")\nprint(f\"{'-' * 22:<22}  {'-' * 14:>14}  {'-' * 14:>14}  {'-' * 9:>9}\")\nprint(\n    f\"{'v(L/2)':<22}  {v_midspan_fe * 1e3:>10.4f} mm  \"\n    f\"{v_midspan_pub * 1e3:>10.4f} mm  {err_mid:>+8.3f}%\"\n)\nprint(\n    f\"{'v_peak':<22}  {v_peak_fe * 1e3:>10.4f} mm  {v_peak_pub * 1e3:>10.4f} mm  {err_peak:>+8.3f}%\"\n)\nprint(f\"  peak location FE = {x_peak_fe:.4f} m  vs  pub = {u_peak_root * L:.4f} m\")\n\nassert abs(err_mid) < 0.1, f\"v(L/2) deviation {err_mid:.3f} % exceeds 0.1 %\"\nassert abs(err_peak) < 0.5, f\"v_peak deviation {err_peak:.3f} % exceeds 0.5 %\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deformed beam\n\nMagnify the displacement so the deflection is visible against\nthe 2-m span.  Mark the three supports and the peak-deflection\npoint in each span.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = m.grid.copy()\ndisp_xyz = np.zeros((g.n_points, 3))\nfor ni in range(g.n_points):\n    disp_xyz[ni, 1] = _node_dof(int(ni + 1), 1)\ng.point_data[\"displacement\"] = disp_xyz\ng.point_data[\"UY\"] = disp_xyz[:, 1]\n\nscale = 1.0 / max(abs(v_peak_fe), 1e-12) * 0.05  # scale so peak \u2248 5 cm visually\nwarped = g.copy()\nwarped.points = np.asarray(g.points) + scale * disp_xyz\n\nplotter = pv.Plotter(off_screen=True, window_size=(900, 360))\nplotter.add_mesh(g, color=\"grey\", opacity=0.4, line_width=2, label=\"undeformed\")\nplotter.add_mesh(warped, scalars=\"UY\", cmap=\"coolwarm\", line_width=4)\nsupport_pts = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0], [2 * L, 0.0, 0.0]])\nplotter.add_points(\n    support_pts,\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"black\",\n    label=\"supports\",\n)\nplotter.add_points(\n    np.array([[x_peak_fe, scale * v_peak_fe, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#d62728\",\n    label=f\"v_peak FE = {v_peak_fe * 1e3:.3f} mm\",\n)\nplotter.view_xy()\nplotter.add_legend()\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Closing notes\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(\"Take-aways:\")\nprint(f\"  \u2022 v(L/2) within {abs(err_mid):.3f} % of the closed form.\")\nprint(f\"  \u2022 v_peak within {abs(err_peak):.3f} % at x \u2248 {x_peak_fe:.4f} m.\")\nprint(\"  \u2022 The Hermite-cubic BEAM2 element interpolates Bernoulli kinematics\")\nprint(\"    *exactly* on the polynomial moment field \u2014 only the consistent-load\")\nprint(\"    quadrature introduces any discretisation error here.\")"
      ]
    }
  ],
  "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
}