#!/usr/bin/env python3 """Reviewed, local, fixed-schema elastic FE worker. No networking or external commands. Input: a LignoSparseFEJob28 JSON produced by the research workspace. Output: a source-bound numerical receipt. This is not signed execution attestation. NumPy and SciPy must be installed by the computer owner. """ import argparse import hashlib import json import math import os from pathlib import Path import sys import tempfile import warnings VERSION = "lingo-sparse-fe28-v1" MAX_INPUT = 2 * 1024 * 1024 MAX_OUTPUT = 24 * 1024 * 1024 MAX_JOB = 4 * 1024 * 1024 + 4096 def digest(value): return hashlib.sha256(value.encode("utf-8") if isinstance(value, str) else value).hexdigest() def strict_pairs(pairs): result = {} for key, value in pairs: if key in result: raise ValueError("Duplicate JSON object field.") result[key] = value return result def strict_json(text): return json.loads(text, object_pairs_hook=strict_pairs, parse_constant=lambda _: (_ for _ in ()).throw(ValueError("Nonfinite JSON number."))) def finite_tree(value, depth=0): if depth > 35: raise ValueError("Source nesting exceeds the supported bound.") if isinstance(value, float) and not math.isfinite(value): raise ValueError("Source contains a nonfinite value.") if isinstance(value, dict): for child in value.values(): finite_tree(child, depth + 1) if isinstance(value, list): for child in value: finite_tree(child, depth + 1) def run_model(model): import numpy as np import scipy from scipy.sparse import coo_matrix, csr_matrix from scipy.sparse.linalg import spsolve, MatrixRankWarning expected = {"format", "version", "title", "illustrative", "units", "nodes", "elements", "materials", "boundary", "loads", "ties", "submodel", "dependsOn", "sourceRefs", "reviewed", "externalProcessing", "nativeReceipt"} if set(model) != expected or model["format"] != "LignoElements28" or model["version"] != 1 or model["reviewed"] is not True: raise ValueError("Use an exact reviewed LignoElements28 source object.") if model["units"] != {"length": "mm", "force": "N", "stress": "MPa", "rotation": "rad"}: raise ValueError("The worker requires explicit N-mm-MPa-rad units.") if model["externalProcessing"] is not True or model["nativeReceipt"] is not None: raise ValueError("Prepared source needs explicit external-processing consent and an empty receipt slot.") if model["submodel"] is not None: raise ValueError("Native sparse submodel transfer is withheld; transfer and verify original parent fields in the browser first.") nodes, elements = model["nodes"], model["elements"] if not 2 <= len(nodes) <= 2000 or not 1 <= len(elements) <= 6000 or not 1 <= len(model["materials"]) <= 64: raise ValueError("Source exceeds the worker's 2000-node/6000-element/64-material limit.") if len(model["loads"]) > 12000 or len(model["boundary"]) > 12000 or len(model["ties"]) > 12000: raise ValueError("Source exceeds the fixed load/constraint budget.") node_index = {n["id"]: i for i, n in enumerate(nodes)} materials = {m["id"]: m for m in model["materials"]} if len(node_index) != len(nodes) or len(materials) != len(model["materials"]) or len({e["id"] for e in elements}) != len(elements): raise ValueError("Original node, element and material IDs must be distinct.") def numeric(value): return type(value) in (int, float) and math.isfinite(value) def vector(value): return isinstance(value, list) and len(value) == 3 and all(numeric(v) for v in value) for node in nodes: if set(node) != {"id", "xyz"} or not isinstance(node["id"], str) or not vector(node["xyz"]): raise ValueError("Invalid original node record.") for material in materials.values(): keys = {"id", "E1", "E2", "E3", "nu12", "nu23", "nu13", "G12", "G23", "G13", "axis1", "axis2"} if set(material) != keys or not vector(material["axis1"]) or not vector(material["axis2"]) or not all(numeric(material[k]) for k in keys - {"id", "axis1", "axis2"}): raise ValueError("Invalid orthotropic engineering constants.") for entry in model["boundary"] + model["loads"]: if set(entry) != {"node", "dof", "value"} or not numeric(entry["value"]): raise ValueError("Invalid original load or prescribed value.") for tie in model["ties"]: if set(tie) != {"slave", "terms", "offset"} or not numeric(tie["offset"]) or abs(tie["offset"]) > 1e6: raise ValueError("Invalid multipoint offset.") if set(tie["slave"]) != {"node", "dof"} or any(set(t) != {"node", "dof", "coefficient"} or not numeric(t["coefficient"]) for t in tie["terms"]): raise ValueError("Invalid multipoint source term.") for element in elements: allowed = {"id", "type", "nodes", "materialId", "plyId", "thicknessMm", "plies", "section"} if not {"id", "type", "nodes", "materialId"}.issubset(element) or set(element) - allowed: raise ValueError("Unsupported element source fields.") if "thicknessMm" in element and (not numeric(element["thicknessMm"]) or element["thicknessMm"] <= 0): raise ValueError("Invalid shell thickness.") for ply in element.get("plies", []): if set(ply) != {"materialId", "thicknessMm", "angleDeg"} or not numeric(ply["thicknessMm"]) or not numeric(ply["angleDeg"]) or ply["thicknessMm"] <= 0 or abs(ply["angleDeg"]) > 360: raise ValueError("Invalid original ply definition.") if "section" in element: section = element["section"] if set(section) != {"areaMm2", "iyMm4", "izMm4", "torsionMm4", "localY"} or not vector(section["localY"]) or not all(numeric(section[k]) and section[k] > 0 for k in set(section) - {"localY"}): raise ValueError("Invalid original beam section.") coords = np.asarray([n["xyz"] for n in nodes], dtype=float) if coords.shape != (len(nodes), 3) or np.max(np.abs(coords)) > 1e7: raise ValueError("Invalid or excessive original coordinate.") def axes(a, b): x = np.asarray(a, dtype=float) if np.linalg.norm(x) < 1e-12: raise ValueError("Degenerate material/element direction.") x = x / np.linalg.norm(x) z = np.cross(x, np.asarray(b, dtype=float)) if np.linalg.norm(z) < 1e-12: raise ValueError("Parallel material/element directions.") z /= np.linalg.norm(z) return np.asarray([x, np.cross(z, x), z]) def elasticity(m): moduli = [m[k] for k in ["E1", "E2", "E3", "G23", "G13", "G12"]] if min(moduli) <= 1e-9 or max(moduli) > 1e7: raise ValueError("Invalid elastic modulus.") if any(abs(m[k]) >= .99 for k in ["nu12", "nu23", "nu13"]): raise ValueError("Invalid Poisson ratio.") s = np.diag(1 / np.asarray(moduli)) s[0, 1] = s[1, 0] = -m["nu12"] / m["E1"] s[0, 2] = s[2, 0] = -m["nu13"] / m["E1"] s[1, 2] = s[2, 1] = -m["nu23"] / m["E2"] np.linalg.cholesky(s) return np.linalg.inv(s) def strain_rotation(r): t = np.zeros((6, 6)) for i in range(6): e = np.zeros(6) e[i] = 1 tensor = np.array([[e[0], e[5] / 2, e[4] / 2], [e[5] / 2, e[1], e[3] / 2], [e[4] / 2, e[3] / 2, e[2]]]) local = r @ tensor @ r.T t[:, i] = [local[0, 0], local[1, 1], local[2, 2], 2 * local[1, 2], 2 * local[0, 2], 2 * local[0, 1]] return t signs = np.array([[-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1], [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]]) g = 1 / math.sqrt(3) minimum_jacobian = math.inf coo_i, coo_j, coo_v = [], [], [] touched = set() def shell_shape(x, y, a, b): ss = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1]]) n = (1 + ss[:, 0] * x) * (1 + ss[:, 1] * y) / 4 dx = ss[:, 0] * (1 + ss[:, 1] * y) / (4 * a) dy = ss[:, 1] * (1 + ss[:, 0] * x) / (4 * b) return n, dx, dy def raw_shear(x, y, a, b): n, dx, dy = shell_shape(x, y, a, b) s = np.zeros((2, 24)) for i in range(4): s[0, 6 * i + 2], s[0, 6 * i + 4] = dx[i], n[i] s[1, 6 * i + 2], s[1, 6 * i + 3] = dy[i], -n[i] return s for element in elements: kind = element["type"] if kind not in {"tet4", "hex8", "shell4", "beam2"} or len(element["nodes"]) != {"tet4": 4, "hex8": 8, "shell4": 4, "beam2": 2}[kind]: raise ValueError("Unsupported element type or connectivity.") if len(set(element["nodes"])) != len(element["nodes"]): raise ValueError("An element repeats a node.") ids = [node_index[n] for n in element["nodes"]] touched.update(ids) xyz = coords[ids] dofs = np.array([6 * i + j for i in ids for j in range(6)]) m = materials[element["materialId"]] ke = np.zeros((len(dofs), len(dofs))) if kind in {"hex8", "tet4"}: c = elasticity(m) rotation = strain_rotation(axes(m["axis1"], m["axis2"])) def derivatives(natural): if kind == "tet4": return np.array([[-1, -1, -1], [1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float) x, y, z = natural return np.array([[s[0] * (1 + s[1] * y) * (1 + s[2] * z), s[1] * (1 + s[0] * x) * (1 + s[2] * z), s[2] * (1 + s[0] * x) * (1 + s[1] * y)] for s in signs]) / 8 checks = list(signs) + [[0, 0, 0]] if kind == "hex8" else [[.25, .25, .25]] for natural in checks: det = float(np.linalg.det(xyz.T @ derivatives(natural))) minimum_jacobian = min(minimum_jacobian, det) if det <= 1e-12: raise ValueError("Zero or inverted solid Jacobian.") integration = [([.25, .25, .25], 1 / 6)] if kind == "tet4" else [(s * g, 1) for s in signs] for natural, weight in integration: dn = derivatives(natural) jac = xyz.T @ dn det = float(np.linalg.det(jac)) if det <= 1e-12: raise ValueError("Zero or inverted solid quadrature Jacobian.") minimum_jacobian = min(minimum_jacobian, det) grad = dn @ np.linalg.inv(jac) b = np.zeros((6, len(dofs))) for i, (dx, dy, dz) in enumerate(grad): b[0, 6 * i], b[1, 6 * i + 1], b[2, 6 * i + 2] = dx, dy, dz b[3, 6 * i + 1], b[3, 6 * i + 2] = dz, dy b[4, 6 * i], b[4, 6 * i + 2] = dz, dx b[5, 6 * i], b[5, 6 * i + 1] = dy, dx b = rotation @ b ke += b.T @ c @ b * det * weight elif kind == "shell4": ex, ey = xyz[1] - xyz[0], xyz[3] - xyz[0] r = axes(ex, ey) length, height = np.linalg.norm(ex), np.linalg.norm(ey) if abs(np.dot(ex / length, ey / height)) > 1e-8 or np.linalg.norm(xyz[2] - xyz[0] - ex - ey) > 1e-8 * max(length, height, 1): raise ValueError("Only flat rectangular assumed-shear shell4 panels are supported.") a, b = length / 2, height / 2 minimum_jacobian = min(minimum_jacobian, a * b) transform = np.zeros((24, 24)) for i in range(8): transform[3 * i:3 * i + 3, 3 * i:3 * i + 3] = r plies = element.get("plies") or [{"materialId": element["materialId"], "thicknessMm": element["thicknessMm"], "angleDeg": 0}] if not 1 <= len(plies) <= 16 or any(p["thicknessMm"] <= 0 for p in plies): raise ValueError("Invalid shell laminate thickness.") thickness = sum(p["thicknessMm"] for p in plies) z0 = -thickness / 2 for ply in plies: material = materials[ply["materialId"]] c0 = elasticity(material) c = c0 - np.outer(c0[:, 2], c0[2, :]) / c0[2, 2] c[2, :] = c[:, 2] = 0 c[3, 3] *= 5 / 6 c[4, 4] *= 5 / 6 angle = math.radians(ply["angleDeg"]) cs, sn = math.cos(angle), math.sin(angle) rotation = strain_rotation(np.array([[cs, sn, 0], [-sn, cs, 0], [0, 0, 1]])) for x in [-g, g]: for y in [-g, g]: _, dx, dy = shell_shape(x, y, a, b) shear_x = ((1 - y) * raw_shear(0, -1, a, b)[0] + (1 + y) * raw_shear(0, 1, a, b)[0]) / 2 shear_y = ((1 - x) * raw_shear(-1, 0, a, b)[1] + (1 + x) * raw_shear(1, 0, a, b)[1]) / 2 for zz in [-g, g]: z = z0 + ply["thicknessMm"] * (1 + zz) / 2 strain = np.zeros((6, 24)) for i in range(4): strain[0, 6 * i], strain[0, 6 * i + 4] = dx[i], z * dx[i] strain[1, 6 * i + 1], strain[1, 6 * i + 3] = dy[i], -z * dy[i] strain[5, 6 * i], strain[5, 6 * i + 1] = dy[i], dx[i] strain[5, 6 * i + 4], strain[5, 6 * i + 3] = z * dy[i], -z * dx[i] strain[3, :], strain[4, :] = shear_y, shear_x strain = rotation @ strain @ transform ke += strain.T @ c @ strain * a * b * ply["thicknessMm"] / 2 z0 += ply["thicknessMm"] drill_modulus = min(materials[p["materialId"]]["G12"] for p in plies) * thickness * 1e-4 for x in [-g, g]: for y in [-g, g]: shape, dx, dy = shell_shape(x, y, a, b) drill = np.zeros(24) for i in range(4): drill[6 * i], drill[6 * i + 1], drill[6 * i + 5] = dy[i] / 2, -dx[i] / 2, shape[i] drill = drill @ transform ke += np.outer(drill, drill) * a * b * drill_modulus else: section = element["section"] ex = xyz[1] - xyz[0] length = np.linalg.norm(ex) if length <= 1e-9 or min(section[k] for k in ["areaMm2", "iyMm4", "izMm4", "torsionMm4"]) <= 0: raise ValueError("Invalid beam geometry or section.") r = axes(ex, section["localY"]) transform = np.zeros((12, 12)) for i in range(4): transform[3 * i:3 * i + 3, 3 * i:3 * i + 3] = r local = np.zeros((12, 12)) for i, j, stiffness in [(0, 6, m["E1"] * section["areaMm2"] / length), (3, 9, m["G12"] * section["torsionMm4"] / length)]: local[i, i] += stiffness local[j, j] += stiffness local[i, j] -= stiffness local[j, i] -= stiffness for ix, inertia, sign in [([1, 5, 7, 11], section["izMm4"], 1), ([2, 4, 8, 10], section["iyMm4"], -1)]: block = np.array([[12, sign * 6 * length, -12, sign * 6 * length], [sign * 6 * length, 4 * length**2, -sign * 6 * length, 2 * length**2], [-12, -sign * 6 * length, 12, -sign * 6 * length], [sign * 6 * length, 2 * length**2, -sign * 6 * length, 4 * length**2]]) * m["E1"] * inertia / length**3 local[np.ix_(ix, ix)] += block ke = transform.T @ local @ transform minimum_jacobian = min(minimum_jacobian, length) if not np.all(np.isfinite(ke)) or np.max(np.abs(ke - ke.T)) > 1e-7 * max(1, np.max(np.abs(ke))): raise ValueError("Element matrix is nonfinite or nonsymmetric.") for i, row in enumerate(dofs): for j, col in enumerate(dofs): if ke[i, j] != 0: coo_i.append(int(row)); coo_j.append(int(col)); coo_v.append(float(ke[i, j])) if len(touched) != len(nodes): raise ValueError("Unconnected source nodes are not accepted.") size = len(nodes) * 6 stiffness = coo_matrix((coo_v, (coo_i, coo_j)), shape=(size, size)).tocsr() load = np.zeros(size) def at(entry): if type(entry["dof"]) is not int or not 0 <= entry["dof"] <= 5: raise ValueError("Invalid exact node degree of freedom.") return 6 * node_index[entry["node"]] + entry["dof"] if len({at(q) for q in model["loads"]}) != len(model["loads"]): raise ValueError("Duplicate loaded DOF.") for q in model["loads"]: if abs(q["value"]) > 1e12: raise ValueError("Load exceeds the declared bound.") load[at(q)] = q["value"] equations = {} for q in model["boundary"]: i = at(q) if i in equations or abs(q["value"]) > 1e6: raise ValueError("Duplicate or excessive prescribed DOF.") equations[i] = (q["value"], []) for q in model["ties"]: i = at(q["slave"]) if i in equations or not 1 <= len(q["terms"]) <= 24: raise ValueError("Invalid or conflicting dependent DOF.") masters = [(at(t), t["coefficient"]) for t in q["terms"]] if len({j for j, _ in masters}) != len(masters) or any(j == i or abs(v) > 1e7 for j, v in masters): raise ValueError("Invalid multipoint master.") equations[i] = (q["offset"], masters) diagonal = stiffness.diagonal() unused = np.asarray(stiffness.getnnz(axis=1) == 0) for i in np.where(unused)[0]: if abs(load[i]) > 1e-12 or i in equations and equations[i][0] != 0: raise ValueError("A load or motion acts on an unused solid rotational DOF.") equations.setdefault(int(i), (0, [])) free = [i for i in range(size) if i not in equations] free_index = {n: i for i, n in enumerate(free)} expanded, visiting = {}, set() def expand(i): if i in expanded: return expanded[i] if i in visiting: raise ValueError("Multipoint constraint dependency cycle.") if i in free_index: return 0, {free_index[i]: 1} if len(visiting) > 128: raise ValueError("Multipoint constraint depth exceeds the reviewed bound.") visiting.add(i) offset, masters = equations[i] terms = {} for j, coefficient in masters: value, dependencies = expand(j) offset += coefficient * value for k, factor in dependencies.items(): terms[k] = terms.get(k, 0) + coefficient * factor if len(terms) > 512: raise ValueError("Multipoint expansion exceeds the reviewed sparse bound.") visiting.remove(i) expanded[i] = (offset, terms) return expanded[i] u0 = np.zeros(size) ti, tj, tv = [], [], [] for i in range(size): offset, terms = expand(i) u0[i] = offset for j, value in terms.items(): ti.append(i); tj.append(j); tv.append(value) transform = csr_matrix((tv, (ti, tj)), shape=(size, len(free))) if free: reduced = transform.T @ stiffness @ transform rhs = transform.T @ (load - stiffness @ u0) with warnings.catch_warnings(): warnings.simplefilter("error", MatrixRankWarning) unknowns = spsolve(reduced.tocsc(), rhs) displacement = u0 + transform @ unknowns else: displacement = u0 if not np.all(np.isfinite(displacement)) or np.max(np.abs(displacement)) > 1e9: raise ValueError("Sparse solution is singular or nonfinite.") internal = stiffness @ displacement reaction = internal - load residual = transform.T @ reaction max_residual = float(np.max(np.abs(residual))) if len(free) else 0 scale = max(1, float(np.max(np.abs(load))), float(np.max(np.abs(internal)))) if max_residual > 1e-7 * scale + 1e-7: raise ValueError("Sparse solution fails original free-DOF equilibrium.") energy = float(displacement @ internal / 2) force = internal.reshape((-1, 6))[:, :3].sum(axis=0) moment = (np.cross(coords, internal.reshape((-1, 6))[:, :3]) + internal.reshape((-1, 6))[:, 3:]).sum(axis=0) return {"nodes": [{"id": n["id"], "displacement": displacement[6 * i:6 * i + 6].tolist()} for i, n in enumerate(nodes)], "checks": {"nodes": len(nodes), "elements": len(elements), "activeDofs": int(np.sum(~unused)), "freeDofs": len(free), "matrixNonzeros": int(stiffness.nnz), "maxFreeResidual": max_residual, "strainEnergyNmm": energy, "minimumJacobian": minimum_jacobian, "forceBalanceN": force.tolist(), "momentBalanceNmm": moment.tolist()}, "solverVersion": {"worker": VERSION, "numpy": np.__version__, "scipy": scipy.__version__}, "scope": "Small-strain linear tet4/hex8, flat rectangular assumed-shear laminate shell4 and Euler-Bernoulli beam2; exact prescribed/MPC elimination; SciPy sparse solve. No external solver or physical validation is claimed."} def write_private(path, value): encoded = json.dumps(value, separators=(",", ":"), allow_nan=False).encode("utf-8") if len(encoded) > MAX_OUTPUT: raise ValueError("Receipt exceeds the 24 MiB standalone output bound.") target = Path(path).resolve() target.parent.mkdir(parents=True, exist_ok=True) if target.exists(): raise ValueError("Choose a new receipt filename; the existing output is preserved.") handle, temporary = tempfile.mkstemp(prefix=".lingo-fe28-", dir=target.parent) try: os.fchmod(handle, 0o600) with os.fdopen(handle, "wb") as stream: stream.write(encoded) stream.flush() os.fsync(stream.fileno()) os.link(temporary, target) finally: if os.path.exists(temporary): os.unlink(temporary) def main(): parser = argparse.ArgumentParser(description="Run an exact reviewed LignoSparseFEJob28 locally with NumPy/SciPy; no network or device control.") parser.add_argument("--input", required=True, help="Prepared job JSON path, or - for stdin") parser.add_argument("--output", required=True, help="Private source-bound receipt JSON path") args = parser.parse_args() if args.input == "-": raw = sys.stdin.buffer.read(MAX_JOB + 1) else: with Path(args.input).open("rb") as stream: raw = stream.read(MAX_JOB + 1) if len(raw) > MAX_JOB: raise ValueError("Job transport exceeds the reviewed bound.") job = strict_json(raw.decode("utf-8")) if set(job) != {"format", "version", "inputJSON", "inputSHA256", "workerSourceSHA256"} or job["format"] != "LignoSparseFEJob28" or job["version"] != 1: raise ValueError("Unsupported prepared worker job.") if not isinstance(job["inputJSON"], str) or len(job["inputJSON"].encode("utf-8")) > MAX_INPUT: raise ValueError("Original source exceeds the 2 MiB bound.") source_hash = digest(Path(__file__).read_bytes()) if job["workerSourceSHA256"] != source_hash or digest(job["inputJSON"]) != job["inputSHA256"]: raise ValueError("Worker or exact original input fingerprint differs from the reviewed job.") model = strict_json(job["inputJSON"]) finite_tree(model) result = run_model(model) result_json = json.dumps(result, separators=(",", ":"), allow_nan=False) log = "Completed fixed-schema original elastic FE source with NumPy/SciPy sparse assembly and exact constraint elimination." receipt = {"format": "LignoSparseFEReceipt28", "version": 1, "inputSHA256": job["inputSHA256"], "workerSourceSHA256": source_hash, "outputJSON": result_json, "outputSHA256": digest(result_json), "log": log, "logSHA256": digest(log)} write_private(args.output, receipt) print("Completed original elastic finite-element job; private receipt written.") if __name__ == "__main__": try: main() except Exception as error: print("Finite-element worker withheld a result: " + str(error), file=sys.stderr) sys.exit(2)