r""" Cook's membrane — mesh-distortion benchmark (1974) ================================================== Robert D. Cook's 1974 *International Journal for Numerical Methods in Engineering* paper introduced what became the universal stress- test for low-order solid elements under **mesh distortion**: a tapered plate clamped on its short edge and loaded by a uniform transverse shear on its slanted edge. The non-rectangular geometry forces every interior element to be a parallelogram (or worse) — the canonical setting where a fully-integrated 8-node hex (plain 2 × 2 × 2 Gauss) **shear-locks** badly even though the material is moderately compressible (:math:`\nu = 1/3`) and no volumetric pathology is present. The reference quantity is the **vertical displacement of the upper- right corner** :math:`C = (48, 60)`. Successive refinements of QUAD4 / HEX8 with progressively better formulations track the mesh-converged value :math:`u_y(C) \approx 23.96` (Cook 1989 §6.13; Belytschko et al. 2014 §8.4 Table 8.1): * **Plain 2 × 2 × 2 Gauss** locks badly even at :math:`8 \times 8`, reporting :math:`\approx 21` — about 88 % of the converged value. * **B-bar** (selective dilatational integration, Hughes 1980) helps marginally on this problem because the lock is **shear**, not volumetric — recovers :math:`\approx 22`. * **Enhanced assumed strain** (EAS, Simo-Rifai 1990) is the textbook cure: recovers :math:`\approx 23.6`-:math:`23.9` on the same coarse mesh. This is the same hierarchy :ref:`sphx_glr_gallery_verification_example_verify_shear_locking_demo.py` shows on a slender prismatic cantilever — the present example extends that lesson to a **distorted-mesh** scenario, which is where production meshes actually live. Reference geometry (Cook 1974 Fig. 1) ------------------------------------- Trapezoidal plate, plane-stress assumption (we model it as a thin 3D slab with the through-thickness :math:`u_z` free): .. code-block:: text C = (48, 60) ╱│ ╱ │ D ─────╱ │ (0,44) │ │ │ │ ← uniform vertical shear │ │ F_total = 1 │ │ A │ │ (0, 0)─┴──B = (48, 44) * Clamped on edge :math:`AD` (:math:`x = 0`). * Vertical shear traction on edge :math:`BC` (:math:`x = 48`), total force :math:`F = 1` (unit). * Material: :math:`E = 1`, :math:`\nu = 1/3`. References ---------- * Cook, R. D. (1974) "Improved two-dimensional finite element," *Journal of the Structural Division (ASCE)* 100 (ST9), 1851–1863 — the original benchmark introduction. * Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley, §6.13 (selective-reduced integration), §6.14 (incompatible modes). * Belytschko, T., Liu, W. K., Moran, B. and Elkhodary, K. (2014) *Nonlinear Finite Elements for Continua and Structures*, 2nd ed., Wiley, §8.4 Table 8.1 (Cook's-membrane reference values). * Simo, J. C. and Rifai, M. S. (1990) "A class of mixed assumed strain methods and the method of incompatible modes," *International Journal for Numerical Methods in Engineering* 29 (8), 1595–1638 — the EAS formulation. * Hughes, T. J. R. (1980) "Generalization of selective integration procedures to anisotropic and nonlinear media," *International Journal for Numerical Methods in Engineering* 15 (9), 1413–1418 — B-bar. """ from __future__ import annotations import matplotlib.pyplot as plt import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS # %% # Reference geometry & material — Cook 1974 # ----------------------------------------- E = 1.0 NU = 1.0 / 3.0 RHO = 1.0 # irrelevant for the static solve THICKNESS = 1.0 F_TOTAL = 1.0 # total transverse force on the loaded edge # Cook 1974 reference value for the corner UY: the mesh-converged # QUAD4-EAS *plane-stress* limit reported in Belytschko et al. 2014 # Table 8.1 is 23.96. Earlier sources (Cook 1989 §6.13) round to # 23.91. We model the membrane as a thin 3D slab with the through- # thickness u_z free (a 3D analogue of plane stress), so the HEX8- # EAS converged answer overshoots the 2D-Q4 reference slightly # (~+3 %) — the membrane's through-thickness freedom adds a Poisson # contribution that 2D plane stress collapses out. UY_C_REF = 23.96 print("Cook's membrane — Cook 1974 mesh-distortion benchmark") print(" geometry vertices A=(0,0), B=(48,44), C=(48,60), D=(0,44)") print(f" E = {E}, ν = {NU}, t = {THICKNESS}, F_total = {F_TOTAL}") print(f" reference u_y(C) (Belytschko 2014 Table 8.1) = {UY_C_REF:.3f}") def build_grid(n: int) -> pv.UnstructuredGrid: """Bilinearly-mapped n × n × 1 hex slab for Cook's membrane.""" A = np.array([0.0, 0.0]) B = np.array([48.0, 44.0]) C = np.array([48.0, 60.0]) D = np.array([0.0, 44.0]) # Bilinear xi-eta to xy mapping over the trapezoid. ξ runs A→B # along the bottom edge, η runs A→D along the left edge. xi = np.linspace(0.0, 1.0, n + 1) eta = np.linspace(0.0, 1.0, n + 1) pts: list[list[float]] = [] for z in (0.0, THICKNESS): for j in range(n + 1): for i in range(n + 1): # Bilinear blend of the four corners. p = ( (1 - xi[i]) * (1 - eta[j]) * A + xi[i] * (1 - eta[j]) * B + xi[i] * eta[j] * C + (1 - xi[i]) * eta[j] * D ) pts.append([p[0], p[1], z]) pts_arr = np.array(pts, dtype=np.float64) nx_plane = (n + 1) * (n + 1) n_cells = n * n cells = np.empty((n_cells, 9), dtype=np.int64) cells[:, 0] = 8 c = 0 for j in range(n): for i in range(n): n00b = j * (n + 1) + i n10b = j * (n + 1) + (i + 1) n11b = (j + 1) * (n + 1) + (i + 1) n01b = (j + 1) * (n + 1) + i cells[c, 1:] = ( n00b, n10b, n11b, n01b, n00b + nx_plane, n10b + nx_plane, n11b + nx_plane, n01b + nx_plane, ) c += 1 return pv.UnstructuredGrid( cells.ravel(), np.full(n_cells, 12, dtype=np.uint8), pts_arr, ) def run_one(n: int, integration: str | None) -> float: """Solve Cook's membrane on an n × n hex slab. Returns u_y(C).""" grid = build_grid(n) m = femorph_solver.Model.from_grid(grid) spec = ELEMENTS.HEX8(integration=integration) if integration else ELEMENTS.HEX8 m.assign(spec, material={"EX": E, "PRXY": NU, "DENS": RHO}) pts = np.asarray(m.grid.points) node_nums = np.asarray(m.grid.point_data["ansys_node_num"]) tol = 1e-9 # Clamp left edge AD (x = 0) — pin all three displacement DOFs. left_mask = pts[:, 0] < tol for nn in node_nums[left_mask]: m.fix(nodes=[int(nn)], dof="ALL") # Plane stress is approximated by a thin slab with the through- # thickness u_z left free; nothing extra to do. # Distributed transverse shear on the right edge BC (x = 48). # Total force F_TOTAL spread evenly over the right-edge nodes # at both z = 0 and z = THICKNESS. right_mask = np.abs(pts[:, 0] - 48.0) < tol right_nodes = node_nums[right_mask] fy_per = F_TOTAL / len(right_nodes) for nn in right_nodes: m.apply_force(int(nn), fy=fy_per) res = m.solve_static() g = femorph_solver.io.static_result_to_grid(m, res) pts_g = np.asarray(g.points) # Reference corner C = (48, 60) at any z; pick the closest node. target = np.array([48.0, 60.0]) dist = np.linalg.norm(pts_g[:, :2] - target, axis=1) c_idx = int(np.argmin(dist)) return float(g.point_data["displacement"][c_idx, 1]) # %% # Run all three integration variants on a single coarse mesh # ---------------------------------------------------------- # # We pick :math:`n = 8` (8 × 8 × 1 mesh, 162 nodes total). Cook's # original 1974 paper runs at :math:`2 \times 2`, :math:`4 \times 4`, # and :math:`8 \times 8`; the third is enough to expose the locking # hierarchy without burning too much CPU. N_MESH = 8 print() print(f"On the {N_MESH} × {N_MESH} × 1 mesh ({(N_MESH + 1) ** 2 * 2} nodes):") print() print(f" {'integration':<22} {'u_y(C)':>10} {'u_y / u_ref':>12} {'verdict':<28}") print(f" {'-' * 22:<22} {'-' * 10:>10} {'-' * 12:>12} {'-' * 28:<28}") variants = ( ("plain_gauss", "shear-locks badly"), (None, "default B-bar — partial cure"), ("enhanced_strain", "EAS — fully converged"), ) ratios: dict[str, float] = {} uy_results: dict[str, float] = {} for integ, verdict in variants: uy_c = run_one(N_MESH, integ) ratio = uy_c / UY_C_REF label = integ if integ else "default (B-bar)" ratios[label] = ratio uy_results[label] = uy_c print(f" {label:<22} {uy_c:>10.4f} {ratio:>12.4f} {verdict:<28}") # %% # Verify the textbook ordering # ---------------------------- # # Plain Gauss is known to recover only ~85-92 % on an 8 × 8 mesh # (Cook 1989 Table 6.13.1 / Belytschko 2014 §8.4); B-bar gives a # small bump to ~92-94 % because the lock here is *shear*, not # volumetric. EAS lifts the result to ≥ 98 % (and slightly # overshoots the 2D plane-stress Q4 reference because the 3D slab # is plane-stress-like, not pure plane stress). assert ratios["plain_gauss"] < 0.94, ( f"plain Gauss must lock on Cook's membrane — got {ratios['plain_gauss']:.4f}" ) assert ratios["enhanced_strain"] > 0.98, ( f"EAS must recover ≥ 98 % — got {ratios['enhanced_strain']:.4f}" ) assert ratios["enhanced_strain"] > ratios["plain_gauss"], "EAS must beat plain Gauss" # %% # Refinement ladder — show EAS converges fastest # ---------------------------------------------- ladder = (2, 4, 8, 16) plain: list[float] = [] bbar: list[float] = [] eas: list[float] = [] print() print("Refinement ladder:") print(f" {'n':>4} {'plain_gauss':>12} {'B-bar':>10} {'EAS':>10}") print(f" {'-' * 4:>4} {'-' * 12:>12} {'-' * 10:>10} {'-' * 10:>10}") for n in ladder: pg = run_one(n, "plain_gauss") bb = run_one(n, None) en = run_one(n, "enhanced_strain") plain.append(pg / UY_C_REF) bbar.append(bb / UY_C_REF) eas.append(en / UY_C_REF) print(f" {n:>4} {pg:>12.4f} {bb:>10.4f} {en:>10.4f}") # %% # Render the convergence comparison # --------------------------------- fig, ax = plt.subplots(1, 1, figsize=(7.0, 4.0)) ax.plot(ladder, plain, "o-", color="#d62728", lw=2, label="plain_gauss (shear-locks)") ax.plot(ladder, bbar, "s-", color="#9467bd", lw=2, label="default (B-bar)") ax.plot(ladder, eas, "^-", color="#1f77b4", lw=2, label="enhanced_strain (EAS)") ax.axhline(1.0, color="black", ls="--", lw=1.0, label="mesh-converged limit") ax.set_xlabel("mesh refinement n (n × n × 1)") ax.set_ylabel(r"$u_y(C)\, /\, u_y^{\mathrm{ref}}$") ax.set_xscale("log", base=2) ax.set_xticks(ladder) ax.set_xticklabels([str(n) for n in ladder]) ax.set_title("Cook's membrane — corner displacement convergence", fontsize=11) ax.set_ylim(0.55, 1.05) ax.legend(loc="lower right", fontsize=9) ax.grid(True, ls=":", alpha=0.5) fig.tight_layout() fig.show() # %% # Render the deformed shape (EAS, n=8) # ------------------------------------ grid = build_grid(N_MESH) m_eas = femorph_solver.Model.from_grid(grid) m_eas.assign( ELEMENTS.HEX8(integration="enhanced_strain"), material={"EX": E, "PRXY": NU, "DENS": RHO}, ) pts = np.asarray(m_eas.grid.points) node_nums = np.asarray(m_eas.grid.point_data["ansys_node_num"]) tol = 1e-9 for nn in node_nums[pts[:, 0] < tol]: m_eas.fix(nodes=[int(nn)], dof="ALL") right = node_nums[np.abs(pts[:, 0] - 48.0) < tol] for nn in right: m_eas.apply_force(int(nn), fy=F_TOTAL / len(right)) res_eas = m_eas.solve_static() g_eas = femorph_solver.io.static_result_to_grid(m_eas, res_eas) # Magnify the displacement to make the deformed shape readable. warp_factor = 0.5 # a u_y(C) of ~24 over an L=48 plate is already # visible at unit scale; bump to 0.5× for legibility. warped = g_eas.copy() warped.points = np.asarray(g_eas.points) + warp_factor * np.asarray( g_eas.point_data["displacement"] ) warped["|u|"] = np.linalg.norm(g_eas.point_data["displacement"], axis=1) plotter = pv.Plotter(off_screen=True, window_size=(720, 540)) plotter.add_mesh(g_eas, style="wireframe", color="grey", opacity=0.4, line_width=1) plotter.add_mesh( warped, scalars="|u|", cmap="viridis", show_edges=True, scalar_bar_args={"title": "|u| (×0.5 warp, EAS)"}, ) plotter.add_points( np.array([[48.0, 60.0, 0.5]]), render_points_as_spheres=True, point_size=18, color="#d62728", label=f"corner C — u_y = {uy_results['enhanced_strain']:.3f}", ) plotter.add_legend() plotter.view_xy() plotter.camera.zoom(1.05) plotter.show() # %% # Closing note # ------------ print() print("Take-aways:") print(f" • plain_gauss at n = 8 → {ratios['plain_gauss']:.1%} of converged (shear-locked)") print(f" • default B-bar at n = 8 → {ratios['default (B-bar)']:.1%} (only marginal help)") print(f" • enhanced_strain at n = 8 → {ratios['enhanced_strain']:.1%} (EAS is the cure)") print(" • EAS's edge widens with mesh distortion — pick it for any production hex mesh.")