1# Runs inside Pyodide. Bridges meshio to the JS app.
2#
3# parse_quad_mesh reads an uploaded mesh file with meshio, keeps only its
4# quadrilateral cells, and produces three outputs in Pyodide's in-memory
5# filesystem: the canonical Gmsh MSH 2.2 ASCII file the numbl solver reads
6# (full float64 precision), plus float32 positions and uint32 quad indices
7# for the JS-side preview and connectivity checks. See meshio.ts.
9import json
10import os
12import numpy as np
14import meshio
16WORK = "/work"
17OUT_MSH = WORK + "/out.msh"
18POSITIONS_F32 = WORK + "/positions.f32"
19QUADS_U32 = WORK + "/quads.u32"
21os.makedirs(WORK, exist_ok=True)
24def _collect_quads(mesh, warnings):
25 """All quad cells, corner-nodes only; reject meshes without any."""
26 blocks = []
27 found = set()
28 for block in mesh.cells:
29 data = block.data
30 if not isinstance(data, np.ndarray) or data.ndim != 2:
31 continue
32 found.add(block.type)
33 if block.type == "quad":
34 blocks.append(data)
35 elif block.type in ("quad8", "quad9"):
36 blocks.append(data[:, :4])
37 warnings.append(
38 f"{len(data)} higher-order {block.type} cells reduced to corner nodes"
39 )
40 elif block.type == "polygon" and data.shape[1] == 4:
41 blocks.append(data)
42 if not blocks:
43 kinds = ", ".join(sorted(found)) or "none"
44 raise ValueError(
45 "No quadrilateral cells found (cell types in file: "
46 + kinds
47 + "). surfacefun solves on quad meshes; "
48 "convert your mesh to quads before uploading."
49 )
50 return np.ascontiguousarray(np.vstack(blocks).astype(np.int64))
53def _write_msh(path, points, quads):
54 """Canonical Gmsh MSH 2.2 ASCII: sequential ids, type-3 quads, two tags."""
55 lines = ["$MeshFormat", "2.2 0 8", "$EndMeshFormat", "$Nodes", str(len(points))]
56 for i, p in enumerate(points, start=1):
57 lines.append("%d %.16g %.16g %.16g" % (i, p[0], p[1], p[2]))
58 lines.append("$EndNodes")
59 lines.append("$Elements")
60 lines.append(str(len(quads)))
61 for i, q in enumerate(quads, start=1):
62 lines.append("%d 3 2 1 1 %d %d %d %d" % (i, q[0] + 1, q[1] + 1, q[2] + 1, q[3] + 1))
63 lines.append("$EndElements")
64 with open(path, "w") as f:
65 f.write("\n".join(lines) + "\n")
68def parse_quad_mesh(path, file_format=None):
69 warnings = []
70 try:
71 # meshio's read helper exits the interpreter when every candidate
72 # reader fails; turn that into a normal exception
73 mesh = meshio.read(path, file_format)
74 except SystemExit:
75 raise ValueError(f"Could not read file as {file_format or 'any known format'}")
77 points = np.asarray(mesh.points, dtype=np.float64)
78 if points.ndim != 2:
79 raise ValueError(f"Unexpected points array shape {points.shape}")
80 if points.shape[1] == 2:
81 points = np.column_stack([points, np.zeros(len(points))])
82 warnings.append("2D points: added z=0")
83 points = np.ascontiguousarray(points[:, :3])
85 quads = _collect_quads(mesh, warnings)
86 if quads.size and (int(quads.min()) < 0 or int(quads.max()) >= len(points)):
87 raise ValueError("Quad node index out of range")
89 # Drop vertices not referenced by any quad (e.g. triangle-only regions of
90 # a mixed mesh) so the .msh stays minimal and ids stay dense.
91 used = np.unique(quads)
92 if len(used) < len(points):
93 remap = np.full(len(points), -1, dtype=np.int64)
94 remap[used] = np.arange(len(used))
95 points = points[used]
96 quads = remap[quads]
97 warnings.append(f"dropped {len(remap) - len(used)} unused vertices")
99 _write_msh(OUT_MSH, points, quads)
100 with open(POSITIONS_F32, "wb") as f:
101 f.write(points.astype(np.float32).tobytes())
102 with open(QUADS_U32, "wb") as f:
103 f.write(np.ascontiguousarray(quads.astype(np.uint32)).tobytes())
105 return json.dumps(
106 {
107 "numVertices": len(points),
108 "numQuads": len(quads),
109 "warnings": warnings,
110 }
111 )