28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 1"""Generate the bundled sample quad meshes (public/samples/*.msh).
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 3Writes Gmsh MSH 4.1 ASCII files in the same canonical form the in-app
4converter produces: one surface entity block, sequential 1-based node ids,
5and 4-node quadrangle elements (type 3) — the layout surfacemesh.import
6reads.
8Run from the repo root: python3 scripts/make_samples.py
9"""
11import math
12import os
14OUT_DIR = os.path.join(os.path.dirname(__file__), "..", "public", "samples")
17def write_msh(path, points, quads):
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 18 n, m = len(points), len(quads)
19 lines = ["$MeshFormat", "4.1 0 8", "$EndMeshFormat"]
28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 20 lines.append("$Nodes")
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 21 lines.append(f"1 {n} 1 {n}") # numEntityBlocks numNodes minTag maxTag
22 lines.append(f"2 1 0 {n}") # entityDim entityTag parametric numNodes
23 lines.extend(str(i) for i in range(1, n + 1))
24 lines.extend(f"{x:.16g} {y:.16g} {z:.16g}" for x, y, z in points)
28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 25 lines.append("$EndNodes")
26 lines.append("$Elements")
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 27 lines.append(f"1 {m} 1 {m}") # numEntityBlocks numElements minTag maxTag
28 lines.append(f"2 1 3 {m}") # entityDim entityTag elementType(3=quad) numElements
28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 29 for i, q in enumerate(quads, start=1):
30 a, b, c, d = (v + 1 for v in q) # 0-based -> 1-based
43c610bUse surfacemesh.import instead of the hand-rolled mesh readersJeremy Magland 31 lines.append(f"{i} {a} {b} {c} {d}")
28f8ec1mesh-pde-solver: upload a quad mesh, solve PDEs on the surface in-browserJeremy Magland 32 lines.append("$EndElements")
33 with open(path, "w") as f:
34 f.write("\n".join(lines) + "\n")
37def cubed_sphere(m):
38 """Cube [-1,1]^3 with m-by-m quads per face, projected onto the sphere."""
39 points = []
40 index = {}
42 def vertex(x, y, z):
43 # Normalize onto the unit sphere; dedup shared face-boundary vertices.
44 r = math.sqrt(x * x + y * y + z * z)
45 p = (x / r, y / r, z / r)
46 key = tuple(round(v, 12) for v in p)
47 if key not in index:
48 index[key] = len(points)
49 points.append(p)
50 return index[key]
52 # Each face: origin corner + two axis vectors spanning the face.
53 faces = [
54 ((-1, -1, 1), (1, 0, 0), (0, 1, 0)), # +z
55 ((-1, 1, -1), (1, 0, 0), (0, -1, 0)), # -z
56 ((-1, -1, -1), (0, 1, 0), (0, 0, 1)), # -x
57 ((1, -1, -1), (0, 0, 1), (0, 1, 0)), # +x
58 ((-1, -1, -1), (0, 0, 1), (1, 0, 0)), # -y
59 ((-1, 1, -1), (1, 0, 0), (0, 0, 1)), # +y
60 ]
61 quads = []
62 for origin, du, dv in faces:
63 for i in range(m):
64 for j in range(m):
65 corners = []
66 for di, dj in ((0, 0), (1, 0), (1, 1), (0, 1)):
67 s = 2 * (i + di) / m
68 t = 2 * (j + dj) / m
69 corners.append(
70 vertex(
71 origin[0] + s * du[0] + t * dv[0],
72 origin[1] + s * du[1] + t * dv[1],
73 origin[2] + s * du[2] + t * dv[2],
74 )
75 )
76 quads.append(tuple(corners))
77 return points, quads
80def torus(nu, nv, R=1.0, r=0.4):
81 points = []
82 for i in range(nu):
83 a = 2 * math.pi * i / nu
84 for j in range(nv):
85 b = 2 * math.pi * j / nv
86 points.append(
87 (
88 (R + r * math.cos(b)) * math.cos(a),
89 (R + r * math.cos(b)) * math.sin(a),
90 r * math.sin(b),
91 )
92 )
93 quads = []
94 for i in range(nu):
95 for j in range(nv):
96 i2 = (i + 1) % nu
97 j2 = (j + 1) % nv
98 quads.append((i * nv + j, i2 * nv + j, i2 * nv + j2, i * nv + j2))
99 return points, quads
102def main():
103 os.makedirs(OUT_DIR, exist_ok=True)
104 pts, quads = cubed_sphere(6)
105 write_msh(os.path.join(OUT_DIR, "sphere.msh"), pts, quads)
106 print(f"sphere.msh: {len(pts)} nodes, {len(quads)} quads")
107 pts, quads = torus(24, 12)
108 write_msh(os.path.join(OUT_DIR, "torus.msh"), pts, quads)
109 print(f"torus.msh: {len(pts)} nodes, {len(quads)} quads")
112if __name__ == "__main__":
113 main()