{
  "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# Shear-locking demonstration \u2014 HEX8 integration variants\n\nThe textbook trap of low-order solid elements: a fully-\nintegrated 8-node trilinear hex (HEX8 with full 2 \u00d7 2 \u00d7 2\nGauss, **no** B-bar / enhanced strain) **cannot** reproduce\nlinear-bending kinematics, because the constant-strain\nGauss point at the element centre overstates the bending\nstiffness (Cook \u00a76.6, \u00a76.13).  The result is **shear locking**:\neven on a slender cantilever where the Euler-Bernoulli\nanalytical answer is exact, the FE displacement is too small \u2014\nsometimes by orders of magnitude \u2014 and the missing strain\nenergy goes into spurious transverse shear at the Gauss point.\n\nTwo practitioner cures, both shipped with femorph-solver's\n:doc:`/reference/elements/hex8`:\n\n* **B-bar volumetric stabilisation** (Hughes 1980, default\n  ``ELEMENTS.HEX8``): averages the volumetric component of the\n  strain field across the element, suppressing the related\n  *volumetric* locking at $\\nu \\to 0.5$ but not the\n  shear locking on its own.\n\n* **Enhanced assumed strain** (EAS \u2014 Simo\u2013Rifai 1990,\n  ``ELEMENTS.HEX8(integration=\"enhanced_strain\")``): adds nine\n  static-condensed internal modes that fix bending-mode shear\n  locking exactly.  The natural choice for thin-bending solid\n  geometries.\n\nThis example solves the same cantilever-tip-load problem with\n**three integration variants** side by side and prints the tip-\ndeflection ratio against the Euler-Bernoulli closed form\n$\\delta = P L^{3} / (3 E I)$:\n\n* ``ELEMENTS.HEX8(integration=\"plain_gauss\")`` \u2014 full 2\u00d72\u00d72\n  Gauss, no stabilisation (the textbook locking case).\n* ``ELEMENTS.HEX8`` \u2014 full Gauss + B-bar (default).\n* ``ELEMENTS.HEX8(integration=\"enhanced_strain\")`` \u2014 EAS.\n\n## References\n\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, \u00a76.6 (volumetric locking),\n  \u00a76.13 (shear locking + EAS), \u00a717.2 (low-order shell\n  cure).\n* Hughes, T. J. R. (1980) \"Generalization of selective\n  integration procedures to anisotropic and nonlinear media,\"\n  *International Journal for Numerical Methods in\n  Engineering* 15 (9), 1413\u20131418 (B-bar).\n* Simo, J. C. and Rifai, M. S. (1990) \"A class of mixed\n  assumed strain methods and the method of incompatible\n  modes,\" *International Journal for Numerical Methods in\n  Engineering* 29 (8), 1595\u20131638 (EAS).\n* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,\n  \u00a75.3.4 (incompatible modes / EAS).\n* Belytschko, T., Liu, W. K., Moran, B. and Elkhodary,\n  K. (2014) *Nonlinear Finite Elements for Continua and\n  Structures*, 2nd ed., Wiley, \u00a78.4 (locking modes).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport matplotlib.pyplot as plt\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.3\nRHO = 7850.0\nb = 0.01  # cross-section side; L / b = 100 \u2192 very slender \u2192 severe locking\nA = b * b\nI = b**4 / 12.0  # noqa: E741\nL = 1.0\nP_TOTAL = -100.0  # N (small load \u2192 linear regime)\n\nNX, NY, NZ = 20, 1, 1  # 20 elements axial, 1 through each thickness\n\ndelta_eb = P_TOTAL * L**3 / (3.0 * E * I)\nprint(f\"L / b = {L / b:.0f} (slender)\")\nprint(f\"\u03b4_EB = P L^3 / (3 EI) = {delta_eb * 1e6:.3f} \u00b5m\")\n\n\ndef run_with_integration(label: str, kwargs: dict) -> tuple[float, float]:\n    \"\"\"Solve the same cantilever with the given HEX8 spec; return\n    (UY tip mid-surface, ratio_to_EB).\"\"\"\n    xs = np.linspace(0.0, L, NX + 1)\n    ys = np.linspace(0.0, b, NY + 1)\n    zs = np.linspace(0.0, b, NZ + 1)\n    grid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\n\n    m = femorph_solver.Model.from_grid(grid)\n    if kwargs:\n        m.assign(ELEMENTS.HEX8(**kwargs), material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})\n    else:\n        m.assign(ELEMENTS.HEX8, material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})\n\n    pts = np.asarray(m.grid.points)\n    node_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"])\n    fixed = node_nums[pts[:, 0] < 1e-9]\n    m.fix(nodes=fixed.tolist(), dof=\"ALL\")\n    tip_nodes = node_nums[pts[:, 0] > L - 1e-9]\n    fz_per = P_TOTAL / len(tip_nodes)\n    for nn in tip_nodes:\n        m.apply_force(int(nn), fy=fz_per)\n\n    res = m.solve_static()\n    g = femorph_solver.io.static_result_to_grid(m, res)\n    pts_g = np.asarray(g.points)\n    tip_mask = pts_g[:, 0] > L - 1e-9\n    uy_tip = float(g.point_data[\"displacement\"][tip_mask, 1].mean())\n    return uy_tip, uy_tip / delta_eb"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Three integration variants \u2014 side by side\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "variants = [\n    (\"plain_gauss\", {\"integration\": \"plain_gauss\"}, \"#d62728\"),\n    (\"default (B-bar)\", {}, \"#1f77b4\"),\n    (\"enhanced_strain (EAS)\", {\"integration\": \"enhanced_strain\"}, \"#2ca02c\"),\n]\nresults: list[tuple[str, float, float, str]] = []\n\nprint()\nprint(f\"{'integration':>26}  {'UY_tip [\u00b5m]':>12}  {'UY / \u03b4_EB':>10}\")\nprint(f\"{'-' * 26:>26}  {'-' * 12:>12}  {'-' * 10:>10}\")\nfor label, kwargs, colour in variants:\n    uy, ratio = run_with_integration(label, kwargs)\n    results.append((label, uy, ratio, colour))\n    print(f\"{label:>26}  {uy * 1e6:12.4f}  {ratio:10.6f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Verify the textbook ordering\n\nThe locked plain-Gauss solution must underpredict the\nEuler-Bernoulli reference on a slender beam (textbook\nexpectation: locks to ~10\u201330 % of the analytical answer at\nL/b = 100, NX = 20 with one element through the thickness).\nB-bar is **not** a shear-locking cure \u2014 it cures volumetric\nlocking only \u2014 so it should give the same locked deflection\nas plain Gauss in this benchmark.  Only EAS recovers the\nbending mode.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ratios = {label: ratio for label, _, ratio, _ in results}\nassert ratios[\"plain_gauss\"] < 0.5, \"plain_gauss should lock to << 1.0\"\nassert ratios[\"default (B-bar)\"] < 0.5, \"B-bar alone should also lock here\"\nassert ratios[\"enhanced_strain (EAS)\"] > 0.85, \"EAS should recover \u2265 85 % of EB\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the deflection-ratio comparison\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "labels = [r[0] for r in results]\nratios_values = [r[2] for r in results]\ncolours = [r[3] for r in results]\n\nfig, ax = plt.subplots(1, 1, figsize=(7.0, 4.0))\nbars = ax.bar(labels, ratios_values, color=colours)\nax.axhline(1.0, linestyle=\"--\", color=\"black\", lw=1.5, label=\"Euler-Bernoulli (target)\")\nax.set_ylabel(r\"$\\delta^{h}_{\\mathrm{tip}} / \\delta_{\\mathrm{EB}}$\")\nax.set_title(\n    f\"Cantilever HEX8 shear-locking \u2014 L / b = {L / b:.0f}, N_x = {NX}, 1 element through h\",\n    fontsize=10,\n)\nax.set_ylim(0.0, 1.10)\nfor bar, ratio in zip(bars, ratios_values):\n    ax.annotate(\n        f\"{ratio:.3f}\",\n        xy=(bar.get_x() + bar.get_width() / 2, ratio),\n        xytext=(0, 4),\n        textcoords=\"offset points\",\n        ha=\"center\",\n        fontsize=9,\n    )\nax.legend(loc=\"upper left\", fontsize=9)\nax.grid(True, axis=\"y\", ls=\":\", alpha=0.5)\nfig.tight_layout()\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Closing note\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print()\nprint(\"Take-aways:\")\nprint(\n    f\"  \u2022 plain_gauss    \u2192 {ratios['plain_gauss']:.1%} of EB  (locked; spurious shear-strain energy)\"\n)\nprint(\n    f\"  \u2022 B-bar default  \u2192 {ratios['default (B-bar)']:.1%} of EB  (B-bar cures *volumetric* locking, not shear)\"\n)\nprint(\n    f\"  \u2022 EAS variant    \u2192 {ratios['enhanced_strain (EAS)']:.1%} of EB  (the right cure for thin-bending solid kernels)\"\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
}