Hyperbolic-paraboloid cantilever: restored warped-quad consistency

This rectangular hyperbolic paraboloid is a deliberately controlled #7 fixture. Its mid-surface is z = c*x*(y-W/2): it has negative Gaussian curvature, and, unlike the sphere, the mixed-curvature (warp) rate c does not disappear when the structured mesh is refined. A clamp at x=0 and a total global-z tip force give a simple, reproducible scalar response at the centre of the free edge.

The physical target is a converged triangular solve on exactly this geometry. Triangles have planar cells and are objective for each candidate curvature formulation, so this is an independently computed reference rather than an assumed published number. At n = 32, 48, 64 (where the mesh is 2n x n), the triangle answers are 0.92081419, 0.92107296, 0.92117998. The final two differ by only 1.16e-4 relatively, so the n=64 result is used as the reference. The old warped-quad sequence was 1.1166539, 1.1233583, 1.1257320 at n=8,16,32: it converged 22.2% above, and away from, that target. With Hermit’s constant-normal curvature the sequence is 0.9171184, 0.9199326, 0.9208417; the quad solution now converges to the independently computed triangle result.

"""Hyperbolic-paraboloid cantilever: restored warped-quad consistency

This rectangular hyperbolic paraboloid is a deliberately controlled #7 fixture.
Its mid-surface is ``z = c*x*(y-W/2)``: it has negative Gaussian curvature, and,
unlike the sphere, the mixed-curvature (warp) rate ``c`` does not disappear when the
structured mesh is refined.  A clamp at ``x=0`` and a total global-z tip force give
a simple, reproducible scalar response at the centre of the free edge.

The physical target is a converged *triangular* solve on exactly this geometry.
Triangles have planar cells and are objective for each candidate curvature
formulation, so this is an independently computed reference rather than an assumed
published number.  At ``n = 32, 48, 64`` (where the mesh is ``2n x n``), the triangle
answers are ``0.92081419, 0.92107296, 0.92117998``.  The final two differ by only
``1.16e-4`` relatively, so the ``n=64`` result is used as the reference.  The
old warped-quad sequence was ``1.1166539, 1.1233583, 1.1257320`` at
``n=8,16,32``: it converged 22.2% above, and away from, that target. With Hermit's
constant-normal curvature the sequence is ``0.9171184, 0.9199326, 0.9208417``;
the quad solution now converges to the independently computed triangle result.
"""

import pathlib
import sys

import csdl_alpha as csdl
import numpy as np

import hermit as hm

sys.path.insert(0, str(pathlib.Path(__file__).parent))
from _geometry import hyperbolic_paraboloid              # noqa: E402
from _harness import Case, main, node_nearest             # noqa: E402

LENGTH, WIDTH, WARP = 4.0, 2.0, 0.5
E, NU, THICKNESS, LOAD = 1.0e5, 0.25, 0.1, 1.0


def solve_response(n, *, cell, warp=WARP, thickness=THICKNESS):
    """Free-edge-centre z deflection on a ``2n x n`` hypar mesh."""
    mesh = hyperbolic_paraboloid(LENGTH, WIDTH, warp, nx=2 * n, ny=n, cell=cell)
    rec = csdl.Recorder(inline=True)
    rec.start()
    domain = hm.ShellDomain(mesh, element="CG2CG1")
    material = hm.isotropic(domain, E=E, nu=NU, thickness=thickness, density=1.0)
    xyz = np.asarray(domain.node_coords)
    tip = np.flatnonzero(np.isclose(xyz[:, 0], LENGTH))
    loads = sum((hm.point_load(domain, at=xyz[k], force=[0.0, 0.0, LOAD / len(tip)])
                 for k in tip), start=hm.Loads(domain))
    state = hm.solve(domain, material, loads,
                     hm.clamp(domain, where=lambda x: np.isclose(x[0], 0.0)))
    u = hm.nodal_displacement(state).value.reshape(-1, 3)
    rec.stop()
    k, distance = node_nearest(domain, [LENGTH, WIDTH / 2, 0.0])
    if distance > 1e-10:
        raise RuntimeError(f"free-edge sample moved by {distance:.3e}")
    return float(u[k, 2])


def solve_at(n):
    """Warped-quad response, compared with the separately converged triangle target."""
    return solve_response(n, cell="quad")


# Independent triangle result at n=64; see the module docstring for the convergence
# evidence.  This is deliberately a literal result of a separate solve, not fitted
# from the quad sequence or a target adjusted to its answer.
REFERENCE = 0.9211799767764

CASE = Case(
    name="Hyperbolic-paraboloid cantilever (warped quadrilaterals)",
    quantity="global-z deflection at the centre of the free edge",
    reference=REFERENCE,
    tolerance=0.02,
    citation="Computed in this file: triangle n=64 = 0.92117998; n=48 differs by "
             "1.16e-4 relative on the same geometry. Triangles are planar and "
             "objective under the compared curvature formulations",
    levels=(8, 16, 32),
    quick_level=8,
    solve=solve_at,
    monotone=False,
    notes="the former #7 result was 1.1257320 at n=32 (22.2% high); the "
          "constant-normal result is 0.9208417 (0.037% low).",
)

if __name__ == "__main__":
    main(CASE)