# Warped quadrilateral cells preserve rigid-body consistency **This case records a resolved defect and guards the formulation choice.** Hermit's constant-normal bending curvature now gives warped quadrilaterals the same effective flat-facet theory as triangles. Before that change, a doubly-curved surface whose quadrilateral cells are *warped* -- four corners that do not lie in a plane -- converged cleanly and monotonically to a wrong answer. The fixture is the MacNeal-Harder twisted beam, whose reference is known (``5.424e-3``, see ``ex_twisted_beam.py``). A twisted strip is a ruled but *non-developable* surface, so a structured quad mesh on it is necessarily warped, while a triangle mesh on the same nodes is not -- every triangle is planar by construction. Meshing the same geometry both ways isolates the cell shape as the only difference: =========== ================== ================== refinement triangles warped quads =========== ================== ================== 12 x 4 5.3879e-3 (0.993) 7.6097e-3 (1.403) 24 x 6 5.4081e-3 (0.997) 7.6150e-3 (1.404) 48 x 8 5.4133e-3 (0.998) 7.6168e-3 (1.404) =========== ================== ================== Triangles converged to the reference. Quads converged to a value **40 % too large**, and refinement moved it *away* from the reference rather than toward it -- the signature of a consistency error, not a discretisation error. With the constant-normal measure, the same sweep is: =========== ================== ================== ============== refinement triangles warped quads quad / triangle =========== ================== ================== ============== 12 x 4 5.3879e-3 (0.993) 5.4024e-3 (0.996) 1.00270 24 x 6 5.4081e-3 (0.997) 5.4122e-3 (0.998) 1.00075 48 x 8 5.4133e-3 (0.998) 5.4148e-3 (0.998) 1.00026 =========== ================== ================== ============== The quad and triangle sequences now converge together and both approach the MacNeal-Harder reference. The quad answer is not merely inaccurate, it is unphysical: untwisted beam theory gives ``P L^3 / (3 E I)`` with ``I = w t^3 / 12``, or ``6.6125e-3``, for a tip load perpendicular to the wide face. A 90-degree twist rotates the stiff axis into the load path and can only *stiffen* that response, so any correct answer must be below ``6.6125e-3``. The published ``5.424e-3`` is; the quad result, ``7.617e-3``, is 15 % above a bound it cannot legitimately exceed. Root cause of the former result -- established, not conjectured. The old curvature measure (``ElasticModel._bending_curvature``) was:: kappa = sym(gradv_local(grad(cross(E2, theta)), E01)) Under a rigid-body motion the rotation field is a constant vector ``theta = omega``, so by the product rule ``grad(cross(E2, omega)) = cross(grad(E2), omega)``. On a planar cell ``E2`` is constant, ``grad(E2) = 0``, and ``kappa`` vanishes identically -- as it must. On a warped cell ``E2`` varies *within* the cell, so ``kappa`` is **non-zero for a constant rotation field**: the measure differentiates the shell normal along with the rotation. The element therefore stores energy under rigid-body motion, which is a consistency violation, not an accuracy limitation. Measured directly, as spurious elastic energy under an exact rigid motion ``u = c + omega x x``, ``theta = omega`` (normalised by ``E t^3``): ========== =========== =========== warp quad triangle ========== =========== =========== 0 8.5e-30 2.2e-29 1e-4 2.7e-9 2.3e-29 1e-3 2.7e-7 2.3e-29 1e-2 2.7e-5 2.4e-29 ========== =========== =========== It scales as ``warp**2``, exactly as ``kappa`` proportional to ``grad(E2)`` predicts, and splitting the energy by term puts **all** of it in the bending term -- membrane, shear and drilling stay at ~1e-30 for both cell types at every warp. Independently ruled out: quadrature (forcing degrees 2 to 20 changes nothing) and the drilling stabilisation (scaling it over four orders of magnitude moves the result 0.2 %). Hermit now holds ``E2`` constant in this measure and constructs curvature from ``grad(theta)`` alone. A constant rigid rotation therefore gives zero curvature structurally. Triangles are no longer required merely because geometry is doubly curved; quads and triangles now implement the same flat-facet theory. The remaining limitation is physical rather than cell-specific: moderately thick, strongly curved, membrane-loaded shells may need a full covariant formulation (see the user guide). The gate below asserts that the warped-quad and triangle answers agree; the old 1.407 ratio fails it by a wide margin. conda activate hermit python examples/verification/ex_warped_quad_consistency.py ```python """Warped quadrilateral cells preserve rigid-body consistency **This case records a resolved defect and guards the formulation choice.** Hermit's constant-normal bending curvature now gives warped quadrilaterals the same effective flat-facet theory as triangles. Before that change, a doubly-curved surface whose quadrilateral cells are *warped* -- four corners that do not lie in a plane -- converged cleanly and monotonically to a wrong answer. The fixture is the MacNeal-Harder twisted beam, whose reference is known (``5.424e-3``, see ``ex_twisted_beam.py``). A twisted strip is a ruled but *non-developable* surface, so a structured quad mesh on it is necessarily warped, while a triangle mesh on the same nodes is not -- every triangle is planar by construction. Meshing the same geometry both ways isolates the cell shape as the only difference: =========== ================== ================== refinement triangles warped quads =========== ================== ================== 12 x 4 5.3879e-3 (0.993) 7.6097e-3 (1.403) 24 x 6 5.4081e-3 (0.997) 7.6150e-3 (1.404) 48 x 8 5.4133e-3 (0.998) 7.6168e-3 (1.404) =========== ================== ================== Triangles converged to the reference. Quads converged to a value **40 % too large**, and refinement moved it *away* from the reference rather than toward it -- the signature of a consistency error, not a discretisation error. With the constant-normal measure, the same sweep is: =========== ================== ================== ============== refinement triangles warped quads quad / triangle =========== ================== ================== ============== 12 x 4 5.3879e-3 (0.993) 5.4024e-3 (0.996) 1.00270 24 x 6 5.4081e-3 (0.997) 5.4122e-3 (0.998) 1.00075 48 x 8 5.4133e-3 (0.998) 5.4148e-3 (0.998) 1.00026 =========== ================== ================== ============== The quad and triangle sequences now converge together and both approach the MacNeal-Harder reference. The quad answer is not merely inaccurate, it is unphysical: untwisted beam theory gives ``P L^3 / (3 E I)`` with ``I = w t^3 / 12``, or ``6.6125e-3``, for a tip load perpendicular to the wide face. A 90-degree twist rotates the stiff axis into the load path and can only *stiffen* that response, so any correct answer must be below ``6.6125e-3``. The published ``5.424e-3`` is; the quad result, ``7.617e-3``, is 15 % above a bound it cannot legitimately exceed. Root cause of the former result -- established, not conjectured. The old curvature measure (``ElasticModel._bending_curvature``) was:: kappa = sym(gradv_local(grad(cross(E2, theta)), E01)) Under a rigid-body motion the rotation field is a constant vector ``theta = omega``, so by the product rule ``grad(cross(E2, omega)) = cross(grad(E2), omega)``. On a planar cell ``E2`` is constant, ``grad(E2) = 0``, and ``kappa`` vanishes identically -- as it must. On a warped cell ``E2`` varies *within* the cell, so ``kappa`` is **non-zero for a constant rotation field**: the measure differentiates the shell normal along with the rotation. The element therefore stores energy under rigid-body motion, which is a consistency violation, not an accuracy limitation. Measured directly, as spurious elastic energy under an exact rigid motion ``u = c + omega x x``, ``theta = omega`` (normalised by ``E t^3``): ========== =========== =========== warp quad triangle ========== =========== =========== 0 8.5e-30 2.2e-29 1e-4 2.7e-9 2.3e-29 1e-3 2.7e-7 2.3e-29 1e-2 2.7e-5 2.4e-29 ========== =========== =========== It scales as ``warp**2``, exactly as ``kappa`` proportional to ``grad(E2)`` predicts, and splitting the energy by term puts **all** of it in the bending term -- membrane, shear and drilling stay at ~1e-30 for both cell types at every warp. Independently ruled out: quadrature (forcing degrees 2 to 20 changes nothing) and the drilling stabilisation (scaling it over four orders of magnitude moves the result 0.2 %). Hermit now holds ``E2`` constant in this measure and constructs curvature from ``grad(theta)`` alone. A constant rigid rotation therefore gives zero curvature structurally. Triangles are no longer required merely because geometry is doubly curved; quads and triangles now implement the same flat-facet theory. The remaining limitation is physical rather than cell-specific: moderately thick, strongly curved, membrane-loaded shells may need a full covariant formulation (see the user guide). The gate below asserts that the warped-quad and triangle answers agree; the old 1.407 ratio fails it by a wide margin. conda activate hermit python examples/verification/ex_warped_quad_consistency.py """ import pathlib import sys import numpy as np import csdl_alpha as csdl import hermit as hm sys.path.insert(0, str(pathlib.Path(__file__).parent)) from _geometry import twisted_strip # noqa: E402 from _harness import Case, main # noqa: E402 LENGTH, WIDTH, THICKNESS, TWIST = 12.0, 1.1, 0.32, 90.0 E, NU, LOAD = 29.0e6, 0.22, 1.0 PUBLISHED = 5.424e-3 # Tip load perpendicular to the wide face on the UNTWISTED beam. A 90-degree twist # can only stiffen this, so any correct twisted answer lies below it. UNTWISTED_BOUND = LOAD * LENGTH**3 / (3.0 * E * WIDTH * THICKNESS**3 / 12.0) def _tip_deflection(nx, ny, cell): rec = csdl.Recorder(inline=True) rec.start() domain = hm.ShellDomain(twisted_strip(LENGTH, WIDTH, TWIST, nx=nx, ny=ny, cell=cell), element="CG2CG1") material = hm.isotropic(domain, E=E, nu=NU, thickness=THICKNESS, density=1.0) bcs = hm.clamp(domain, where=lambda x: np.isclose(x[0], 0.0), method="strong") xyz = np.asarray(domain.node_coords) tip = np.flatnonzero(np.isclose(xyz[:, 0], LENGTH)) load = None for k in tip: term = hm.point_load(domain, at=xyz[k], force=[0.0, 0.0, LOAD / len(tip)]) load = term if load is None else load + term state = hm.solve(domain, material, load, bcs) u = hm.nodal_displacement(state).value.reshape(-1, 3) rec.stop() k = int(np.argmin(np.abs(xyz[:, 0] - LENGTH) + np.abs(xyz[:, 1]) + np.abs(xyz[:, 2]))) return float(u[k, 2]) def solve_at(n): """Ratio of the warped-quad answer to the triangle answer on the same geometry.""" ny = n // 6 + 2 tri = _tip_deflection(n, ny, "triangle") quad = _tip_deflection(n, ny, "quad") print(f" triangles={tri:.6e} ({tri / PUBLISHED:.4f} of ref) " f"warped quads={quad:.6e} ({quad / PUBLISHED:.4f} of ref)") assert quad <= UNTWISTED_BOUND, ( f"warped-quad answer {quad:.6e} exceeds untwisted bound {UNTWISTED_BOUND:.6e}") return quad / tri CASE = Case( name="Warped quadrilateral rigid-body consistency", quantity="ratio of the warped-quad tip deflection to the triangle tip deflection", reference=1.0, tolerance=0.02, citation="Measured on this fixture; the underlying twisted-beam reference is " "MacNeal & Harder, Finite Elements in Analysis and Design 1(1), 1985", levels=(12, 24, 48), quick_level=24, solve=solve_at, monotone=False, notes="the pre-fix ratio converged to 1.407; the constant-normal curvature makes " "rigid-body objectivity structural and the ratio now converges to 1.", ) if __name__ == "__main__": main(CASE) ```