{
  "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# Cantilever under linearly-varying distributed load (triangular)\n\nA clamped-free prismatic beam of length $L$ carrying a\n**triangular distributed load** that grows from zero at the root to\npeak intensity $q_{0}$ at the tip:\n\n\\begin{align}w(x) \\;=\\; q_{0}\\, \\frac{x}{L},\n    \\qquad 0 \\le x \\le L.\\end{align}\n\nThe total downward load is $Q = q_{0}\\, L / 2$ with centroid\nat $x = 2 L / 3$.  Integrating shear and moment,\n\n\\begin{align}:label: cantilever-triangular-moment\n\n    M(x) \\;=\\; -\\,\\frac{q_{0}}{6\\, L}\\,(L - x)^{2}\\,(2 L + x),\\end{align}\n\nso the moment at the root takes its closed-form value\n\n\\begin{align}M_{\\mathrm{root}} \\;=\\; -\\,\\frac{q_{0}\\, L^{2}}{3}.\\end{align}\n\nTip deflection by virtual work (unit-load applied at the tip,\n$m(x) = -(L - x)$):\n\n\\begin{align}:label: cantilever-triangular-deflection\n\n    \\delta_{\\mathrm{tip}}\n    \\;=\\;\n    \\frac{1}{E\\, I}\n    \\int_{0}^{L} M(x)\\, m(x)\\, \\mathrm{d}x\n    \\;=\\;\n    \\frac{11\\, q_{0}\\, L^{4}}{120\\, E\\, I},\\end{align}\n\n(Timoshenko 1955 \u00a732; Roark Table 8 case 5 \u2014 load distribution\npeaking at the free end).  Tip slope ($\\theta_{\\mathrm{tip}}\n= q_{0}\\, L^{3} / (8\\, E\\, I)$) and the deflection at any interior\npoint follow by the same integration.\n\nComparison to a uniform load: if instead the same total load\n$Q = q_{0} L / 2$ is spread *uniformly*\n($\\bar w = q_{0}/2$), the tip deflection is\n\n\\begin{align}\\delta_{\\mathrm{tip}}^{\\mathrm{UDL}}\n    \\;=\\;\n    \\frac{\\bar w\\, L^{4}}{8\\, E\\, I}\n    \\;=\\;\n    \\frac{q_{0}\\, L^{4}}{16\\, E\\, I},\\end{align}\n\nso $\\delta_{\\mathrm{tip}}^{\\mathrm{UDL}} /\n\\delta_{\\mathrm{tip}}^{\\mathrm{tri}} = 120 / (16 \\cdot 11) =\n15/22 \\approx 0.682$.  Concentrating the load toward the tip\ndeflects the beam **~46 % more** than spreading it uniformly,\neven though the total load is identical \u2014 the moment-arm\ndistribution dominates.\n\n## Implementation\n\nA 40-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam.\nThe triangular load is converted to consistent nodal forces by\ntrapezoidal integration of $w(x)$ over each element edge\n(interior nodes receive $w(x_i)\\, h$, the end nodes\n$w(x_i)\\, h / 2$); the small nodal-discretisation error this\nintroduces falls below 0.1 % on a 40-element mesh.\n\n## References\n\n* Timoshenko, S. P. (1955) *Strength of Materials, Part I:\n  Elementary Theory and Problems*, 3rd ed., Van Nostrand, \u00a732\n  (cantilever with distributed load).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 5\n  (linearly varying distributed load on a cantilever beam).\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); \u00a72.10 (consistent nodal loads).\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_section = 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 [m]\nq_0 = 2.0e3  # peak load intensity [N/m, downward = -y]\n\nEI = E * I_z\n\n# Closed-form references (Roark Table 8 case 5).\nM_root_pub = -q_0 * L**2 / 3.0\ndelta_tip_pub = -11.0 * q_0 * L**4 / (120.0 * EI)\ntheta_tip_pub = -q_0 * L**3 / (8.0 * EI)\n\n# Comparison against the same total load Q = q_0 L / 2 spread uniformly.\ndelta_tip_udl = -q_0 * L**4 / (16.0 * EI)\n\nprint(\"Cantilever under linearly-varying triangular distributed load\")\nprint(f\"  L = {L} m, peak intensity q_0 = {q_0} N/m at the free end\")\nprint(f\"  total load Q = q_0 L / 2 = {q_0 * L / 2:.1f} N (centroid at x = 2 L / 3)\")\nprint(f\"  E = {E:.2e} Pa, I = {I_z:.3e} m^4, EI = {EI:.3e} N m^2\")\nprint()\nprint(\"Closed-form references (Roark Table 8 case 5):\")\nprint(f\"  M_root      = {M_root_pub:+.4e} N m   (= -q_0 L^2 / 3)\")\nprint(f\"  v_tip       = {delta_tip_pub * 1e3:+.4e} mm   (= -11 q_0 L^4 / 120 E I)\")\nprint(f\"  theta_tip   = {theta_tip_pub:+.4e} rad   (= -q_0 L^3 / 8 E I)\")\nprint()\nprint(f\"  v_tip with same Q spread uniformly = {delta_tip_udl * 1e3:+.4e} mm\")\nprint(\n    f\"  \u2192 triangular peak-at-tip deflects \"\n    f\"{abs(delta_tip_pub / delta_tip_udl):.3f}\u00d7 the UDL value (= 22/15 \u2248 1.467)\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build a 40-element BEAM2 line\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_ELEM = 40\nxs = np.linspace(0.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_section, I_z, I_y, J),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary conditions: full clamp at x = 0\n\nOut-of-plane motion (UZ, ROTX, ROTY) suppressed across the line so\nthe response stays strictly 2D in the x-y plane.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m.fix(nodes=1, dof=\"ALL\")\nfor i in range(N_ELEM + 1):\n    m.fix(nodes=int(i + 1), dof=\"UZ\")\n    m.fix(nodes=int(i + 1), dof=\"ROTX\")\n    m.fix(nodes=int(i + 1), dof=\"ROTY\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Apply the triangular load as consistent nodal forces\n\nEach interior node ``i`` receives ``w(x_i) \u00b7 h`` (the trapezoidal\nintegral of ``w`` over its half-element on each side); the two\nextremes get ``w(x) \u00b7 h / 2``.  ``w(x) = q_0 \u00b7 x / L`` so the\nroot contributes zero, and the tip gets the largest force.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "h = L / N_ELEM\nfor i, x in enumerate(xs):\n    w = q_0 * x / L\n    weight = w * h / 2.0 if (i == 0 or i == N_ELEM) else w * h\n    m.apply_force(int(i + 1), fy=-weight)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve + extract tip deflection / slope\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    \"\"\"0=UX, 1=UY, 2=UZ, 3=ROTX, 4=ROTY, 5=ROTZ.\"\"\"\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\nv_tip_fe = _node_dof(N_ELEM + 1, 1)\ntheta_tip_fe = _node_dof(N_ELEM + 1, 5)\n\nerr_v = (abs(v_tip_fe) - abs(delta_tip_pub)) / abs(delta_tip_pub) * 100.0\nerr_theta = (abs(theta_tip_fe) - abs(theta_tip_pub)) / abs(theta_tip_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_tip':<22}  {v_tip_fe * 1e3:>10.4f} mm  {delta_tip_pub * 1e3:>10.4f} mm  {err_v:>+8.4f}%\"\n)\nprint(\n    f\"{'|theta_tip|':<22}  {abs(theta_tip_fe):>11.6e}  \"\n    f\"{abs(theta_tip_pub):>11.6e}  {err_theta:>+8.4f}%\"\n)\n\n# Lumped consistent-load mapping introduces a small discretisation\n# error (O(h\u00b2)); 40 elements lands well within 0.5 %.\nassert abs(err_v) < 0.5, f\"v_tip deviation {err_v:.4f} % exceeds 0.5 %\"\nassert abs(err_theta) < 0.5, f\"theta_tip deviation {err_theta:.4f} % exceeds 0.5 %\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the deflection profile against the closed form\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def v_closed_form(x: float) -> float:\n    \"\"\"Closed-form deflection.\n\n    Integrate :math:`E I\\\\,v'' = M(x)` twice with\n    :math:`v(0) = v'(0) = 0`.  After algebra,\n\n    .. math::\n        E I\\\\,v(x)\n        = -\\\\frac{q_{0}}{120 L}\\\\, x^{2}\\\\,(20 L^{3} - 10 L^{2} x + x^{3}).\n    \"\"\"\n    return -q_0 * x**2 * (20.0 * L**3 - 10.0 * L**2 * x + x**3) / (120.0 * L * EI)\n\n\nfe_uy = np.array([_node_dof(int(i + 1), 1) for i in range(N_ELEM + 1)])\n\nimport matplotlib.pyplot as plt  # noqa: E402\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11.0, 4.0))\n\nxs_dense = np.linspace(0.0, L, 401)\nax1.plot(\n    xs_dense,\n    [v_closed_form(x) * 1e3 for x in xs_dense],\n    \"-\",\n    color=\"#1f77b4\",\n    lw=2.0,\n    label=\"closed form\",\n)\nax1.plot(xs, fe_uy * 1e3, \"o\", color=\"#d62728\", markersize=4, label=\"BEAM2 FE\")\nax1.set_xlabel(\"x [m]\")\nax1.set_ylabel(\"v(x)  [mm]   (downward negative)\")\nax1.set_title(\"Deflection profile\", fontsize=11)\nax1.legend(loc=\"lower left\", fontsize=9)\nax1.grid(True, ls=\":\", alpha=0.5)\n\n# Show the load distribution alongside.\nax2.fill_between(xs_dense, 0.0, q_0 * xs_dense / L, alpha=0.4, color=\"#2ca02c\")\nax2.plot(xs_dense, q_0 * xs_dense / L, color=\"#2ca02c\", lw=2.0, label=\"w(x) = q_0 x / L\")\nax2.set_xlabel(\"x [m]\")\nax2.set_ylabel(\"w(x)  [N/m]\")\nax2.set_title(\"Triangular load distribution\", fontsize=11)\nax2.legend(loc=\"upper left\", fontsize=9)\nax2.grid(True, ls=\":\", alpha=0.5)\nfig.tight_layout()\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deformed beam\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_tip_fe), 1e-12) * 0.05\nwarped = g.copy()\nwarped.points = np.asarray(g.points) + scale * disp_xyz\n\nplotter = pv.Plotter(off_screen=True, window_size=(900, 320))\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)\nplotter.add_points(\n    np.array([[0.0, 0.0, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"black\",\n    label=\"clamp\",\n)\nplotter.add_points(\n    np.array([[L, scale * v_tip_fe, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#d62728\",\n    label=f\"tip \u2014 v_tip = {v_tip_fe * 1e3:.4f} mm\",\n)\nplotter.view_xy()\nplotter.add_legend()\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Take-aways\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(\"Take-aways:\")\nprint(f\"  \u2022 v_tip within {abs(err_v):.4f} % of -11 q_0 L\u2074 / (120 E I).\")\nprint(f\"  \u2022 |theta_tip| within {abs(err_theta):.4f} % of q_0 L\u00b3 / (8 E I).\")\nprint(\n    f\"  \u2022 Same total load Q spread uniformly deflects the tip by only \"\n    f\"{abs(delta_tip_udl / delta_tip_pub):.3f}\u00d7 the triangular case \u2014 \"\n    \"moment-arm matters more than total load magnitude on a slender cantilever.\"\n)"
      ]
    }
  ],
  "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
}