{
  "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# L-shaped frame under tip load \u2014 Castigliano on a two-member portal\n\nTwo prismatic Bernoulli beams welded at a rigid corner form an\nL-shaped portal frame.  A vertical column of length $L_{v}$\nruns from the clamp at the origin to the corner; a horizontal beam\nof length $L_{h}$ cantilevers from the corner.  A\nconcentrated transverse load $P$ acts at the free tip\n$(L_{h}, L_{v})$ in the $-y$ direction (downward).\n\nReading off the moments by superposition:\n\n* **Column** (vertical, $0 \\le y \\le L_{v}$): axial force\n  $N = -P$ (compression); the tip load applied at horizontal\n  offset $L_{h}$ produces a constant moment about $z$,\n\n  .. math::\n     M_{\\mathrm{col}}(y) \\;=\\; -\\, P\\, L_{h}.\n\n* **Horizontal beam** ($0 \\le x \\le L_{h}$): standard\n  cantilever-with-tip-load moment,\n\n  .. math::\n     M_{\\mathrm{beam}}(x) \\;=\\; -\\, P\\, (L_{h} - x).\n\nCastigliano's theorem gives the tip deflection in the load\ndirection (Timoshenko & Young 1968 \u00a780; Roark Table 9 case 6):\n\n\\begin{align}:label: l-frame-tip-deflection\n\n    \\delta_{\\mathrm{tip}}\n    \\;=\\;\n    \\frac{1}{E\\, I}\n    \\Biggl[\n        \\int_{0}^{L_{v}} M_{\\mathrm{col}}^{2}\\, \\mathrm{d}y\n        +\n        \\int_{0}^{L_{h}} M_{\\mathrm{beam}}^{2}\\, \\mathrm{d}x\n    \\Biggr]_{ \\!\\!\\!/ P}\n    \\;=\\;\n    \\frac{P\\, L_{h}^{2}\\, L_{v}}{E\\, I}\n    \\;+\\;\n    \\frac{P\\, L_{h}^{3}}{3\\, E\\, I}\n    \\;+\\;\n    \\frac{P\\, L_{v}}{E\\, A},\\end{align}\n\nwith the axial-strain term $P L_{v} / (E A)$ from the\ncolumn.  For a slender frame $(I / A \\ll L^{2})$ the axial\ncontribution is negligible \u2014 at the default geometry below it\nsits 6 000\u00d7 below the bending term.\n\nTwo limits collapse the formula:\n\n* $L_{v} \\to 0$ \u2014 the column vanishes and the structure\n  reduces to a horizontal cantilever:\n  $\\delta = P L_{h}^{3} / (3 E I)$.  Recovers the\n  `sphx_glr_gallery_verification_example_verify_cantilever_eb.py`\n  result.\n* $L_{h} \\to 0$ \u2014 only axial compression of the column\n  survives: $\\delta = P L_{v} / (E A)$.  Recovers the\n  `sphx_glr_gallery_verification_example_verify_single_hex_uniaxial.py`\n  / Hooke's-law axial response.\n\n## Implementation\n\n40-element BEAM2 (Hermite-cubic Bernoulli) line for the column,\n40 more for the horizontal beam \u2014 sharing a common node at the\nrigid corner so all moment / shear / axial transfer is implicit.\nThe clamp at $(0, 0, 0)$ pins all six DOFs; out-of-plane\nDOFs are pinned across the line so the response stays strictly\n2D in the $x$-$y$ plane.\n\n## References\n\n* Timoshenko, S. P. and Young, D. H. (1968) *Elements of\n  Strength of Materials*, 5th ed., Van Nostrand, \u00a780\n  (Castigliano on a curved-beam / frame).\n* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for\n  Stress and Strain*, 6th ed., McGraw-Hill, Table 9 case 6\n  (right-angle bend, end 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, \u00a716.5 (multi-element frame\n  assembly).\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 \u2014 equal-leg L-frame\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_v = 1.0  # column length [m]\nL_h = 1.0  # beam length [m]\nP = 1.0e3  # tip load [N, downward = -y]\n\nEI = E * I_z\n\n# Closed-form tip deflection (Castigliano, Roark Table 9 case 6).\ndelta_bending = P * L_h**2 * L_v / EI + P * L_h**3 / (3.0 * EI)\ndelta_axial = P * L_v / (E * A_section)\ndelta_tip_pub = -(delta_bending + delta_axial)  # downward in -y\n\nprint(\"L-shaped frame under tip load\")\nprint(f\"  column L_v = {L_v} m, beam L_h = {L_h} m, P = {P} N (-y)\")\nprint(f\"  E = {E:.2e} Pa, I = {I_z:.3e} m^4, A = {A_section:.3e} m^2\")\nprint()\nprint(\"Closed-form references (Castigliano / Roark Table 9 case 6):\")\nprint(f\"  bending \u03b4           = {-delta_bending * 1e3:+.4e} mm\")\nprint(\"    = P L_h\u00b2 L_v / E I + P L_h\u00b3 / (3 E I)\")\nprint(f\"  axial \u03b4 (column)    = {-delta_axial * 1e3:+.4e} mm   (= P L_v / E A)\")\nprint(f\"  total \u03b4_tip         = {delta_tip_pub * 1e3:+.4e} mm\")\nprint(\n    f\"  axial / bending ratio = {delta_axial / delta_bending:.2e}  \"\n    \"(slender-frame regime \u2014 axial negligible)\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build the two-segment BEAM2 frame\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "N_PER_SEG = 40\n\n# Column: 0 \u2192 L_v along +y axis.\ncol_y = np.linspace(0.0, L_v, N_PER_SEG + 1)\ncol_pts = np.column_stack((np.zeros_like(col_y), col_y, np.zeros_like(col_y)))\n\n# Beam: 0 \u2192 L_h along +x axis at y = L_v.  Skip the corner node (it's\n# already in the column) and concatenate the rest.\nbeam_x = np.linspace(0.0, L_h, N_PER_SEG + 1)[1:]\nbeam_pts = np.column_stack((beam_x, np.full_like(beam_x, L_v), np.zeros_like(beam_x)))\n\npts = np.vstack((col_pts, beam_pts))\n\n# Cells \u2014 column elements then beam elements.  Both segments share\n# node ``N_PER_SEG`` (the corner), so beam-element indices have an\n# implicit offset of ``N_PER_SEG`` for their first node.\ncells = []\nfor i in range(N_PER_SEG):\n    cells.append([2, i, i + 1])\n# Beam segment: connect corner (node N_PER_SEG) to first beam point\n# (which is the point right after corner = N_PER_SEG + 1).\ncells.append([2, N_PER_SEG, N_PER_SEG + 1])\nfor i in range(1, N_PER_SEG):\n    cells.append([2, N_PER_SEG + i, N_PER_SEG + i + 1])\ncells_arr = np.array(cells, dtype=np.int64)\n\nn_cells = cells_arr.shape[0]\ngrid = pv.UnstructuredGrid(\n    cells_arr.ravel(),\n    np.full(n_cells, 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)\n\nn_nodes = pts.shape[0]\nprint(f\"\\nL-frame mesh: {n_nodes} nodes, {n_cells} BEAM2 cells\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Boundary conditions\n\nClamp at (0, 0, 0) \u2014 node 1.  Out-of-plane motion suppressed at\nevery node so the response stays strictly 2D in x-y.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m.fix(nodes=1, dof=\"ALL\")\n\nfor i in range(n_nodes):\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": [
        "## Tip load at (L_h, L_v, 0)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_tip = n_nodes - 1  # last node added (the beam tip)\nm.apply_force(int(i_tip + 1), fy=-P)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Solve + extract tip deflection\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(int(i_tip + 1), 1)\nerr = (abs(v_tip_fe) - abs(delta_tip_pub)) / abs(delta_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 (UY)':<22}  {v_tip_fe * 1e3:>10.4f} mm  \"\n    f\"{delta_tip_pub * 1e3:>10.4f} mm  {err:>+8.4f}%\"\n)\n\nassert abs(err) < 0.1, f\"v_tip deviation {err:.4f} % exceeds 0.1 %\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deformed frame\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, 0] = _node_dof(int(ni + 1), 0)\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.1\nwarped = g.copy()\nwarped.points = np.asarray(g.points) + scale * disp_xyz\n\nplotter = pv.Plotter(off_screen=True, window_size=(600, 600))\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([[0.0, L_v, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=14,\n    color=\"#888888\",\n    label=\"rigid corner\",\n)\nplotter.add_points(\n    np.array([[L_h, L_v + 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):.4f} % of P L_h\u00b2 L_v / E I + P L_h\u00b3 / (3 E I) + P L_v / (E A).\")\nprint(\n    \"  \u2022 Bending dominates the slender-frame response; the column's axial \"\n    f\"compression contributes only {abs(delta_axial / delta_bending) * 100:.4f} % \"\n    \"of the tip deflection at this geometry.\"\n)\nprint(\n    \"  \u2022 Limit cases collapse cleanly: L_v \u2192 0 \u21d2 horizontal cantilever P L\u00b3/(3 E I); \"\n    \"L_h \u2192 0 \u21d2 axial column P L / (E A).\"\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
}