{ "cells": [ { "cell_type": "markdown", "id": "82467c68", "metadata": {}, "source": [ "Benchmark: MacNeal-Harder standard problems for the shell elements\n", "===\n", "\n", "Date: 22 of October 2026\n", "\n", "Author: Saullo G. P. Castro\n", "\n", "MacNeal and Harder [1] proposed a set of small problems that excite, with coarse and distorted meshes, the failure modes of shell and solid elements: spurious mechanisms and locking. This notebook runs the problems of that set that a flat shell element can model with the four shell elements of `pyfe3d`, `Quad4`, `Quad4R`, `Tria3R` and `Tria3DSG`, every element with its default attributes, and compares the normalised results with the theoretical values of the paper and with the results of the MSC/NASTRAN QUAD4 element tabulated there.\n", "\n", "| Problem of [1] | Figure, tables | Here |\n", "|---|---|---|\n", "| Patch test for plates | Fig. 2, Tables 2 and 7 | membrane and bending |\n", "| Straight cantilever beam | Fig. 4, Tables 3 and 8 | regular, trapezoidal and parallelogram meshes |\n", "| Curved beam | Fig. 5, Tables 3 and 9 | in-plane and out-of-plane load |\n", "| Twisted beam | Fig. 6, Tables 3 and 10 | in-plane and out-of-plane load |\n", "| Rectangular plate | Fig. 7, Tables 4, 11 and 12 | the two cases of Tables 11 and 12 |\n", "| Scordelis-Lo roof, spherical shell | Figs. 8 and 9, Tables 5, 13 and 14 | in the notebook of the shell obstacle course |\n", "| Patch test for solids, thick-walled cylinder | Figs. 3 and 10 | not modelled: solid and plane strain problems |\n", "\n", "The triangular meshes are obtained by cutting every quadrilateral of the mesh of [1] along its diagonal 1-3, as the paper suggests for triangles (p. 7). Unless said otherwise, a displacement is normalised by the theoretical value of [1], and the letter grades follow the rule of [1], p. 14: A for an error of at most 2%, B up to 10%, C up to 20%, D up to 50% and F above, the absolute errors being averaged when several cases contribute to one grade." ] }, { "cell_type": "markdown", "id": "57b459db", "metadata": {}, "source": [ "Common code\n", "---\n", "\n", "`stiffness()` assembles the linear stiffness matrix `KC0` of a mesh of quadrilaterals for one of the four elements, `solve()` solves the static problem with prescribed degrees-of-freedom, `tip_forces()` distributes a unit end load as in Fig. 11(a) of [1], half of the load of each element edge going to each of its two nodes." ] }, { "cell_type": "code", "execution_count": 1, "id": "d1e0b467", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:37.931234Z", "iopub.status.busy": "2026-10-25T09:19:37.931071Z", "iopub.status.idle": "2026-10-25T09:19:38.280391Z", "shell.execute_reply": "2026-10-25T09:19:38.279650Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from scipy.sparse import coo_matrix\n", "from scipy.sparse.linalg import spsolve\n", "\n", "from pyfe3d.shellprop_utils import isotropic_plate\n", "from pyfe3d import (Quad4, Quad4Data, Quad4Probe, Quad4R, Quad4RData,\n", " Quad4RProbe, Tria3R, Tria3RData, Tria3RProbe, Tria3DSG,\n", " Tria3DSGData, Tria3DSGProbe, INT, DOUBLE, DOF)\n", "\n", "ELEMENTS = {\n", " 'Quad4': (Quad4, Quad4Data, Quad4Probe),\n", " 'Quad4R': (Quad4R, Quad4RData, Quad4RProbe),\n", " 'Tria3R': (Tria3R, Tria3RData, Tria3RProbe),\n", " 'Tria3DSG': (Tria3DSG, Tria3DSGData, Tria3DSGProbe),\n", "}\n", "\n", "\n", "def split(quads):\n", " # two triangles per quadrilateral, cut along the diagonal 1-3\n", " return [t for (n1, n2, n3, n4) in quads for t in ((n1, n2, n3), (n1, n3, n4))]\n", "\n", "\n", "def stiffness(ncoords, quads, name, prop):\n", " Elem, Data, Probe = ELEMENTS[name]\n", " conn = split(quads) if name.startswith('Tria') else quads\n", " data, probe = Data(), Probe()\n", " ncf = np.ascontiguousarray(ncoords, dtype=DOUBLE).flatten()\n", " N = DOF*ncoords.shape[0]\n", " size = data.KC0_SPARSE_SIZE*len(conn)\n", " KC0r = np.zeros(size, dtype=INT)\n", " KC0c = np.zeros(size, dtype=INT)\n", " KC0v = np.zeros(size, dtype=DOUBLE)\n", " areas = np.zeros(len(conn))\n", " for k, nodes in enumerate(conn):\n", " e = Elem(probe)\n", " for i, n in enumerate(nodes):\n", " setattr(e, 'n%d' % (i + 1), n)\n", " setattr(e, 'c%d' % (i + 1), DOF*n)\n", " e.init_k_KC0 = k*data.KC0_SPARSE_SIZE\n", " e.update_rotation_matrix(ncf)\n", " e.update_probe_xe(ncf)\n", " e.update_area()\n", " e.update_KC0(KC0r, KC0c, KC0v, prop)\n", " areas[k] = e.area\n", " K = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc()\n", " return K, conn, areas\n", "\n", "\n", "def solve(K, bk, f, uk=None):\n", " # static solution with the degrees-of-freedom of the mask bk prescribed to uk\n", " u = np.zeros(K.shape[0])\n", " if uk is not None:\n", " u[bk] = uk[bk]\n", " bu = ~bk\n", " u[bu] = spsolve(K[bu][:, bu], f[bu] - K[bu][:, bk] @ u[bk])\n", " return u\n", "\n", "\n", "def area_load(ncoords, conn, areas, q):\n", " # nodal forces of a load per unit area q, lumped equally to the nodes\n", " f = np.zeros(DOF*ncoords.shape[0])\n", " for nodes, A in zip(conn, areas):\n", " for n in nodes:\n", " f[DOF*n:DOF*n + 3] += np.asarray(q)*A/len(nodes)\n", " return f\n", "\n", "\n", "def grid(nx, ny):\n", " # node ids of a structured grid and its quadrilaterals\n", " ids = np.arange((nx + 1)*(ny + 1)).reshape(nx + 1, ny + 1)\n", " quads = [(ids[i, j], ids[i + 1, j], ids[i + 1, j + 1], ids[i, j + 1])\n", " for i in range(nx) for j in range(ny)]\n", " return ids, quads\n", "\n", "\n", "def tip_forces(N, tip, comp, total=1.):\n", " f = np.zeros(N)\n", " for a, b in zip(tip[:-1], tip[1:]):\n", " f[DOF*a + comp] += total/(len(tip) - 1)/2\n", " f[DOF*b + comp] += total/(len(tip) - 1)/2\n", " return f\n", "\n", "\n", "def grade(normalised):\n", " # letter grade of [1], p. 14, from the mean absolute error\n", " err = 100*np.mean(np.abs(np.atleast_1d(normalised) - 1.))\n", " for letter, limit in (('A', 2.), ('B', 10.), ('C', 20.), ('D', 50.)):\n", " if err <= limit:\n", " return letter\n", " return 'F'" ] }, { "cell_type": "markdown", "id": "0ca78cd2", "metadata": {}, "source": [ "Patch tests for plates\n", "---\n", "\n", "The patch of Fig. 2 of [1] is a rectangle of 0.24 by 0.12 divided into five quadrilaterals by four interior nodes at (0.04, 0.02), (0.18, 0.03), (0.16, 0.08) and (0.08, 0.08), with $t = 0.001$, $E = 10^6$ and $\\nu = 0.25$. The displacements of the exact field of Table 2 of [1] are prescribed at the four corner nodes and the interior nodes are free:\n", "\n", "- membrane: $u = 10^{-3}(x + y/2)$, $v = 10^{-3}(y + x/2)$, a state of constant strain $\\epsilon_x = \\epsilon_y = \\gamma_{xy} = 10^{-3}$;\n", "- bending: $w = 10^{-3}(x^2 + xy + y^2)/2$, with the rotations of zero transverse shear, which in the convention of `pyfe3d` are $r_x = \\partial w/\\partial y$ and $r_y = -\\partial w/\\partial x$, a state of constant curvature.\n", "\n", "The out-of-plane degrees-of-freedom are fixed in the membrane test, and $u$, $v$ and $r_z$ in the bending test. The table of [1] prescribes $u$ and $v$ only, for elements without drilling rotations. `Quad4` and `Tria3R` enrich their in-plane displacements with the drilling rotations [2], so that the displacement of a boundary edge depends on $r_z$ at its nodes, and the drilling rotation of the boundary nodes is prescribed as well, to the in-plane rotation of the exact field, $(v_{,x} - u_{,y})/2 = 0$. Left free, it lets the boundary edges deform beyond the prescribed field.\n", "\n", "[1] records the maximum error of the stresses (Table 7). Stresses are recovered differently by each element of `pyfe3d`, so the measures here are the maximum error of the displacements of the interior nodes, relative to the largest value of the exact field, and the error of the strain energy $U = \\frac{1}{2} \\boldsymbol{u}^\\top \\boldsymbol{K} \\boldsymbol{u}$ against the exact $\\frac{1}{2} \\boldsymbol{\\epsilon}^\\top \\boldsymbol{A} \\boldsymbol{\\epsilon}$ (membrane) or $\\frac{1}{2} \\boldsymbol{\\kappa}^\\top \\boldsymbol{D} \\boldsymbol{\\kappa}$ (bending) times the area. Both are zero, to round-off, when the element passes." ] }, { "cell_type": "code", "execution_count": 2, "id": "683a1d3c", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:38.283228Z", "iopub.status.busy": "2026-10-25T09:19:38.282946Z", "iopub.status.idle": "2026-10-25T09:19:38.309513Z", "shell.execute_reply": "2026-10-25T09:19:38.308731Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "element test interior displ. error energy error check\n", "Quad4 membrane 2.82e-16 -1.11e-16 PASS\n", "Quad4 bending 7.39e-01 5.80e+00 FAIL\n", "Quad4R membrane 4.52e-16 -1.11e-16 PASS\n", "Quad4R bending 7.24e-01 1.34e+01 FAIL\n", "Tria3R membrane 6.78e-16 -4.44e-16 PASS\n", "Tria3R bending 8.41e-02 2.43e-01 FAIL\n", "Tria3DSG membrane 1.08e-15 -1.11e-16 PASS\n", "Tria3DSG bending 3.82e-14 -2.38e-11 PASS\n" ] } ], "source": [ "PATCH_XY = np.array([[0., 0.], [0.24, 0.], [0.24, 0.12], [0., 0.12],\n", " [0.04, 0.02], [0.18, 0.03], [0.16, 0.08], [0.08, 0.08]])\n", "PATCH_QUADS = [(4, 5, 6, 7), (0, 1, 5, 4), (1, 2, 6, 5), (2, 3, 7, 6), (3, 0, 4, 7)]\n", "\n", "\n", "def patch_test(name, form):\n", " ncoords = np.column_stack((PATCH_XY, np.zeros(8)))\n", " prop = isotropic_plate(thickness=0.001, E=1.e6, nu=0.25)\n", " K, conn, areas = stiffness(ncoords, PATCH_QUADS, name, prop)\n", " N = K.shape[0]\n", " x, y = PATCH_XY.T\n", " boundary = np.arange(8) < 4\n", " uex = np.zeros(N)\n", " bk = np.zeros(N, dtype=bool)\n", " if form == 'membrane':\n", " uex[0::DOF] = 1e-3*(x + y/2)\n", " uex[1::DOF] = 1e-3*(y + x/2)\n", " checked, fixed, prescribed = (0, 1), (2, 3, 4), (0, 1, 5)\n", " C = np.array([[prop.A11, prop.A12, prop.A16],\n", " [prop.A12, prop.A22, prop.A26],\n", " [prop.A16, prop.A26, prop.A66]])\n", " else:\n", " uex[2::DOF] = 1e-3*(x**2 + x*y + y**2)/2\n", " uex[3::DOF] = 1e-3*(y + x/2)\n", " uex[4::DOF] = -1e-3*(x + y/2)\n", " checked, fixed, prescribed = (2, 3, 4), (0, 1, 5), (2, 3, 4)\n", " C = np.array([[prop.D11, prop.D12, prop.D16],\n", " [prop.D12, prop.D22, prop.D26],\n", " [prop.D16, prop.D26, prop.D66]])\n", " for d in fixed:\n", " bk[d::DOF] = True\n", " for d in prescribed:\n", " bk[d::DOF] |= boundary\n", " u = solve(K, bk, np.zeros(N), uex)\n", " err_u = max(np.abs(u[d::DOF][~boundary] - uex[d::DOF][~boundary]).max()\n", " /np.abs(uex[d::DOF]).max() for d in checked)\n", " # the three strains, or the three curvatures, are all equal to 1e-3 in magnitude\n", " e = np.full(3, 1e-3)\n", " U_exact = 0.5*e @ C @ e*0.24*0.12\n", " err_U = 0.5*u @ (K @ u)/U_exact - 1.\n", " return err_u, err_U\n", "\n", "\n", "patch = {}\n", "print('element test interior displ. error energy error check')\n", "for name in ELEMENTS:\n", " for form in ('membrane', 'bending'):\n", " err_u, err_U = patch[name, form] = patch_test(name, form)\n", " ok = err_u < 1e-8 and abs(err_U) < 1e-8\n", " print('%-10s %-10s %20.2e %14.2e %s' % (name, form, err_u, err_U,\n", " 'PASS' if ok else 'FAIL'))" ] }, { "cell_type": "markdown", "id": "adb967f8", "metadata": {}, "source": [ "All four elements pass the membrane patch test. Only `Tria3DSG` passes the bending patch test. [1] reports 0.0% for the QUAD4 of MSC/NASTRAN in both tests (Table 7).\n", "\n", "The failure of `Quad4` and `Quad4R` in bending is in the transverse shear. The exact transverse displacement is quadratic and both quadrilaterals interpolate it bilinearly, with the transverse shear sampled at the centroid; on these distorted elements the gradient of the bilinear interpolant at the centroid differs from the prescribed rotation there, on every element of the patch, which leaves a transverse shear strain and, with $\\ell/t \\approx 100$, a large parasitic shear energy. `Tria3R` takes its transverse shear at the centroid as well, softened by its `alpha_shear_locking` factor. The discrete shear gap of `Tria3DSG` makes the shear strain of a state of constant curvature vanish exactly, which is the property this test checks." ] }, { "cell_type": "markdown", "id": "ce5abf7a", "metadata": {}, "source": [ "Straight cantilever beam\n", "---\n", "\n", "Fig. 4 of [1]: length 6.0, width 0.2, thickness 0.1, $E = 10^7$, $\\nu = 0.3$, a 6 by 1 mesh clamped at one end, with regular, trapezoidal and parallelogram elements, the interior node lines of the distorted meshes being slanted at 45 degrees. The unit loads are applied at the free end, the twist as a couple of transverse forces $\\pm 1/0.2$ on the two tip nodes (Fig. 11(b) of [1]). The extension and the shear results are the mean tip displacement in the direction of the load, the twist result is the tip rotation $(w_2 - w_1)/0.2$." ] }, { "cell_type": "code", "execution_count": 3, "id": "2d8328e2", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:38.311844Z", "iopub.status.busy": "2026-10-25T09:19:38.311648Z", "iopub.status.idle": "2026-10-25T09:19:38.412181Z", "shell.execute_reply": "2026-10-25T09:19:38.411297Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normalised tip displacement\n", "mesh load Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "regular extension 0.999 0.998 0.992 0.992 0.995\n", "regular in-plane shear 0.904 0.119 0.145 0.027 0.904\n", "regular out-of-plane shear 0.980 0.973 0.876 0.971 0.986\n", "regular twist 377.368 4.780 1.669 0.675 0.941\n", "trapezoidal extension 0.998 0.999 0.992 0.991 0.996\n", "trapezoidal in-plane shear 0.761 0.033 0.149 0.013 0.071\n", "trapezoidal out-of-plane shear 0.981 0.233 0.589 0.959 0.968\n", "trapezoidal twist 109.715 4.907 3.244 0.658 0.951\n", "parallelogram extension 0.998 0.999 0.992 0.991 0.996\n", "parallelogram in-plane shear 0.838 0.068 0.102 0.010 0.080\n", "parallelogram out-of-plane shear 0.982 0.772 0.831 0.958 0.977\n", "parallelogram twist 137.968 5.258 2.709 0.638 0.945\n" ] } ], "source": [ "STRAIGHT_REF = {'extension': 3.0e-5, 'in-plane shear': 0.1081,\n", " 'out-of-plane shear': 0.4321, 'twist': 0.03208} # Table 3 of [1]\n", "MH_QUAD4_STRAIGHT = { # Table 8 of [1]\n", " 'regular': {'extension': 0.995, 'in-plane shear': 0.904, 'out-of-plane shear': 0.986, 'twist': 0.941},\n", " 'trapezoidal': {'extension': 0.996, 'in-plane shear': 0.071, 'out-of-plane shear': 0.968, 'twist': 0.951},\n", " 'parallelogram': {'extension': 0.996, 'in-plane shear': 0.080, 'out-of-plane shear': 0.977, 'twist': 0.945},\n", "}\n", "\n", "\n", "def straight_beam(name, shape, load):\n", " L, w, nx, ny = 6., 0.2, 6, 1\n", " ids, quads = grid(nx, ny)\n", " ys = np.linspace(0., w, ny + 1)\n", " ncoords = np.zeros(((nx + 1)*(ny + 1), 3))\n", " for i in range(nx + 1):\n", " d = 0.\n", " if 0 < i < nx and shape == 'parallelogram':\n", " d = 0.1\n", " elif 0 < i < nx and shape == 'trapezoidal':\n", " d = 0.1 if i % 2 else -0.1\n", " ncoords[ids[i, :], 0] = i*L/nx + d*(2*ys/w - 1)\n", " ncoords[ids[i, :], 1] = ys\n", " prop = isotropic_plate(thickness=0.1, E=1.e7, nu=0.3)\n", " K, conn, areas = stiffness(ncoords, quads, name, prop)\n", " N = K.shape[0]\n", " bk = np.zeros(N, dtype=bool)\n", " for d in range(DOF):\n", " bk[d::DOF] |= np.isclose(ncoords[:, 0], 0.)\n", " tip = ids[-1, :]\n", " if load == 'twist':\n", " f = np.zeros(N)\n", " f[DOF*tip[0] + 2] = -1/w\n", " f[DOF*tip[-1] + 2] = +1/w\n", " u = solve(K, bk, f)\n", " return (u[DOF*tip[-1] + 2] - u[DOF*tip[0] + 2])/w\n", " comp = {'extension': 0, 'in-plane shear': 1, 'out-of-plane shear': 2}[load]\n", " u = solve(K, bk, tip_forces(N, tip, comp))\n", " return u[DOF*tip + comp].mean()\n", "\n", "\n", "straight = {}\n", "print('normalised tip displacement')\n", "print('%-14s %-19s' % ('mesh', 'load') + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')\n", "for shape in ('regular', 'trapezoidal', 'parallelogram'):\n", " for load, ref in STRAIGHT_REF.items():\n", " for name in ELEMENTS:\n", " straight[name, shape, load] = straight_beam(name, shape, load)/ref\n", " print('%-14s %-19s' % (shape, load)\n", " + ''.join('%10.3f' % straight[name, shape, load] for name in ELEMENTS)\n", " + '%10.3f' % MH_QUAD4_STRAIGHT[shape][load])" ] }, { "cell_type": "markdown", "id": "a5da3d82", "metadata": {}, "source": [ "Extension is reproduced by every element on every mesh. The other cases show the following failures:\n", "\n", "- **Twist of `Quad4`.** The one-point transverse shear of `Quad4`, from Hughes et al. [3], leaves the hourglass mode $w = (x - x_c)(y - y_c)$ without stiffness, and a strip one element wide loaded in torsion deforms in that mode, so the twist is two orders of magnitude too large. `Quad4R` controls that mode with its hourglass stiffness, and the result depends on the amount of it, here about five times too flexible. Both are documented in the Notes of the two modules, and the tests `test_quad4_spurious_shear_mode.py` and `test_drilling_benchmarks.py` show that two elements across the width recover the twist with `Quad4`.\n", "- **In-plane shear, i.e. in-plane bending.** `Quad4`, with the drilling enrichment of Allman [2], reaches 90% of the reference on the regular mesh, like the QUAD4 of [1], and keeps 76% to 84% on the distorted meshes, where the QUAD4 of [1] locks (7% to 8%). `Quad4R` locks on all three meshes: its one-point membrane gives no strain to the in-plane bending pattern of the bilinear element, which is resisted only by its hourglass stabilisation and by its drilling term, and on this one-element-deep strip these make it far too stiff. The triangles lock too: `Tria3DSG` has the constant-strain membrane of the linear triangle, and `Tria3R`, which carries the Allman enrichment, reaches 10% to 15%.\n", "- **Out-of-plane shear on distorted meshes.** `Quad4R` drops to 23% on the trapezoidal mesh and `Tria3R` to 59%.\n", "- **Twist of the triangles.** `Tria3R` is 1.7 to 3.2 times too flexible and `Tria3DSG` 32% to 36% too stiff." ] }, { "cell_type": "markdown", "id": "61877a74", "metadata": {}, "source": [ "Curved beam\n", "---\n", "\n", "Fig. 5 of [1]: a 90 degree arc of inner radius 4.12 and outer radius 4.32, thickness 0.1, $E = 10^7$, $\\nu = 0.25$, a 6 by 1 mesh clamped at one end and loaded by unit forces at the tip, in the plane of the arc (radially at the tip) or normal to it." ] }, { "cell_type": "code", "execution_count": 4, "id": "67fd86de", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:38.415125Z", "iopub.status.busy": "2026-10-25T09:19:38.414880Z", "iopub.status.idle": "2026-10-25T09:19:38.439467Z", "shell.execute_reply": "2026-10-25T09:19:38.438636Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normalised tip displacement\n", "load Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "in-plane 0.908 0.103 0.119 0.022 0.833\n", "out-of-plane 26.828 1.587 1.100 0.760 0.951\n" ] } ], "source": [ "CURVED_REF = {'in-plane': 0.08734, 'out-of-plane': 0.5022} # Table 3 of [1]\n", "MH_QUAD4_CURVED = {'in-plane': 0.833, 'out-of-plane': 0.951} # Table 9 of [1]\n", "\n", "\n", "def curved_beam(name, load):\n", " nx, ny = 6, 1\n", " th = np.linspace(0., np.pi/2, nx + 1)\n", " r = np.linspace(4.12, 4.32, ny + 1)\n", " T, Rr = np.meshgrid(th, r, indexing='ij')\n", " ncoords = np.column_stack(((Rr*np.cos(T)).ravel(), (Rr*np.sin(T)).ravel(),\n", " np.zeros(T.size)))\n", " ids, quads = grid(nx, ny)\n", " prop = isotropic_plate(thickness=0.1, E=1.e7, nu=0.25)\n", " K, conn, areas = stiffness(ncoords, quads, name, prop)\n", " N = K.shape[0]\n", " bk = np.zeros(N, dtype=bool)\n", " for d in range(DOF):\n", " bk[d::DOF] |= np.isclose(ncoords[:, 1], 0.)\n", " tip = ids[-1, :]\n", " comp = 1 if load == 'in-plane' else 2\n", " u = solve(K, bk, tip_forces(N, tip, comp))\n", " return u[DOF*tip + comp].mean()\n", "\n", "\n", "curved = {}\n", "print('normalised tip displacement')\n", "print('%-14s' % 'load' + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')\n", "for load, ref in CURVED_REF.items():\n", " for name in ELEMENTS:\n", " curved[name, load] = curved_beam(name, load)/ref\n", " print('%-14s' % load + ''.join('%10.3f' % curved[name, load] for name in ELEMENTS)\n", " + '%10.3f' % MH_QUAD4_CURVED[load])" ] }, { "cell_type": "markdown", "id": "4e54f3ac", "metadata": {}, "source": [ "The in-plane case repeats the in-plane bending result of the straight beam: `Quad4` reaches 91%, better than the 83% of the QUAD4 of [1], and the other three elements lock. The out-of-plane load twists the strip, and `Quad4` is again far too flexible through the hourglass mode of its transverse shear." ] }, { "cell_type": "markdown", "id": "f51a76bd", "metadata": {}, "source": [ "Twisted beam\n", "---\n", "\n", "Fig. 6 of [1]: length 12.0, width 1.1, thickness 0.32, twisted by 90 degrees from root to tip, $E = 29 \\times 10^6$, $\\nu = 0.22$, a 12 by 2 mesh clamped at the root, with unit forces at the tip. The reference surface is a helicoid, so every quadrilateral is warped, by 7.5 degrees according to [1] (p. 11). At the tip the width lies along $z$: the in-plane load is along $z$ and the out-of-plane load along $y$.\n", "\n", "The quadrilaterals of `pyfe3d` take their normal from the cross product of the two diagonals and work with the coordinates of the nodes in that frame (`update_rotation_matrix()` and `update_probe_xe()`), with no correction for warping; a triangle is flat by construction. The mesh is refined beyond the one of [1] to show the effect." ] }, { "cell_type": "code", "execution_count": 5, "id": "30f03bf7", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:38.442405Z", "iopub.status.busy": "2026-10-25T09:19:38.442155Z", "iopub.status.idle": "2026-10-25T09:19:58.538695Z", "shell.execute_reply": "2026-10-25T09:19:58.538043Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normalised tip displacement\n", "mesh load Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "12x2 in-plane 0.252 0.234 0.910 0.969 0.993\n", "12x2 out-of-plane 0.332 0.325 0.884 0.817 0.985\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "48x4 in-plane 0.659 0.691 0.967 0.994\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "48x4 out-of-plane 0.711 0.743 0.965 0.954\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "192x8 in-plane 0.899 0.946 0.994 0.998\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "192x8 out-of-plane 0.917 0.957 0.994 0.990\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "768x16 in-plane 0.973 0.992 0.998 0.999\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "768x16 out-of-plane 0.978 0.994 0.998 0.997\n" ] } ], "source": [ "TWISTED_REF = {'in-plane': 0.005424, 'out-of-plane': 0.001754} # Table 3 of [1]\n", "MH_QUAD4_TWISTED = {'in-plane': 0.993, 'out-of-plane': 0.985} # Table 10 of [1]\n", "\n", "\n", "def twisted_beam(name, load, nx=12, ny=2):\n", " L, w = 12., 1.1\n", " X, S = np.meshgrid(np.linspace(0., L, nx + 1), np.linspace(-w/2, w/2, ny + 1),\n", " indexing='ij')\n", " beta = np.pi/2*X/L\n", " ncoords = np.column_stack((X.ravel(), (S*np.cos(beta)).ravel(),\n", " (S*np.sin(beta)).ravel()))\n", " ids, quads = grid(nx, ny)\n", " prop = isotropic_plate(thickness=0.32, E=29.e6, nu=0.22)\n", " K, conn, areas = stiffness(ncoords, quads, name, prop)\n", " N = K.shape[0]\n", " bk = np.zeros(N, dtype=bool)\n", " for d in range(DOF):\n", " bk[d::DOF] |= np.isclose(ncoords[:, 0], 0.)\n", " tip = ids[-1, :]\n", " comp = 2 if load == 'in-plane' else 1\n", " u = solve(K, bk, tip_forces(N, tip, comp))\n", " return u[DOF*tip + comp].mean()\n", "\n", "\n", "twisted = {}\n", "print('normalised tip displacement')\n", "print('%-9s %-14s' % ('mesh', 'load') + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')\n", "for nx, ny in ((12, 2), (48, 4), (192, 8), (768, 16)):\n", " for load, ref in TWISTED_REF.items():\n", " for name in ELEMENTS:\n", " twisted[name, nx, load] = twisted_beam(name, load, nx, ny)/ref\n", " print('%-9s %-14s' % ('%dx%d' % (nx, ny), load)\n", " + ''.join('%10.3f' % twisted[name, nx, load] for name in ELEMENTS)\n", " + ('%10.3f' % MH_QUAD4_TWISTED[load] if nx == 12 else ''))" ] }, { "cell_type": "markdown", "id": "10832ba3", "metadata": {}, "source": [ "On the mesh of [1] the two quadrilaterals are three to four times too stiff, while the QUAD4 of [1] is within 2%. The triangles, which are flat by construction, are within 20%. The error of the quadrilaterals decreases as the warping of each element decreases with the refinement, and all four elements approach the reference of [1] on the finest mesh. MacNeal [4], one of the two references given for the QUAD4 in Table 6 of [1], modifies the matrices of the plane element for warped quadrilaterals; nothing equivalent is implemented in the quadrilaterals of `pyfe3d`." ] }, { "cell_type": "markdown", "id": "85328937", "metadata": {}, "source": [ "Rectangular plate\n", "---\n", "\n", "Fig. 7 of [1]: $a = 2.0$, $b = 2.0$ or $10.0$, thickness 0.0001, $E = 1.7472 \\times 10^7$, $\\nu = 0.3$, an $N \\times N$ mesh on one quarter of the plate with symmetry conditions on the two inner edges. The two cases tabulated by [1] for the QUAD4 are computed: simply supported edges with a uniform pressure $q = 10^{-4}$ (Table 11) and clamped edges with a central load $P = 4 \\times 10^{-4}$ (Table 12), the quarter model carrying $P/4$. Simply supported means $w = 0$ on the edge, with free rotations; clamped means $w = r_x = r_y = 0$. The in-plane degrees-of-freedom $u$, $v$ and $r_z$ are fixed, the plate being loaded in bending only. The pressure is lumped to the nodes, one quarter (one third for triangles) of the load of each element going to each node. The result is the deflection at the centre normalised by Table 4 of [1]." ] }, { "cell_type": "code", "execution_count": 6, "id": "778ea97c", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:58.540263Z", "iopub.status.busy": "2026-10-25T09:19:58.540153Z", "iopub.status.idle": "2026-10-25T09:19:58.728958Z", "shell.execute_reply": "2026-10-25T09:19:58.728510Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "simple supports, uniform load, b/a = 1, normalised deflection at the centre\n", "N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "2 1.065 0.781 1.410 0.883 0.981\n", "4 1.016 0.930 1.069 0.989 1.004\n", "6 1.007 0.967 0.995 1.009 1.003\n", "8 1.004 0.981 0.965 1.017 1.002\n", "\n", "simple supports, uniform load, b/a = 5, normalised deflection at the centre\n", "N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "2 0.826 0.875 4.583 0.922 1.052\n", "4 0.978 0.974 1.868 1.115 0.991\n", "6 0.991 0.988 1.391 1.155 0.997\n", "8 0.995 0.993 1.225 1.170 0.998\n", "\n", "clamped supports, central load, b/a = 1, normalised deflection at the centre\n", "N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "2 0.865 0.533 2.286 0.000 0.934\n", "4 0.965 0.819 1.401 0.000 1.010\n", "6 0.985 0.907 1.170 0.002 1.012\n", "8 0.992 0.944 1.074 0.009 1.010\n", "\n", "clamped supports, central load, b/a = 5, normalised deflection at the centre\n", "N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "2 0.318 0.287 8.366 0.000 0.519\n", "4 0.825 0.666 3.942 0.000 0.863\n", "6 0.911 0.799 2.591 0.001 0.940\n", "8 0.946 0.867 1.994 0.003 0.972\n", "\n" ] } ], "source": [ "PLATE_REF = {('simple', 1.): 4.062, ('simple', 5.): 12.97, # Table 4 of [1]\n", " ('clamped', 1.): 5.60, ('clamped', 5.): 7.23}\n", "MH_QUAD4_PLATE = {('simple', 1.): [0.981, 1.004, 1.003, 1.002], # Table 11 of [1]\n", " ('simple', 5.): [1.052, 0.991, 0.997, 0.998],\n", " ('clamped', 1.): [0.934, 1.010, 1.012, 1.010], # Table 12 of [1]\n", " ('clamped', 5.): [0.519, 0.863, 0.940, 0.972]}\n", "NS = (2, 4, 6, 8)\n", "\n", "\n", "def rect_plate(name, n, ratio, support):\n", " a, t, E, nu = 2., 0.0001, 1.7472e7, 0.3\n", " b = ratio*a\n", " X, Y = np.meshgrid(np.linspace(0., a/2, n + 1), np.linspace(0., b/2, n + 1),\n", " indexing='ij')\n", " ncoords = np.column_stack((X.ravel(), Y.ravel(), np.zeros(X.size)))\n", " ids, quads = grid(n, n)\n", " prop = isotropic_plate(thickness=t, E=E, nu=nu)\n", " K, conn, areas = stiffness(ncoords, quads, name, prop)\n", " N = K.shape[0]\n", " x, y = ncoords[:, 0], ncoords[:, 1]\n", " bk = np.zeros(N, dtype=bool)\n", " for d in (0, 1, 5):\n", " bk[d::DOF] = True\n", " edge = np.isclose(x, 0.) | np.isclose(y, 0.)\n", " bk[2::DOF] |= edge\n", " if support == 'clamped':\n", " bk[3::DOF] |= edge\n", " bk[4::DOF] |= edge\n", " bk[4::DOF] |= np.isclose(x, a/2) # symmetry, rotation about y\n", " bk[3::DOF] |= np.isclose(y, b/2) # symmetry, rotation about x\n", " centre = ids[n, n]\n", " if support == 'simple':\n", " f = area_load(ncoords, conn, areas, (0., 0., 1.e-4))\n", " else:\n", " f = np.zeros(N)\n", " f[DOF*centre + 2] = 4.e-4/4\n", " u = solve(K, bk, f)\n", " return u[DOF*centre + 2]\n", "\n", "\n", "plate = {}\n", "for (support, ratio), ref in PLATE_REF.items():\n", " print('%s, b/a = %g, normalised deflection at the centre'\n", " % ('simple supports, uniform load' if support == 'simple'\n", " else 'clamped supports, central load', ratio))\n", " print('%-4s' % 'N' + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')\n", " for i, n in enumerate(NS):\n", " for name in ELEMENTS:\n", " plate[name, support, ratio, n] = rect_plate(name, n, ratio, support)/ref\n", " print('%-4d' % n + ''.join('%10.3f' % plate[name, support, ratio, n] for name in ELEMENTS)\n", " + '%10.3f' % MH_QUAD4_PLATE[support, ratio][i])\n", " print()" ] }, { "cell_type": "markdown", "id": "cb36680b", "metadata": {}, "source": [ "`Quad4` converges in all four cases, more slowly than the QUAD4 of [1] for the clamped plate with $b/a = 5$. `Quad4R` converges from below, more slowly. The two triangles are the interesting cases at this thickness, $\\ell/t$ between 1250 and 25000:\n", "\n", "- `Tria3R` is too flexible on the coarse meshes, up to eight times on the elongated elements of $b/a = 5$, and is not converged at $N = 8$. This is the transverse shear stiffness divided by $1 + \\alpha \\ell^2/t^2$ with the default `alpha_shear_locking = 0.7`, whose error behaves as $O(1/N^2)$ (Notes of `pyfe3d.tria3r`).\n", "- `Tria3DSG` locks with clamped supports and the central load, at less than 1% of the reference up to $N = 8$, and with simple supports and the uniform load it is too flexible on the elongated elements of $b/a = 5$, by 17% at $N = 8$. Both converge under further refinement, computed below, but slowly." ] }, { "cell_type": "code", "execution_count": 7, "id": "40859bfb", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:19:58.730932Z", "iopub.status.busy": "2026-10-25T09:19:58.730795Z", "iopub.status.idle": "2026-10-25T09:20:00.829531Z", "shell.execute_reply": "2026-10-25T09:20:00.828894Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "case N Quad4 Tria3DSG\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "clamped, b/a = 1 16 0.999 0.328\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "clamped, b/a = 1 32 1.001 0.913\n", "simple, b/a = 5 16 0.999 1.177\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "simple, b/a = 5 32 1.000 1.033\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "simple, b/a = 5 64 1.000 1.000\n" ] } ], "source": [ "print('%-36s %4s %10s %10s' % ('case', 'N', 'Quad4', 'Tria3DSG'))\n", "for support, ratio, ns in (('clamped', 1., (16, 32)), ('simple', 5., (16, 32, 64))):\n", " for n in ns:\n", " print('%-36s %4d' % ('%s, b/a = %g' % (support, ratio), n)\n", " + ''.join('%10.3f' % (rect_plate(name, n, ratio, support)/PLATE_REF[support, ratio])\n", " for name in ('Quad4', 'Tria3DSG')))" ] }, { "cell_type": "markdown", "id": "9cd6b5c6", "metadata": {}, "source": [ "Summary\n", "---\n", "\n", "The grades of Table 16 of [1] for the tests modelled here, the absolute errors of the cases of each test being averaged. For the patch tests the error is the larger of the two measures above, in per cent. As a check of the grading, the same function applied to the QUAD4 results of Tables 8 to 12 of [1] must give the QUAD4 column of their Table 16." ] }, { "cell_type": "code", "execution_count": 8, "id": "0ddc9992", "metadata": { "execution": { "iopub.execute_input": "2026-10-25T09:20:00.831042Z", "iopub.status.busy": "2026-10-25T09:20:00.830933Z", "iopub.status.idle": "2026-10-25T09:20:00.841030Z", "shell.execute_reply": "2026-10-25T09:20:00.840243Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "grading of the QUAD4 results of [1] reproduces their Table 16: PASS\n", "\n", "test Quad4 Quad4R Tria3R Tria3DSG MH QUAD4\n", "(1) Patch test, membrane A A A A A\n", "(2) Patch test, bending F F D A A\n", "(3) Straight beam, extension A A A A A\n", "(4) Straight beam, in-plane, regular B F F F B\n", "(5) Straight beam, in-plane, irregular D F F F F\n", "(6) Straight beam, out-of-plane, regular A B C B A\n", "(7) Straight beam, out-of-plane, irregular A D D B B\n", "(8) Straight beam, twist F F F D B\n", "(9) Curved beam, in-plane B F F F C\n", "(10) Curved beam, out-of-plane F F B D B\n", "(11) Twisted beam F F C C A\n", "(12) Rectangular plate (N = 4) B C F F B\n", "failed tests (D and F) 5 8 7 6 1\n" ] } ], "source": [ "def summary(r):\n", " # r(test) returns the normalised results of the cases of each test\n", " tests = [\n", " ('(1) Patch test, membrane', r('patch', 'membrane')),\n", " ('(2) Patch test, bending', r('patch', 'bending')),\n", " ('(3) Straight beam, extension', r('straight', 'extension', ('regular', 'trapezoidal', 'parallelogram'))),\n", " ('(4) Straight beam, in-plane, regular', r('straight', 'in-plane shear', ('regular',))),\n", " ('(5) Straight beam, in-plane, irregular', r('straight', 'in-plane shear', ('trapezoidal', 'parallelogram'))),\n", " ('(6) Straight beam, out-of-plane, regular', r('straight', 'out-of-plane shear', ('regular',))),\n", " ('(7) Straight beam, out-of-plane, irregular', r('straight', 'out-of-plane shear', ('trapezoidal', 'parallelogram'))),\n", " ('(8) Straight beam, twist', r('straight', 'twist', ('regular', 'trapezoidal', 'parallelogram'))),\n", " ('(9) Curved beam, in-plane', r('curved', 'in-plane')),\n", " ('(10) Curved beam, out-of-plane', r('curved', 'out-of-plane')),\n", " ('(11) Twisted beam', r('twisted')),\n", " ('(12) Rectangular plate (N = 4)', r('plate')),\n", " ]\n", " return [(label, grade(v) if v is not None else '-') for label, v in tests]\n", "\n", "\n", "def ours(name):\n", " def r(test, *args):\n", " if test == 'patch':\n", " # error in per cent turned into a normalised value\n", " return 1. + max(abs(e) for e in patch[name, args[0]])\n", " if test == 'straight':\n", " load, shapes = args\n", " return [straight[name, s, load] for s in shapes]\n", " if test == 'curved':\n", " return curved[name, args[0]]\n", " if test == 'twisted':\n", " return [twisted[name, 12, l] for l in TWISTED_REF]\n", " if test == 'plate':\n", " return [plate[name, s, ra, 4] for (s, ra) in PLATE_REF]\n", " return r\n", "\n", "\n", "def macneal_harder_quad4(test, *args):\n", " if test == 'patch':\n", " return 1. # Table 7 of [1]: 0.0% in both tests\n", " if test == 'straight':\n", " load, shapes = args\n", " return [MH_QUAD4_STRAIGHT[s][load] for s in shapes]\n", " if test == 'curved':\n", " return MH_QUAD4_CURVED[args[0]]\n", " if test == 'twisted':\n", " return list(MH_QUAD4_TWISTED.values())\n", " if test == 'plate':\n", " return [MH_QUAD4_PLATE[k][1] for k in PLATE_REF]\n", "\n", "\n", "# QUAD4 column of Table 16 of [1], tests (1) to (12)\n", "TABLE16_QUAD4 = list('AAABFABBCBAB')\n", "check = [g for _, g in summary(macneal_harder_quad4)]\n", "print('grading of the QUAD4 results of [1] reproduces their Table 16:',\n", " 'PASS' if check == TABLE16_QUAD4 else 'FAIL %s' % check)\n", "print()\n", "grades = {name: summary(ours(name)) for name in ELEMENTS}\n", "labels = [label for label, _ in grades['Quad4']]\n", "print('%-44s' % 'test' + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')\n", "for i, label in enumerate(labels):\n", " print('%-44s' % label + ''.join('%10s' % grades[n][i][1] for n in ELEMENTS)\n", " + '%10s' % TABLE16_QUAD4[i])\n", "print('%-44s' % 'failed tests (D and F)'\n", " + ''.join('%10d' % sum(g in 'DF' for _, g in grades[n]) for n in ELEMENTS)\n", " + '%10d' % sum(g in 'DF' for g in TABLE16_QUAD4))" ] }, { "cell_type": "markdown", "id": "b2df7be9", "metadata": {}, "source": [ "With their default attributes, the four elements of `pyfe3d` fail between five and eight of these twelve tests, against one for the QUAD4 of [1]. The failures of `Quad4` are the bending patch test, the in-plane bending of the distorted straight beam (D, where the QUAD4 of [1] gets an F), the two tests that twist a strip one element wide, and the warped elements of the twisted beam. `Quad4R` and the two triangles lock in in-plane bending on these one-element-deep meshes, and `Tria3DSG` is the only element that passes both patch tests. Several failures are confined to the coarse meshes that [1] prescribes on purpose: the refinements above show the twisted beam converging for every element and the locked plates of `Tria3DSG` recovering.\n", "\n", "The Scordelis-Lo roof and the spherical shell of [1] are run, with Tables 13 and 14 of [1], in the notebook of the shell obstacle course.\n", "\n", "References\n", "---\n", "\n", "[1] MacNeal, R. H., and Harder, R. L., A proposed standard set of problems to test finite element accuracy, Finite Elements in Analysis and Design, 1(1), 3-20, 1985. https://doi.org/10.1016/0168-874X(85)90003-4\n", "\n", "[2] Allman, D. J., A compatible triangular element including vertex rotations for plane elasticity analysis, Computers & Structures, 19(1-2), 1-8, 1984. https://doi.org/10.1016/0045-7949(84)90197-4\n", "\n", "[3] Hughes, T. J. R., Taylor, R. L., and Kanoknukulchai, W., A simple and efficient finite element for plate bending, International Journal for Numerical Methods in Engineering, 11(10), 1529-1543, 1977. https://doi.org/10.1002/nme.1620111005\n", "\n", "[4] MacNeal, R. H., A simple quadrilateral shell element, Computers & Structures, 8(2), 175-183, 1978. https://doi.org/10.1016/0045-7949(78)90020-2" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.13.13" } }, "nbformat": 4, "nbformat_minor": 5 }