1# Runs inside Pyodide. Bridges meshio to the JS app.
2#
3# parse_mesh reads an uploaded mesh file with meshio, keeps its triangle and
4# quadrilateral cells, and produces three outputs in Pyodide's in-memory
5# filesystem: the canonical Gmsh MSH 4.1 ASCII file the numbl solver reads
6# with surfacemesh.import (full float64 precision), plus float32 positions
7# and uint32 cell indices for the JS-side preview and connectivity checks.
8# surfacefun cannot mix patch types, so a mesh containing both kinds has its
9# quads split into triangles; the output is always homogeneous (all cells
10# 3 nodes or all 4). See meshio.ts.
12import json
13import os
15import numpy as np
17import meshio
19WORK = "/work"
20OUT_MSH = WORK + "/out.msh"
21POSITIONS_F32 = WORK + "/positions.f32"
22CELLS_U32 = WORK + "/cells.u32"
24os.makedirs(WORK, exist_ok=True)
27class _ObjTolerantMesh(meshio.Mesh):
28 """OBJ faces index texture coordinates and normals independently of
29 vertex positions, so a file with UV seams legally has more vt (or vn)
30 entries than v entries. meshio shoehorns those into point_data, whose
31 per-vertex length check then rejects the whole file; drop the unmappable
32 arrays instead and remember what was dropped so callers can warn."""
34 def __init__(self, points, cells, point_data=None, **kwargs):
35 point_data = point_data or {}
36 self.dropped_point_data = {
37 key: len(value)
38 for key, value in point_data.items()
39 if len(value) != len(points)
40 }
41 point_data = {
42 key: value
43 for key, value in point_data.items()
44 if key not in self.dropped_point_data
45 }
46 super().__init__(points, cells, point_data=point_data, **kwargs)
49# the reader binds Mesh at module level, so this rebinding scopes the
50# tolerance to OBJ reads only (elsewhere a mismatch means real corruption)
51meshio.obj._obj.Mesh = _ObjTolerantMesh
53_OBJ_POINT_DATA_NAMES = {"obj:vt": "texture coordinates", "obj:vn": "vertex normals"}
56def _collect_cells(mesh, warnings):
57 """Triangle and quad cells, corner-nodes only; mixed meshes are reduced
58 to all-triangle; reject meshes with neither kind."""
59 tris = []
60 quads = []
61 found = set()
62 for block in mesh.cells:
63 data = block.data
64 if not isinstance(data, np.ndarray) or data.ndim != 2:
65 continue
66 found.add(block.type)
67 if block.type == "triangle":
68 tris.append(data)
69 elif block.type in ("triangle6", "triangle7"):
70 tris.append(data[:, :3])
71 warnings.append(
72 f"{len(data)} higher-order {block.type} cells reduced to corner nodes"
73 )
74 elif block.type == "quad":
75 quads.append(data)
76 elif block.type in ("quad8", "quad9"):
77 quads.append(data[:, :4])
78 warnings.append(
79 f"{len(data)} higher-order {block.type} cells reduced to corner nodes"
80 )
81 elif block.type == "polygon" and data.shape[1] == 3:
82 tris.append(data)
83 elif block.type == "polygon" and data.shape[1] == 4:
84 quads.append(data)
85 if not tris and not quads:
86 kinds = ", ".join(sorted(found)) or "none"
87 raise ValueError(
88 "No triangle or quadrilateral cells found (cell types in file: "
89 + kinds
90 + "). surfacefun solves on triangle or quad meshes."
91 )
92 if tris and quads:
93 nq = sum(len(q) for q in quads)
94 for q in quads:
95 tris.append(q[:, [0, 1, 2]])
96 tris.append(q[:, [0, 2, 3]])
97 quads = []
98 warnings.append(
99 f"{nq} quads split into triangles (surfacefun cannot mix cell types)"
100 )
101 blocks = tris or quads
102 return np.ascontiguousarray(np.vstack(blocks).astype(np.int64))
105def _write_msh(path, points, cells):
106 """Canonical Gmsh MSH 4.1 ASCII: one surface entity block, sequential
107 1-based node ids, 3-node triangles (type 2) or 4-node quads (type 3) —
108 what surfacemesh.import reads."""
109 n, m = len(points), len(cells)
110 etype = 2 if cells.shape[1] == 3 else 3
111 lines = ["$MeshFormat", "4.1 0 8", "$EndMeshFormat"]
112 lines.append("$Nodes")
113 lines.append("1 %d 1 %d" % (n, n)) # numEntityBlocks numNodes minTag maxTag
114 lines.append("2 1 0 %d" % n) # entityDim entityTag parametric numNodes
115 lines.extend(str(i) for i in range(1, n + 1))
116 lines.extend("%.16g %.16g %.16g" % (p[0], p[1], p[2]) for p in points)
117 lines.append("$EndNodes")
118 lines.append("$Elements")
119 lines.append("1 %d 1 %d" % (m, m)) # numEntityBlocks numElements minTag maxTag
120 lines.append("2 1 %d %d" % (etype, m)) # entityDim entityTag elementType numElements
121 for i, c in enumerate(cells, start=1):
122 lines.append(str(i) + " " + " ".join(str(v + 1) for v in c))
123 lines.append("$EndElements")
124 with open(path, "w") as f:
125 f.write("\n".join(lines) + "\n")
128def parse_mesh(path, file_format=None):
129 warnings = []
130 try:
131 # meshio's read helper exits the interpreter when every candidate
132 # reader fails; turn that into a normal exception
133 mesh = meshio.read(path, file_format)
134 except SystemExit:
135 raise ValueError(f"Could not read file as {file_format or 'any known format'}")
137 for key, count in getattr(mesh, "dropped_point_data", {}).items():
138 name = _OBJ_POINT_DATA_NAMES.get(key, key)
139 warnings.append(
140 f"dropped {name}: {count} entries for {len(mesh.points)} vertices"
141 )
143 points = np.asarray(mesh.points, dtype=np.float64)
144 if points.ndim != 2:
145 raise ValueError(f"Unexpected points array shape {points.shape}")
146 if points.shape[1] == 2:
147 points = np.column_stack([points, np.zeros(len(points))])
148 warnings.append("2D points: added z=0")
149 points = np.ascontiguousarray(points[:, :3])
151 cells = _collect_cells(mesh, warnings)
152 if cells.size and (int(cells.min()) < 0 or int(cells.max()) >= len(points)):
153 raise ValueError("Cell node index out of range")
155 # Drop vertices not referenced by any cell (e.g. stray line elements)
156 # so the .msh stays minimal and ids stay dense.
157 used = np.unique(cells)
158 if len(used) < len(points):
159 remap = np.full(len(points), -1, dtype=np.int64)
160 remap[used] = np.arange(len(used))
161 points = points[used]
162 cells = remap[cells]
163 warnings.append(f"dropped {len(remap) - len(used)} unused vertices")
165 _write_msh(OUT_MSH, points, cells)
166 with open(POSITIONS_F32, "wb") as f:
167 f.write(points.astype(np.float32).tobytes())
168 with open(CELLS_U32, "wb") as f:
169 f.write(np.ascontiguousarray(cells.astype(np.uint32)).tobytes())
171 return json.dumps(
172 {
173 "numVertices": len(points),
174 "numCells": len(cells),
175 "cellSize": int(cells.shape[1]),
176 "warnings": warnings,
177 }
178 )