#!/usr/bin/env python3 """Inspect exact STEP geometry with installed OpenCASCADE; no network or solver calls.""" import argparse import hashlib import json import math from pathlib import Path import re import sys import tempfile VERSION = 'LignoCAD28-1' MAX_BYTES = 1024 * 1024 def digest(value): return hashlib.sha256(value).hexdigest() def inspect_job(raw): if len(raw) > 1500000: raise ValueError('The complete CAD job exceeds 1.5 MB.') def no_duplicates(pairs): result = {} for key, value in pairs: if key in result: raise ValueError('Duplicate JSON keys are not accepted.') result[key] = value return result job = json.loads(raw.decode('utf-8'), object_pairs_hook=no_duplicates, parse_constant=lambda value: (_ for _ in ()).throw(ValueError('Non-finite JSON.'))) expected = {'format', 'version', 'originalName', 'sourceSTEP', 'sourceSHA256', 'repair', 'deflectionMm', 'maxTriangles', 'volumeMeshing', 'meshSizeMm', 'elementOrder'} if not isinstance(job, dict) or set(job) != expected or job['format'] != 'LignoCADJob28' or job['version'] != 1: raise ValueError('Use the exact supported CAD job format.') source = job['sourceSTEP'] if not isinstance(source, str) or len(source.encode('utf-8')) > MAX_BYTES or 'ISO-10303-21;' not in source: raise ValueError('Provide an original STEP Part21 source up to 1 MiB.') source_bytes = source.encode('utf-8') if digest(source_bytes) != job['sourceSHA256']: raise ValueError('The original STEP source digest changed.') if re.search(r'\b(EXTERNAL_SOURCE|EXTERNALLY_DEFINED_ITEM|DOCUMENT_FILE|EXTERNALLY_DEFINED_REPRESENTATION)\s*\(', source, re.I): raise ValueError('Resolve and flatten external document references before this offline inspection.') deflection = job['deflectionMm'] limit = job['maxTriangles'] if isinstance(deflection, bool) or not isinstance(deflection, (int, float)) or not math.isfinite(deflection) or not .01 <= deflection <= 100: raise ValueError('Meshing deflection must be between 0.01 and 100 mm.') if isinstance(limit, bool) or not isinstance(limit, int) or not 12 <= limit <= 6000 or not isinstance(job['repair'], bool): raise ValueError('Use a bounded triangle limit and explicit repair flag.') if not isinstance(job['volumeMeshing'], bool) or type(job['elementOrder']) is not int or job['elementOrder'] not in (1, 2) or isinstance(job['meshSizeMm'], bool) or not isinstance(job['meshSizeMm'], (int, float)) or not math.isfinite(job['meshSizeMm']) or not .1 <= job['meshSizeMm'] <= 10000: raise ValueError('Review the bounded volume mesh size and element order.') if not isinstance(job['originalName'], str) or not 1 <= len(job['originalName']) <= 180: raise ValueError('Provide a bounded original filename.') from OCP.STEPControl import STEPControl_Reader from OCP.IFSelect import IFSelect_RetDone from OCP.Interface import Interface_Static from OCP.BRepCheck import BRepCheck_Analyzer from OCP.ShapeFix import ShapeFix_Shape from OCP.BRepGProp import BRepGProp from OCP.GProp import GProp_GProps from OCP.Bnd import Bnd_Box from OCP.BRepBndLib import BRepBndLib from OCP.BRepMesh import BRepMesh_IncrementalMesh from OCP.TopExp import TopExp_Explorer from OCP.TopAbs import TopAbs_SOLID, TopAbs_FACE, TopAbs_REVERSED from OCP.TopoDS import TopoDS from OCP.BRep import BRep_Tool from OCP.TopLoc import TopLoc_Location Interface_Static.SetCVal_s('xstep.cascade.unit', 'MM') with tempfile.TemporaryDirectory(prefix='lingo-cad28-') as directory: original = Path(directory) / 'original.step' original.write_bytes(source_bytes) reader = STEPControl_Reader() if reader.ReadFile(str(original)) != IFSelect_RetDone or reader.TransferRoots() < 1: raise ValueError('OpenCASCADE could not transfer the original STEP roots.') shape = reader.OneShape() if shape.IsNull(): raise ValueError('The transferred CAD shape is empty.') before = bool(BRepCheck_Analyzer(shape).IsValid()) repaired = False if job['repair'] and not before: fixer = ShapeFix_Shape(shape) fixer.Perform() shape = fixer.Shape() repaired = True after = bool(BRepCheck_Analyzer(shape).IsValid()) volume = GProp_GProps() area = GProp_GProps() BRepGProp.VolumeProperties_s(shape, volume) BRepGProp.SurfaceProperties_s(shape, area) box = Bnd_Box() BRepBndLib.Add_s(shape, box) bounds = list(box.Get()) def count(kind): explorer = TopExp_Explorer(shape, kind) total = 0 while explorer.More(): total += 1 explorer.Next() return total solids, faces = count(TopAbs_SOLID), count(TopAbs_FACE) if not after or solids < 1 or volume.Mass() <= 0 or faces > 12000: raise ValueError('The CAD source has no valid positive-volume bounded solid after inspection.') mesher = BRepMesh_IncrementalMesh(shape, float(deflection), False, .35, False) mesher.Perform() if not mesher.IsDone(): raise ValueError('The CAD surface mesher did not complete.') triangles = [] explorer = TopExp_Explorer(shape, TopAbs_FACE) while explorer.More(): face = TopoDS.Face_s(explorer.Current()) location = TopLoc_Location() triangulation = BRep_Tool.Triangulation_s(face, location) if triangulation is None: raise ValueError('A CAD face has no surface triangulation.') transform = location.Transformation() for index in range(1, triangulation.NbTriangles() + 1): if len(triangles) >= limit: raise ValueError('The complete CAD preview exceeds the triangle budget; increase deflection and reissue the reviewed job.') ids = list(triangulation.Triangle(index).Get()) if face.Orientation() == TopAbs_REVERSED: ids[1], ids[2] = ids[2], ids[1] points = [] for node_id in ids: point = triangulation.Node(node_id).Transformed(transform) points.append([point.X(), point.Y(), point.Z()]) triangles.append(points) explorer.Next() values = [volume.Mass(), area.Mass(), *bounds, *(v for triangle in triangles for point in triangle for v in point)] if not all(math.isfinite(v) for v in values): raise ValueError('The CAD calculation returned non-finite geometry.') volume_mesh = generate_volume_mesh(original, job, volume.Mass()) if job['volumeMeshing'] else None return { 'format': 'LignoCADReceipt28', 'jobSHA256': digest(raw), 'sourceSHA256': digest(source_bytes), 'workerSHA256': digest(Path(__file__).read_bytes()), 'workerVersion': VERSION, 'units': 'mm', 'validBefore': before, 'validAfter': after, 'repaired': repaired, 'solids': solids, 'faces': faces, 'volumeMm3': volume.Mass(), 'areaMm2': area.Mass(), 'boundsMm': bounds, 'triangles': triangles, 'volumeMesh': volume_mesh, 'diagnostics': [ 'OpenCASCADE transferred the original STEP source into millimetres.', 'Bounding boxes include the CAD model tolerance; preview triangles approximate the original boundary.', 'Shape healing is attempted only when validity fails; the original source bytes are preserved.', 'Intersections between independent solids, semantic materials and FE boundary conditions require separate review.', 'A requested Gmsh volume mesh uses unassigned CAD material identities; review materials, supports, loads and interfaces before FE analysis.', 'This receipt establishes portable source consistency, not signed execution attestation or calibrated timber validity.' ] } def generate_volume_mesh(original, job, exact_volume): import gmsh import numpy as np if gmsh.isInitialized(): raise ValueError('Run the reviewed CAD worker in an isolated process; an existing Gmsh session is refused.') gmsh.initialize() try: gmsh.option.setNumber('General.Terminal', 0) gmsh.option.setNumber('General.NumThreads', 1) gmsh.option.setNumber('Mesh.RandomSeed', 1) gmsh.option.setString('Geometry.OCCTargetUnit', 'MM') gmsh.model.add('exact-original') gmsh.model.occ.importShapes(str(original), highestDimOnly=True) gmsh.model.occ.synchronize() if not gmsh.model.getEntities(3): raise ValueError('The original CAD has no meshable volume.') volumes = gmsh.model.getEntities(3) if len(volumes) > 1: gmsh.model.occ.fragment(volumes[:1], volumes[1:]) gmsh.model.occ.synchronize() gmsh.option.setNumber('Mesh.MinimumCircleNodes', 16) gmsh.option.setNumber('Mesh.MeshSizeMin', job['meshSizeMm']) gmsh.option.setNumber('Mesh.MeshSizeMax', job['meshSizeMm']) gmsh.model.mesh.generate(3) if job['elementOrder'] == 2: gmsh.model.mesh.setOrder(2) gmsh.model.mesh.optimize('HighOrder') tags, flat, _ = gmsh.model.mesh.getNodes() if len(tags) > 2200: raise ValueError('The complete CAD mesh exceeds 2,200 nodes; coarsen and issue a new reviewed job.') nodes = [{'id': 'cad-n' + str(int(tag)), 'xyz': [float(q) for q in xyz]} for tag, xyz in zip(tags, np.asarray(flat).reshape(-1, 3))] by_tag = {int(tag): np.asarray(xyz, dtype=float) for tag, xyz in zip(tags, np.asarray(flat).reshape(-1, 3))} owners = {} for _, volume_tag in gmsh.model.getEntities(3): _, groups, _ = gmsh.model.mesh.getElements(3, volume_tag) for group in groups: owners.update({int(tag):volume_tag for tag in group}) element_types, element_tags, node_tags = gmsh.model.mesh.getElements(3) if sum(len(group) for group in element_tags) > 2400: raise ValueError('The complete CAD mesh exceeds 2,400 tetrahedra; coarsen and issue a new reviewed job.') elements, all_tags = [], [] integrated_volume = 0.0 minimum_jacobian = float('inf') for element_type, group_tags, group_nodes in zip(element_types, element_tags, node_tags): if int(element_type) not in (4, 11): raise ValueError('Only actual tetra4/tetra10 volume cells are supported by this CAD receipt.') size = 4 if int(element_type) == 4 else 10 for tag, cell in zip(group_tags, np.asarray(group_nodes).reshape(-1, size)): ids = [int(v) for v in cell] # Gmsh tetra10 edges end with (2,3),(1,3); the platform uses (1,3),(2,3). if size == 10: ids[8], ids[9] = ids[9], ids[8] xyz = np.asarray([by_tag[v] for v in ids]) def jacobian(r, t, u): if size == 4: return np.column_stack((xyz[1]-xyz[0], xyz[2]-xyz[0], xyz[3]-xyz[0])) L = np.asarray([1-r-t-u, r, t, u]) dL = np.asarray([[-1,-1,-1],[1,0,0],[0,1,0],[0,0,1]]) dN = [(4*L[i]-1)*dL[i] for i in range(4)] for a, b in [(0,1),(1,2),(2,0),(0,3),(1,3),(2,3)]: dN.append(4*(dL[a]*L[b]+L[a]*dL[b])) return xyz.T @ np.asarray(dN) if size == 4: cell_volume = float(np.linalg.det(jacobian(.25,.25,.25)))/6 else: cell_volume = -2/15*float(np.linalg.det(jacobian(.25,.25,.25))) for q in [(.5,1/6,1/6),(1/6,.5,1/6),(1/6,1/6,.5),(1/6,1/6,1/6)]: cell_volume += 3/40*float(np.linalg.det(jacobian(*q))) minimum_jacobian = min(minimum_jacobian, *(float(np.linalg.det(jacobian(*q))) for q in [(.25,.25,.25),(.5,1/6,1/6),(1/6,.5,1/6),(1/6,1/6,.5),(1/6,1/6,1/6)])) if not math.isfinite(cell_volume) or cell_volume <= 0: raise ValueError('The generated CAD mesh contains an invalid volume cell.') integrated_volume += cell_volume elements.append({'id':'cad-e'+str(int(tag)), 'type':'tet4' if size == 4 else 'tet10', 'nodes':['cad-n'+str(v) for v in ids], 'materialId':'cad-material', 'plyId':'cad-unassigned', 'partId':'cad-v'+str(owners[int(tag)]), 'localBounds':[float(np.min(xyz[:,0])),float(np.max(xyz[:,0])),float(np.min(xyz[:,1])),float(np.max(xyz[:,1])),float(np.min(xyz[:,2])),float(np.max(xyz[:,2]))], 'volumeMm3':cell_volume}) all_tags.append(int(tag)) if not elements: raise ValueError('The volume mesher produced no supported complete elements.') gmsh_minimum = float(np.min(gmsh.model.mesh.getElementQualities(all_tags, 'minDetJac'))) min_scaled = float(np.min(gmsh.model.mesh.getElementQualities(all_tags, 'minSJ'))) if gmsh_minimum <= 0 or minimum_jacobian <= 0 or min_scaled <= 0: raise ValueError('The actual CAD volume mesh failed its Jacobian checks.') return {'engine':'Gmsh '+gmsh.__version__, 'nodes':nodes, 'elements':elements, 'minimumJacobianMm3':min(minimum_jacobian,gmsh_minimum), 'minimumScaledJacobian':min_scaled, 'volumeMm3':integrated_volume, 'relativeVolumeError':abs(integrated_volume-exact_volume)/exact_volume, 'unassignedMaterials':True} finally: gmsh.finalize() def envelope_result(result): def values(value): if isinstance(value, dict): return 1 + sum(values(v) for v in value.values()) if isinstance(value, list): return 1 + sum(values(v) for v in value) return 1 if values(result) > 95000: raise ValueError('The complete CAD payload exceeds its retained JSON-value budget; coarsen the mesh/preview and reissue the job.') payload = json.dumps(result, ensure_ascii=False, separators=(',', ':'), allow_nan=False) if len(payload.encode('utf-8')) > 1800000: raise ValueError('The complete CAD payload exceeds 1.8 MB; coarsen and reissue the job.') return {'format':'LignoCADReceiptEnvelope28', 'payloadJSON':payload, 'outputSHA256':digest(payload.encode('utf-8'))} def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--job', required=True, help='Exact downloaded CAD job JSON.') parser.add_argument('--output', required=True, help='New receipt filename; existing files are refused.') args = parser.parse_args() target = Path(args.output) if target.exists(): raise ValueError('Choose a new receipt filename; existing output is preserved.') result = inspect_job(Path(args.job).read_bytes()) envelope = envelope_result(result) with target.open('x', encoding='utf-8') as output: json.dump(envelope, output, ensure_ascii=False, separators=(',', ':'), allow_nan=False) output.write('\n') print(json.dumps({'status': 'completed', 'solids': result['solids'], 'faces': result['faces'], 'triangles': len(result['triangles'])})) if __name__ == '__main__': try: main() except Exception as error: print('CAD inspection failed: ' + str(error), file=sys.stderr) sys.exit(1)