"""Generate a CAD-conforming display volume. This is not a CFD solve."""
import sys,json,time
from pathlib import Path
import gmsh,numpy as np
p=Path(sys.argv[1]);start=time.perf_counter()
gmsh.initialize();gmsh.option.setNumber('General.Terminal',0)
try:
 gmsh.model.occ.importShapes(str(p/'input.step'));gmsh.model.occ.synchronize()
 box=np.array(gmsh.model.getBoundingBox(-1,-1));extent=box[3:]-box[:3];h=float(max(extent)/28)
 gmsh.option.setNumber('Mesh.MeshSizeMax',h);gmsh.option.setNumber('Mesh.MeshSizeMin',h/3)
 gmsh.option.setNumber('Mesh.MeshSizeFromCurvature',12);gmsh.option.setNumber('Mesh.ElementOrder',1)
 gmsh.model.mesh.generate(3)
 ids,coords,_=gmsh.model.mesh.getNodes();v=np.asarray(coords).reshape(-1,3);lookup={int(n):i for i,n in enumerate(ids)}
 def elements(dim,typ):
  kinds,_,nodes=gmsh.model.mesh.getElements(dim);x=nodes[list(kinds).index(typ)];return np.array([lookup[int(i)] for i in x]).reshape(-1,dim+1)
 tet=elements(3,4);faces=elements(2,2)
 reference=json.loads((p/'geometry.json').read_text());actual=np.ptp(v,axis=0)
 if not np.allclose(actual,reference['dimensions_mm'],rtol=.02,atol=.3):raise ValueError('Display mesh and STEP reader disagree on dimensions or units')
 corners=v[tet];volume=float(np.abs(np.linalg.det(corners[:,1:]-corners[:,:1])).sum()/6);error=abs(volume-reference['volume_mm3'])/reference['volume_mm3']
 if error>.1:raise ValueError('Display mesh volume differs from the CAD solid by over 10%; a finer mesh is required')
 if len(v)>80000 or len(tet)>400000:raise ValueError('Display volume exceeds memory budget')
 np.savez_compressed(p/'volume.npz',vertices=v,tetrahedra=tet,triangles=faces)
 (p/'volume.json').write_text(json.dumps({'nodes':len(v),'tetrahedra':len(tet),'surface_triangles':len(faces),'method':'CAD-conforming first-order Gmsh tetrahedra, display only','volume_error_fraction':error,'seconds':time.perf_counter()-start}))
finally:gmsh.finalize()