{
  "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# Plate with a circular hole \u2014 Kirsch stress concentration\n\nClassical 2D elasticity benchmark: a thin rectangular plate\nwith a small circular hole, loaded by uniform remote tension\n$\\sigma_{\\infty}$ along the $x$-axis.  Kirsch's\n1898 closed form for the **infinite plate** gives a stress\nfield that, on the hole edge $r = a$, peaks at the\npoints perpendicular to the loading direction\n($\\theta = \\pm\\pi/2$):\n\n\\begin{align}:label: kirsch-kt\n\n    K_{t} \\;\\equiv\\; \\frac{\\sigma_{yy}^{(\\max)}}{\\sigma_{\\infty}}\n    \\;=\\; 3\n    \\qquad \\text{(infinite plate, hole radius } a \\text{)}.\\end{align}\n\nThis is the textbook **3:1 stress concentration**, the\nfoundational example for fatigue-design notch factors and\nthe canonical first-principles check for stress-recovery\ncorrectness.\n\n## Implementation\n\nQuarter-symmetry plane-stress model\n(:class:`~femorph_solver.validation.problems.PlateWithHole`)\non a $16 \\times 12$ QUAD4 mapped mesh that radially\nclusters nodes near the hole where the stress gradient is\nsteep.  Symmetry pins $u_y = 0$ on $y = 0$ and\n$u_x = 0$ on $x = 0$; the right edge\n($x = W/2$) carries a consistent traction load\nequivalent to a uniform $\\sigma_{\\infty}$ remote tension.\nFinite-width correction (Howland 1930; Peterson Table 4-1)\ngives $K_t \\approx 3.04$ at $a/W = 0.1$, so the\nexpected FE result lies between the Kirsch infinite-plate\nvalue (3.0) and the finite-width correction (~3.04) plus a\nsmall mesh-discretisation overhead.\n\n## References\n\n* Kirsch, G. (1898) \"Die Theorie der Elastizit\u00e4t und die\n  Bed\u00fcrfnisse der Festigkeitslehre,\" *Zeitschrift des\n  Vereins Deutscher Ingenieure* 42, 797\u2013807 \u2014 the original\n  closed-form derivation.\n* Timoshenko, S. P. and Goodier, J. N. (1970) *Theory of\n  Elasticity*, 3rd ed., McGraw-Hill, \u00a735 (plate with a\n  circular hole).\n* Howland, R. C. J. (1930) \"On the stresses in the\n  neighbourhood of a circular hole in a strip under\n  tension,\" *Phil. Trans. R. Soc. A* 229, 49\u201386 \u2014 finite-\n  width correction.\n* Peterson, R. E. (1974) *Stress Concentration Factors*,\n  Wiley, Table 4-1 (uniaxial stress, hole in finite plate).\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, \u00a75.5 (stress concentration\n  benchmark).\n\n## Vendor cross-references\n\n======================================  =========================  ============================================\nSource                                   \u03c3_xx at hole top [MPa]     Problem ID / location\n======================================  =========================  ============================================\nClosed form (Kirsch 1898)                30.00                      K_t = 3 infinite plate\nTimoshenko & Goodier (1970)              30.00                      Theory of Elasticity \u00a735\nPeterson (2008) \u2014 finite width           \u2248 30.4                     K_t(a/W=0.1) \u2248 3.04\nMAPDL Verification Manual                \u2248 30.0                     VM-74 Stress concentration around a hole\nAbaqus Verification Manual               \u2248 30.0                     AVM 1.3.6 plate-with-hole family\nNAFEMS                                   \u2248 30.0                     NL2 plate-with-hole (linear-elastic limit)\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 PlateWithHole"
      ]
    },
    {
      "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 = PlateWithHole()\nm = problem.build_model()\n\nprint(\n    f\"plate-with-hole mesh: {m.grid.n_points} nodes, {m.grid.n_cells} QUAD4 cells; \"\n    f\"a = {problem.a} m, W = {problem.W} m, \u03c3_\u221e = {problem.sigma_infty / 1e6:.1f} MPa\"\n)\nprint(f\"E = {problem.E / 1e9:.0f} GPa, \u03bd = {problem.nu}, a/W = {problem.a / problem.W:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Static solve\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "res = m.solve_static()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Extract \u03c3_yy at the hole edge \u2014 the Kirsch peak\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sigma_yy_max = problem.extract(m, res, \"kt_at_hole\")\nsigma_yy_kirsch = 3.0 * problem.sigma_infty\nkt_computed = sigma_yy_max / problem.sigma_infty\nerr_pct = (sigma_yy_max - sigma_yy_kirsch) / sigma_yy_kirsch * 100.0\n\nprint()\nprint(f\"\u03c3_yy at hole edge   computed = {sigma_yy_max / 1e6:8.3f} MPa\")\nprint(f\"\u03c3_yy at hole edge   Kirsch   = {sigma_yy_kirsch / 1e6:8.3f} MPa  (3 \u00b7 \u03c3_\u221e)\")\nprint(f\"K_t computed                 = {kt_computed:8.4f}\")\nprint(\"K_t Kirsch (infinite plate)  = 3.0000\")\nprint(\"K_t Peterson (a/W = 0.1)     \u2248 3.04 (Howland 1930 finite-width correction)\")\nprint(f\"relative error               = {err_pct:+.2f} %\")\n\n# 10% tolerance \u2014 Kirsch's reference is for an infinite plate; the\n# finite-width Howland correction adds ~1 % at a/W = 0.1, and a\n# QUAD4 mesh at 16 \u00d7 12 hits the gradient with ~3-4 % discretisation\n# error on top.\nassert abs(err_pct) < 10.0, f\"K_t deviation {err_pct:.2f}% exceeds 10%\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the \u03c3_yy field on the deformed mesh\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, 2)\ndisp3 = np.zeros((grid.n_points, 3), dtype=np.float64)\ndisp3[:, :2] = u\ngrid.point_data[\"displacement\"] = disp3\ngrid.point_data[\"|u|\"] = np.linalg.norm(u, axis=1)\n\n# \u03c3_yy approximation across the mesh: derive from the displacement\n# field via a finite-difference of u_y in y at every interior node.\n# Adequate for a contour render; the validation problem's `extract`\n# method computes \u03c3_yy(hole edge) more carefully via the plane-stress\n# C-matrix at a single key point.\npts = np.asarray(grid.points)\nsigma_yy_field = np.zeros(grid.n_points, dtype=np.float64)\n# Use the same scaling factor from the analytical field at the hole\n# edge so the colour bar's range covers [0, K_t \u00b7 \u03c3_\u221e] cleanly.\nsigma_yy_field[:] = problem.sigma_infty  # stays at \u03c3_\u221e far from the hole\n# Spike \u03c3_yy near the hole edge (r \u2248 a) to the analytical\n# 3\u00b7\u03c3_\u221e value at \u03b8 = \u03c0/2.  This is for visualisation only \u2014 the\n# rigorous value comes from problem.extract.\nr = np.linalg.norm(pts[:, :2], axis=1)\ntheta = np.arctan2(pts[:, 1], pts[:, 0])\n# Kirsch closed-form \u03c3_yy(r, \u03b8) on the y-axis (cos 2\u03b8 = -1 at \u03b8 = \u03c0/2):\nnear_hole = r > 0\nratio = problem.a**2 / r[near_hole] ** 2\nratio4 = problem.a**4 / r[near_hole] ** 4\nsigma_yy_field[near_hole] = problem.sigma_infty * (\n    1.0 + 0.5 * ratio + (1.5 * ratio4 - 0.5 * ratio) * np.cos(2.0 * theta[near_hole])\n)\ngrid.point_data[\"\u03c3_yy [Pa]\"] = sigma_yy_field\n\nwarped = grid.warp_by_vector(\"displacement\", factor=200.0)\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 480))\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    warped,\n    scalars=\"\u03c3_yy [Pa]\",\n    cmap=\"plasma\",\n    show_edges=True,\n    scalar_bar_args={\"title\": \"\u03c3_yy [Pa]  (deformed \u00d7200, Kirsch closed form)\"},\n)\n# Mark the Kirsch peak point at (0, a) \u2014 the locus of K_t = 3.\nplotter.add_points(\n    np.array([[0.0, problem.a, 0.0]]),\n    render_points_as_spheres=True,\n    point_size=18,\n    color=\"#d62728\",\n    label=f\"K_t = {kt_computed:.3f} at (0, a)\",\n)\nplotter.view_xy()\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
}