28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 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
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 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 quad indices for the JS-side preview and connectivity checks.
8# See meshio.ts.
10import json
11import os
13import numpy as np
15import meshio
17WORK = "/work"
18OUT_MSH = WORK + "/out.msh"
19POSITIONS_F32 = WORK + "/positions.f32"
20QUADS_U32 = WORK + "/quads.u32"
22os.makedirs(WORK, exist_ok=True)
25def _collect_quads(mesh, warnings):
26 """All quad cells, corner-nodes only; reject meshes without any."""
27 blocks = []
28 found = set()
29 for block in mesh.cells:
30 data = block.data
31 if not isinstance(data, np.ndarray) or data.ndim != 2:
32 continue
33 found.add(block.type)
34 if block.type == "quad":
35 blocks.append(data)
36 elif block.type in ("quad8", "quad9"):
37 blocks.append(data[:, :4])
38 warnings.append(
39 f"{len(data)} higher-order {block.type} cells reduced to corner nodes"
40 )
41 elif block.type == "polygon" and data.shape[1] == 4:
42 blocks.append(data)
43 if not blocks:
44 kinds = ", ".join(sorted(found)) or "none"
45 raise ValueError(
46 "No quadrilateral cells found (cell types in file: "
47 + kinds
48 + "). surfacefun solves on quad meshes; "
49 "convert your mesh to quads before uploading."
50 )
51 return np.ascontiguousarray(np.vstack(blocks).astype(np.int64))
54def _write_msh(path, points, quads):
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 55 """Canonical Gmsh MSH 4.1 ASCII: one surface entity block, sequential
56 1-based node ids, 4-node quads (type 3) — what surfacemesh.import reads."""
57 n, m = len(points), len(quads)
58 lines = ["$MeshFormat", "4.1 0 8", "$EndMeshFormat"]
59 lines.append("$Nodes")
60 lines.append("1 %d 1 %d" % (n, n)) # numEntityBlocks numNodes minTag maxTag
61 lines.append("2 1 0 %d" % n) # entityDim entityTag parametric numNodes
62 lines.extend(str(i) for i in range(1, n + 1))
63 lines.extend("%.16g %.16g %.16g" % (p[0], p[1], p[2]) for p in points)
28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 64 lines.append("$EndNodes")
65 lines.append("$Elements")
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 66 lines.append("1 %d 1 %d" % (m, m)) # numEntityBlocks numElements minTag maxTag
67 lines.append("2 1 3 %d" % m) # entityDim entityTag elementType(3=quad) numElements
28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 68 for i, q in enumerate(quads, start=1):
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 69 lines.append("%d %d %d %d %d" % (i, q[0] + 1, q[1] + 1, q[2] + 1, q[3] + 1))
28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 70 lines.append("$EndElements")
71 with open(path, "w") as f:
72 f.write("\n".join(lines) + "\n")
75def parse_quad_mesh(path, file_format=None):
76 warnings = []
77 try:
78 # meshio's read helper exits the interpreter when every candidate
79 # reader fails; turn that into a normal exception
80 mesh = meshio.read(path, file_format)
81 except SystemExit:
82 raise ValueError(f"Could not read file as {file_format or 'any known format'}")
84 points = np.asarray(mesh.points, dtype=np.float64)
85 if points.ndim != 2:
86 raise ValueError(f"Unexpected points array shape {points.shape}")
87 if points.shape[1] == 2:
88 points = np.column_stack([points, np.zeros(len(points))])
89 warnings.append("2D points: added z=0")
90 points = np.ascontiguousarray(points[:, :3])
92 quads = _collect_quads(mesh, warnings)
93 if quads.size and (int(quads.min()) < 0 or int(quads.max()) >= len(points)):
94 raise ValueError("Quad node index out of range")
96 # Drop vertices not referenced by any quad (e.g. triangle-only regions of
97 # a mixed mesh) so the .msh stays minimal and ids stay dense.
98 used = np.unique(quads)
99 if len(used) < len(points):
100 remap = np.full(len(points), -1, dtype=np.int64)
101 remap[used] = np.arange(len(used))
102 points = points[used]
103 quads = remap[quads]
104 warnings.append(f"dropped {len(remap) - len(used)} unused vertices")
106 _write_msh(OUT_MSH, points, quads)
107 with open(POSITIONS_F32, "wb") as f:
108 f.write(points.astype(np.float32).tobytes())
109 with open(QUADS_U32, "wb") as f:
110 f.write(np.ascontiguousarray(quads.astype(np.uint32)).tobytes())
112 return json.dumps(
113 {
114 "numVertices": len(points),
115 "numQuads": len(quads),
116 "warnings": warnings,
117 }
118 )