""" TorsoCAE Validation -- geometric_nonlinear_dynamics case 1: flexible pendulum. *** SUBSTITUTE CASE: swinging plate pendulum vs. independent rigid compound- pendulum ODE, in place of the originally assigned Simo & Vu-Quoc (1986) "flexible whipping cantilever" benchmark. *** Why substituted: the Simo & Vu-Quoc benchmark is a geometrically-exact BEAM-element result whose reference tip-response data lives only in the paper's figures (a WebSearch turned up only a *different*, unrelated Simo-Vu-Quoc tabulated benchmark -- the 45-degree bent cantilever static case, not the dynamic whipping-cantilever one). Digitizing figure curves and mapping them onto TorsoCAE's 3D solid elements (a different discretization than the original geometrically-exact beam) would not give a trustworthy comparison. Substituted, per the same "simpler but still nontrivial" contingency already used for case_02 (double pendulum, see ../../flexible_multibody_dynamics/case_02_slider_crank/README.md): a swinging plate pendulum released from horizontal, geometry and RBE2 pin setup taken directly from the team's own known-good reference journal (journals/plate_pendulum_rbe2_geomnl_dyn1.py), backend switched to dolfinx to match the rest of this validation suite's structural cases. What this case validates: `geometric_nonlinear_dynamics` (Total-Lagrangian, Newmark time integration) combined with a finite-rotation RBE2 pin (`bc("rbe2", ...)`, only rotation about z free) under gravity -- a genuinely large-rotation (0 to ~180 degree swing), nonlinear elastodynamic problem. Compared against an independent rigid COMPOUND-PENDULUM ODE (same mass, same moment of inertia about the pin, same pivot-to-COM distance, computed directly from the plate's own geometry -- not fit), integrated from the same released-from-rest-horizontal initial condition. This validates the finite-rotation kinematics + inertial/gravitational coupling that IS this submodel's defining feature. It does NOT validate flexural/elastic response in isolation -- that is already covered by the modal_analysis and structural_dynamics cases elsewhere in this suite. The steel plate is stiff relative to the gravity loading (self-weight static droop is negligible vs. the ~1m swing amplitude -- see results), so its gross rigid-body-like swing is expected to track the rigid ODE closely, with only a small residual attributable to actual flexibility -- that residual is itself reported, not hidden. Run: HWLOC_COMPONENTS=-gl python3 validation/structural/geometric_nonlinear_dynamics/case_01_flexible_pendulum/run_validation.py """ import json import struct import sys import time from pathlib import Path import numpy as np CASE_DIR = Path(__file__).resolve().parent sys.path.insert(0, str(CASE_DIR.parents[3] / "server")) from pytorsocae import TorsoCAESession import result_frames # ── Geometry (matches journals/plate_pendulum_rbe2_geomnl_dyn1.py exactly) ── LX, LY, LZ = 1.0, 0.2, 0.02 # plate x in [0,1], y in [-0.1,0.1], z in [-0.01,0.01] HOLE_R = 0.04 HOLE_X, HOLE_Y = 0.08, 0.0 # hole (pivot) center E, NU, RHO = 210e9, 0.3, 7800.0 # steel G = 9.81 DT = 0.005 NUM_STEPS = 300 # 1.5 s -- covers the full first half-swing (horizontal to horizontal) with margin # ── Plate (minus hole) mass properties, computed directly from geometry ───── def _mass_properties(): A_full = LX * LY A_hole = np.pi * HOLE_R ** 2 A_net = A_full - A_hole x_cm_full, y_cm_full = 0.5 * LX, 0.0 x_cm = (x_cm_full * A_full - HOLE_X * A_hole) / A_net y_cm = (y_cm_full * A_full - HOLE_Y * A_hole) / A_net mass = A_net * LZ * RHO N = 2000 xs = np.linspace(0.0, LX, N) ys = np.linspace(-0.5 * LY, 0.5 * LY, N) dx, dy = xs[1] - xs[0], ys[1] - ys[0] X, Y = np.meshgrid(xs, ys, indexing="ij") mask = (X - HOLE_X) ** 2 + (Y - HOLE_Y) ** 2 > HOLE_R ** 2 r2_pivot = (X - HOLE_X) ** 2 + (Y - HOLE_Y) ** 2 I_pivot = np.sum(r2_pivot[mask]) * dx * dy * LZ * RHO d_pivot_com = float(np.hypot(x_cm - HOLE_X, y_cm - HOLE_Y)) return {"mass_kg": float(mass), "I_pivot_kg_m2": float(I_pivot), "d_pivot_com_m": d_pivot_com, "x_cm": float(x_cm), "y_cm": float(y_cm)} MP = _mass_properties() # ── independent reference: rigid compound-pendulum ODE ────────────────────── # theta measured from the downward vertical; released from horizontal (pi/2), # at rest -- matches the FEM initial condition exactly. def reference_ode_solution(t_max: float = 1.5): from scipy.integrate import solve_ivp m, I_p, d = MP["mass_kg"], MP["I_pivot_kg_m2"], MP["d_pivot_com_m"] def rhs(_t, y): th, w = y return [w, -(m * G * d / I_p) * np.sin(th)] sol = solve_ivp(rhs, [0, t_max], [np.pi / 2, 0.0], method="DOP853", rtol=1e-12, atol=1e-14, dense_output=True) return sol def decode_scvz(blob: bytes) -> np.ndarray: magic, version, n_verts, has_disp, disp_max, n_fields = struct.unpack_from("<4sIIIfI", blob, 0) if magic != b"SCVZ": raise ValueError(f"Not an SCVZ frame: {magic!r}") off = 24 + n_fields * 80 parts = [("pos", np.float32, (3,))] if has_disp: parts.append(("disp", np.float32, (3,))) if n_fields: parts.append(("scalars", np.float32, (n_fields,))) if version >= 2: parts.append(("mesh_index", np.uint32)) parts.append(("surf_tag", np.uint32)) return np.frombuffer(blob, dtype=np.dtype(parts), count=n_verts, offset=off) def _face_centroid(vbuf: np.ndarray, tag: int) -> np.ndarray: mask = vbuf["surf_tag"] == tag pos = vbuf["pos"][mask].astype(np.float64) disp = vbuf["disp"][mask].astype(np.float64) return (pos + disp).mean(axis=0) def _identify_root_tip_tags(vbuf: np.ndarray) -> tuple[int, int]: """Identify the x=0 (root) and x=LX (tip) face tags by their reference (undeformed) x-extent, rather than assuming a fixed tag numbering -- CSG cut operations do not guarantee stable tag assignment across different geometries/box parameterizations, and picking the wrong tags silently produces a wrong-but-plausible-looking marker (this was caught during development: a corner-vs-center box-offset bug shifted the hole to the plate's edge, and hardcoded tags 1/6 picked up an unrelated small patch near the hole instead of the true end faces).""" tags = sorted(set(vbuf["surf_tag"].tolist())) best_root, best_root_x = None, None best_tip, best_tip_x = None, None for tag in tags: pos = vbuf["pos"][vbuf["surf_tag"] == tag] x_span = float(pos[:, 0].max() - pos[:, 0].min()) if x_span > 1e-9: continue # not a pure x=const face x_val = float(pos[:, 0].mean()) if best_root_x is None or x_val < best_root_x: best_root, best_root_x = tag, x_val if best_tip_x is None or x_val > best_tip_x: best_tip, best_tip_x = tag, x_val if best_root is None or best_tip is None or best_root == best_tip: raise RuntimeError(f"Could not identify distinct root/tip end-cap tags among {tags}") return best_root, best_tip def extract_marker_history(npz_path: str) -> dict: """Root marker = x=0 end-cap face centroid, tip marker = x=LX end-cap face centroid -- any two material points on the rigid plate suffice to recover its rotation angle via atan2 of their relative vector; the actual pivot (hole face) need not be one of them. Tags identified dynamically (see _identify_root_tip_tags) rather than assumed.""" frames = result_frames.frames_from_npz(npz_path) frames = sorted(frames, key=lambda fr: fr["meta"].get("time", 0.0)) root_tag, tip_tag = _identify_root_tip_tags(decode_scvz(frames[0]["bytes"])) rows = [] for fr in frames: t = fr["meta"].get("time") if t is None: continue vbuf = decode_scvz(fr["bytes"]) rows.append((float(t), _face_centroid(vbuf, root_tag), _face_centroid(vbuf, tip_tag))) rows.sort(key=lambda r: r[0]) times = np.array([r[0] for r in rows]) root = np.array([r[1] for r in rows]) tip = np.array([r[2] for r in rows]) return {"time": times, "root": root, "tip": tip, "root_tag": root_tag, "tip_tag": tip_tag} def compare_to_reference(history: dict, sol, windows=(0.5, 1.0, 1.5)) -> dict: t_fem = history["time"] root, tip = history["root"], history["tip"] # atan2(dy,dx) measures from +x (horizontal = 0 rad); convert to the ODE's # convention (theta from the downward vertical, horizontal = +pi/2) via a # constant +pi/2 shift, applied after unwrapping the raw atan2 sequence. angle_fem = np.unwrap(np.arctan2(tip[:, 1] - root[:, 1], tip[:, 0] - root[:, 0])) + np.pi / 2 th_ref = sol.sol(t_fem)[0] rows = [] for w in windows: mask = t_fem <= w rms = float(np.degrees(np.sqrt(np.mean((angle_fem[mask] - th_ref[mask]) ** 2)))) rows.append({"window_s": w, "rms_angle_deg": rms}) return { "windows": rows, "angle_fem_deg": np.degrees(angle_fem).tolist(), "angle_ref_deg": np.degrees(th_ref).tolist(), "time_s": t_fem.tolist(), } def run_case() -> dict: session = TorsoCAESession() session.csg_make_shape("box", {"dx": LX, "dy": LY, "dz": LZ, "cx": 0.0, "cy": -0.5 * LY, "cz": -0.5 * LZ}, name="Plate", body_id="body_plate") session.csg_make_shape("cylinder", {"r": HOLE_R, "h": 5 * LZ, "cx": HOLE_X, "cy": HOLE_Y, "cz": -2.5 * LZ}, name="HoleTool", body_id="body_hole") session.csg_boolean("cut", "body_plate", "body_hole", name="SwingPlate", body_id="body_0") session.csg_select("body_0") session.mesh(mesh_id="mesh_0", algo_id="delaunay", size_factor=0.3, size_min=0, size_max=0, order=1, dim=3) session.solid("SwingPlate").material(E=E, nu=NU, rho=RHO) session.surface(7, scope_id="body_0").bc( "rbe2", reference=[HOLE_X, HOLE_Y, 0.0], components=["x", "y", "z"], independent_components=["rz"], ) session.set_physics("structural", submodel="geometric_nonlinear_dynamics", backend="dolfinx") session.set_model_options(gravity=[0.0, -G, 0.0]) session.set_solver_options( algo="direct", precond="none", tol=1e-06, max_iter=900, num_steps=NUM_STEPS, dt=DT, time_scheme="newmark", max_inner_iter=25, inner_tol=1e-04, n_cores=2, device="cpu", ) t0 = time.time() result = session.compute(mesh_ids=["mesh_0"]) elapsed = time.time() - t0 history = extract_marker_history(result["npz_path"]) sol = reference_ode_solution(t_max=NUM_STEPS * DT) comparison = compare_to_reference(history, sol, windows=(0.5, 1.0, 1.5)) primary = comparison["windows"][1] # 1.0 s: covers the full first half-swing rms_deg = primary["rms_angle_deg"] return { "case": "geometric_nonlinear_dynamics / swinging plate pendulum vs. independent " "rigid compound-pendulum ODE (SUBSTITUTE for Simo & Vu-Quoc 1986 whipping " "cantilever -- see module docstring for why)", "substitution_reason": ( "Simo & Vu-Quoc (1986) 'flexible whipping cantilever' reference data exists " "only in the paper's figures (a different, unrelated SVQ benchmark -- the " "45-degree bent cantilever STATIC case -- is the one with tabulated numbers " "commonly cited online); digitizing figure curves for a geometrically-exact " "beam-element result and mapping onto TorsoCAE's 3D solid discretization would " "not give a trustworthy comparison. Substituted with a swinging plate pendulum " "(same geometry/RBE2 setup as the team's own known-good reference journal, " "journals/plate_pendulum_rbe2_geomnl_dyn1.py), validated via trajectory " "comparison against an independent rigid-body ODE -- same substitution pattern " "already used and accepted for flexible_multibody_dynamics/case_02." ), "geometry": {"plate_m": [LX, LY, LZ], "hole_r_m": HOLE_R, "hole_center_xy": [HOLE_X, HOLE_Y]}, "material": {"E_Pa": E, "nu": NU, "rho_kg_m3": RHO}, "mass_properties": MP, "solver": {"time_scheme": "newmark", "dt_s": DT, "num_steps": NUM_STEPS, "backend": "dolfinx"}, "method": "Root/tip markers (x=0 and x=1 face centroids) tracked over time from the " "SCVZ result-frame history; rotation angle recovered via atan2 of their " "relative vector (any two rigid-body material points suffice -- the pivot " "itself, tag 7, need not be one of the markers). Compared against an " "independent nonlinear rigid compound-pendulum ODE (scipy solve_ivp DOP853, " "rtol=1e-12) using the plate's own mass/I_pivot/pivot-to-COM distance, " "computed directly from its geometry (not fit), released from the same " "rest-horizontal initial condition.", "rms_angle_error_by_window": comparison["windows"], "primary_window_s": primary["window_s"], "rms_angle_error_deg": rms_deg, "final_error_pct": rms_deg, # degrees, not %, kept for the shared results.json key convention "pass": bool(rms_deg < 5.0), # deg, over the 1.0s primary window (full first half-swing) "elapsed_s": elapsed, "time_history": { "time_s": history["time"].tolist(), "angle_fem_deg": comparison["angle_fem_deg"], "angle_ref_deg": comparison["angle_ref_deg"], }, } def plot_response(summary: dict) -> None: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt th = summary["time_history"] t = np.array(th["time_s"]) fig, ax = plt.subplots(figsize=(7.5, 4.8), dpi=150) ax.plot(t, th["angle_fem_deg"], "-", color="#2563eb", linewidth=1.5, label="TorsoCAE FEM (plate, RBE2 pin)") ax.plot(t, th["angle_ref_deg"], "--", color="#dc2626", linewidth=1.3, label="Rigid compound-pendulum ODE (reference)") primary_w = summary["primary_window_s"] ax.axvline(primary_w, color="#9ca3af", linestyle=":", linewidth=1.0, label=f"RMS window (<= {primary_w:.1f}s): {summary['rms_angle_error_deg']:.3f} deg") ax.set_xlabel("Time (s)") ax.set_ylabel("Swing angle from vertical (deg)") ax.set_title("Flexible plate pendulum: TorsoCAE FEM vs. independent rigid-body ODE") ax.legend(loc="best", fontsize=8) ax.grid(alpha=0.3) fig.tight_layout() plots_dir = CASE_DIR / "plots" plots_dir.mkdir(exist_ok=True) fig.savefig(plots_dir / "pendulum_swing.png") def main(): summary = run_case() history = summary.pop("time_history") (CASE_DIR / "results.json").write_text(json.dumps(summary, indent=2)) print(json.dumps(summary, indent=2), flush=True) summary["time_history"] = history plot_response(summary) if __name__ == "__main__": main()