{ "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# Cook's membrane \u2014 mesh-distortion benchmark (1974)\n\nRobert D. Cook's 1974 *International Journal for Numerical Methods\nin Engineering* paper introduced what became the universal stress-\ntest for low-order solid elements under **mesh distortion**: a\ntapered plate clamped on its short edge and loaded by a uniform\ntransverse shear on its slanted edge. The non-rectangular geometry\nforces every interior element to be a parallelogram (or worse) \u2014\nthe canonical setting where a fully-integrated 8-node hex (plain\n2 \u00d7 2 \u00d7 2 Gauss) **shear-locks** badly even though the material\nis moderately compressible ($\\nu = 1/3$) and no volumetric\npathology is present.\n\nThe reference quantity is the **vertical displacement of the upper-\nright corner** $C = (48, 60)$. Successive refinements of\nQUAD4 / HEX8 with progressively better formulations track the\nmesh-converged value $u_y(C) \\approx 23.96$ (Cook 1989 \u00a76.13;\nBelytschko et al. 2014 \u00a78.4 Table 8.1):\n\n* **Plain 2 \u00d7 2 \u00d7 2 Gauss** locks badly even at $8 \\times 8$,\n reporting $\\approx 21$ \u2014 about 88 % of the converged value.\n* **B-bar** (selective dilatational integration, Hughes 1980) helps\n marginally on this problem because the lock is **shear**, not\n volumetric \u2014 recovers $\\approx 22$.\n* **Enhanced assumed strain** (EAS, Simo-Rifai 1990) is the\n textbook cure: recovers $\\approx 23.6$-$23.9$ on the\n same coarse mesh.\n\nThis is the same hierarchy\n`sphx_glr_gallery_verification_example_verify_shear_locking_demo.py`\nshows on a slender prismatic cantilever \u2014 the present example\nextends that lesson to a **distorted-mesh** scenario, which is\nwhere production meshes actually live.\n\n## Reference geometry (Cook 1974 Fig. 1)\n\nTrapezoidal plate, plane-stress assumption (we model it as a\nthin 3D slab with the through-thickness $u_z$ free):\n\n```text\nC = (48, 60)\n```\n \u2571\u2502\n \u2571 \u2502\n D \u2500\u2500\u2500\u2500\u2500\u2571 \u2502\n (0,44) \u2502 \u2502\n \u2502 \u2502 \u2190 uniform vertical shear\n \u2502 \u2502 F_total = 1\n \u2502 \u2502\n A \u2502 \u2502\n (0, 0)\u2500\u2534\u2500\u2500B = (48, 44)\n\n* Clamped on edge $AD$ ($x = 0$).\n* Vertical shear traction on edge $BC$ ($x = 48$),\n total force $F = 1$ (unit).\n* Material: $E = 1$, $\\nu = 1/3$.\n\n## References\n\n* Cook, R. D. (1974) \"Improved two-dimensional finite element,\"\n *Journal of the Structural Division (ASCE)* 100 (ST9),\n 1851\u20131863 \u2014 the original benchmark introduction.\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.13 (selective-reduced\n integration), \u00a76.14 (incompatible modes).\n* Belytschko, T., Liu, W. K., Moran, B. and Elkhodary, K. (2014)\n *Nonlinear Finite Elements for Continua and Structures*, 2nd\n ed., Wiley, \u00a78.4 Table 8.1 (Cook's-membrane reference values).\n* Simo, J. C. and Rifai, M. S. (1990) \"A class of mixed assumed\n strain methods and the method of incompatible modes,\"\n *International Journal for Numerical Methods in Engineering*\n 29 (8), 1595\u20131638 \u2014 the EAS formulation.\n* Hughes, T. J. R. (1980) \"Generalization of selective integration\n procedures to anisotropic and nonlinear media,\"\n *International Journal for Numerical Methods in Engineering*\n 15 (9), 1413\u20131418 \u2014 B-bar.\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": [ "## Reference geometry & material \u2014 Cook 1974\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "E = 1.0\nNU = 1.0 / 3.0\nRHO = 1.0 # irrelevant for the static solve\nTHICKNESS = 1.0\nF_TOTAL = 1.0 # total transverse force on the loaded edge\n\n# Cook 1974 reference value for the corner UY: the mesh-converged\n# QUAD4-EAS *plane-stress* limit reported in Belytschko et al. 2014\n# Table 8.1 is 23.96. Earlier sources (Cook 1989 \u00a76.13) round to\n# 23.91. We model the membrane as a thin 3D slab with the through-\n# thickness u_z free (a 3D analogue of plane stress), so the HEX8-\n# EAS converged answer overshoots the 2D-Q4 reference slightly\n# (~+3 %) \u2014 the membrane's through-thickness freedom adds a Poisson\n# contribution that 2D plane stress collapses out.\nUY_C_REF = 23.96\n\nprint(\"Cook's membrane \u2014 Cook 1974 mesh-distortion benchmark\")\nprint(\" geometry vertices A=(0,0), B=(48,44), C=(48,60), D=(0,44)\")\nprint(f\" E = {E}, \u03bd = {NU}, t = {THICKNESS}, F_total = {F_TOTAL}\")\nprint(f\" reference u_y(C) (Belytschko 2014 Table 8.1) = {UY_C_REF:.3f}\")\n\n\ndef build_grid(n: int) -> pv.UnstructuredGrid:\n \"\"\"Bilinearly-mapped n \u00d7 n \u00d7 1 hex slab for Cook's membrane.\"\"\"\n A = np.array([0.0, 0.0])\n B = np.array([48.0, 44.0])\n C = np.array([48.0, 60.0])\n D = np.array([0.0, 44.0])\n\n # Bilinear xi-eta to xy mapping over the trapezoid. \u03be runs A\u2192B\n # along the bottom edge, \u03b7 runs A\u2192D along the left edge.\n xi = np.linspace(0.0, 1.0, n + 1)\n eta = np.linspace(0.0, 1.0, n + 1)\n pts: list[list[float]] = []\n for z in (0.0, THICKNESS):\n for j in range(n + 1):\n for i in range(n + 1):\n # Bilinear blend of the four corners.\n p = (\n (1 - xi[i]) * (1 - eta[j]) * A\n + xi[i] * (1 - eta[j]) * B\n + xi[i] * eta[j] * C\n + (1 - xi[i]) * eta[j] * D\n )\n pts.append([p[0], p[1], z])\n pts_arr = np.array(pts, dtype=np.float64)\n\n nx_plane = (n + 1) * (n + 1)\n n_cells = n * n\n cells = np.empty((n_cells, 9), dtype=np.int64)\n cells[:, 0] = 8\n c = 0\n for j in range(n):\n for i in range(n):\n n00b = j * (n + 1) + i\n n10b = j * (n + 1) + (i + 1)\n n11b = (j + 1) * (n + 1) + (i + 1)\n n01b = (j + 1) * (n + 1) + i\n cells[c, 1:] = (\n n00b,\n n10b,\n n11b,\n n01b,\n n00b + nx_plane,\n n10b + nx_plane,\n n11b + nx_plane,\n n01b + nx_plane,\n )\n c += 1\n return pv.UnstructuredGrid(\n cells.ravel(),\n np.full(n_cells, 12, dtype=np.uint8),\n pts_arr,\n )\n\n\ndef run_one(n: int, integration: str | None) -> float:\n \"\"\"Solve Cook's membrane on an n \u00d7 n hex slab. Returns u_y(C).\"\"\"\n grid = build_grid(n)\n m = femorph_solver.Model.from_grid(grid)\n spec = ELEMENTS.HEX8(integration=integration) if integration else ELEMENTS.HEX8\n m.assign(spec, 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 tol = 1e-9\n\n # Clamp left edge AD (x = 0) \u2014 pin all three displacement DOFs.\n left_mask = pts[:, 0] < tol\n for nn in node_nums[left_mask]:\n m.fix(nodes=[int(nn)], dof=\"ALL\")\n\n # Plane stress is approximated by a thin slab with the through-\n # thickness u_z left free; nothing extra to do.\n\n # Distributed transverse shear on the right edge BC (x = 48).\n # Total force F_TOTAL spread evenly over the right-edge nodes\n # at both z = 0 and z = THICKNESS.\n right_mask = np.abs(pts[:, 0] - 48.0) < tol\n right_nodes = node_nums[right_mask]\n fy_per = F_TOTAL / len(right_nodes)\n for nn in right_nodes:\n m.apply_force(int(nn), fy=fy_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 # Reference corner C = (48, 60) at any z; pick the closest node.\n target = np.array([48.0, 60.0])\n dist = np.linalg.norm(pts_g[:, :2] - target, axis=1)\n c_idx = int(np.argmin(dist))\n return float(g.point_data[\"displacement\"][c_idx, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run all three integration variants on a single coarse mesh\n\nWe pick $n = 8$ (8 \u00d7 8 \u00d7 1 mesh, 162 nodes total). Cook's\noriginal 1974 paper runs at $2 \\times 2$, $4 \\times 4$,\nand $8 \\times 8$; the third is enough to expose the locking\nhierarchy without burning too much CPU.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N_MESH = 8\n\nprint()\nprint(f\"On the {N_MESH} \u00d7 {N_MESH} \u00d7 1 mesh ({(N_MESH + 1) ** 2 * 2} nodes):\")\nprint()\nprint(f\" {'integration':<22} {'u_y(C)':>10} {'u_y / u_ref':>12} {'verdict':<28}\")\nprint(f\" {'-' * 22:<22} {'-' * 10:>10} {'-' * 12:>12} {'-' * 28:<28}\")\n\nvariants = (\n (\"plain_gauss\", \"shear-locks badly\"),\n (None, \"default B-bar \u2014 partial cure\"),\n (\"enhanced_strain\", \"EAS \u2014 fully converged\"),\n)\n\nratios: dict[str, float] = {}\nuy_results: dict[str, float] = {}\nfor integ, verdict in variants:\n uy_c = run_one(N_MESH, integ)\n ratio = uy_c / UY_C_REF\n label = integ if integ else \"default (B-bar)\"\n ratios[label] = ratio\n uy_results[label] = uy_c\n print(f\" {label:<22} {uy_c:>10.4f} {ratio:>12.4f} {verdict:<28}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verify the textbook ordering\n\nPlain Gauss is known to recover only ~85-92 % on an 8 \u00d7 8 mesh\n(Cook 1989 Table 6.13.1 / Belytschko 2014 \u00a78.4); B-bar gives a\nsmall bump to ~92-94 % because the lock here is *shear*, not\nvolumetric. EAS lifts the result to \u2265 98 % (and slightly\novershoots the 2D plane-stress Q4 reference because the 3D slab\nis plane-stress-like, not pure plane stress).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "assert ratios[\"plain_gauss\"] < 0.94, (\n f\"plain Gauss must lock on Cook's membrane \u2014 got {ratios['plain_gauss']:.4f}\"\n)\nassert ratios[\"enhanced_strain\"] > 0.98, (\n f\"EAS must recover \u2265 98 % \u2014 got {ratios['enhanced_strain']:.4f}\"\n)\nassert ratios[\"enhanced_strain\"] > ratios[\"plain_gauss\"], \"EAS must beat plain Gauss\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Refinement ladder \u2014 show EAS converges fastest\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ladder = (2, 4, 8, 16)\nplain: list[float] = []\nbbar: list[float] = []\neas: list[float] = []\nprint()\nprint(\"Refinement ladder:\")\nprint(f\" {'n':>4} {'plain_gauss':>12} {'B-bar':>10} {'EAS':>10}\")\nprint(f\" {'-' * 4:>4} {'-' * 12:>12} {'-' * 10:>10} {'-' * 10:>10}\")\nfor n in ladder:\n pg = run_one(n, \"plain_gauss\")\n bb = run_one(n, None)\n en = run_one(n, \"enhanced_strain\")\n plain.append(pg / UY_C_REF)\n bbar.append(bb / UY_C_REF)\n eas.append(en / UY_C_REF)\n print(f\" {n:>4} {pg:>12.4f} {bb:>10.4f} {en:>10.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the convergence comparison\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(7.0, 4.0))\nax.plot(ladder, plain, \"o-\", color=\"#d62728\", lw=2, label=\"plain_gauss (shear-locks)\")\nax.plot(ladder, bbar, \"s-\", color=\"#9467bd\", lw=2, label=\"default (B-bar)\")\nax.plot(ladder, eas, \"^-\", color=\"#1f77b4\", lw=2, label=\"enhanced_strain (EAS)\")\nax.axhline(1.0, color=\"black\", ls=\"--\", lw=1.0, label=\"mesh-converged limit\")\nax.set_xlabel(\"mesh refinement n (n \u00d7 n \u00d7 1)\")\nax.set_ylabel(r\"$u_y(C)\\, /\\, u_y^{\\mathrm{ref}}$\")\nax.set_xscale(\"log\", base=2)\nax.set_xticks(ladder)\nax.set_xticklabels([str(n) for n in ladder])\nax.set_title(\"Cook's membrane \u2014 corner displacement convergence\", fontsize=11)\nax.set_ylim(0.55, 1.05)\nax.legend(loc=\"lower right\", fontsize=9)\nax.grid(True, ls=\":\", alpha=0.5)\nfig.tight_layout()\nfig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the deformed shape (EAS, n=8)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "grid = build_grid(N_MESH)\nm_eas = femorph_solver.Model.from_grid(grid)\nm_eas.assign(\n ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n)\npts = np.asarray(m_eas.grid.points)\nnode_nums = np.asarray(m_eas.grid.point_data[\"ansys_node_num\"])\ntol = 1e-9\nfor nn in node_nums[pts[:, 0] < tol]:\n m_eas.fix(nodes=[int(nn)], dof=\"ALL\")\nright = node_nums[np.abs(pts[:, 0] - 48.0) < tol]\nfor nn in right:\n m_eas.apply_force(int(nn), fy=F_TOTAL / len(right))\nres_eas = m_eas.solve_static()\ng_eas = femorph_solver.io.static_result_to_grid(m_eas, res_eas)\n\n# Magnify the displacement to make the deformed shape readable.\nwarp_factor = 0.5 # a u_y(C) of ~24 over an L=48 plate is already\n# visible at unit scale; bump to 0.5\u00d7 for legibility.\nwarped = g_eas.copy()\nwarped.points = np.asarray(g_eas.points) + warp_factor * np.asarray(\n g_eas.point_data[\"displacement\"]\n)\nwarped[\"|u|\"] = np.linalg.norm(g_eas.point_data[\"displacement\"], axis=1)\n\nplotter = pv.Plotter(off_screen=True, window_size=(720, 540))\nplotter.add_mesh(g_eas, style=\"wireframe\", color=\"grey\", opacity=0.4, line_width=1)\nplotter.add_mesh(\n warped,\n scalars=\"|u|\",\n cmap=\"viridis\",\n show_edges=True,\n scalar_bar_args={\"title\": \"|u| (\u00d70.5 warp, EAS)\"},\n)\nplotter.add_points(\n np.array([[48.0, 60.0, 0.5]]),\n render_points_as_spheres=True,\n point_size=18,\n color=\"#d62728\",\n label=f\"corner C \u2014 u_y = {uy_results['enhanced_strain']:.3f}\",\n)\nplotter.add_legend()\nplotter.view_xy()\nplotter.camera.zoom(1.05)\nplotter.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(f\" \u2022 plain_gauss at n = 8 \u2192 {ratios['plain_gauss']:.1%} of converged (shear-locked)\")\nprint(f\" \u2022 default B-bar at n = 8 \u2192 {ratios['default (B-bar)']:.1%} (only marginal help)\")\nprint(f\" \u2022 enhanced_strain at n = 8 \u2192 {ratios['enhanced_strain']:.1%} (EAS is the cure)\")\nprint(\" \u2022 EAS's edge widens with mesh distortion \u2014 pick it for any production hex mesh.\")" ] } ], "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 }