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).
3Writes Gmsh MSH 2.2 ASCII files in the same canonical form the in-app
4converter produces: sequential 1-based node ids and 4-node quadrangle
5elements (type 3) with two tags.
7Run from the repo root: python3 scripts/make_samples.py
8"""
10import math
11import os
13OUT_DIR = os.path.join(os.path.dirname(__file__), "..", "public", "samples")
16def write_msh(path, points, quads):
17 lines = ["$MeshFormat", "2.2 0 8", "$EndMeshFormat"]
18 lines.append("$Nodes")
19 lines.append(str(len(points)))
20 for i, (x, y, z) in enumerate(points, start=1):
21 lines.append(f"{i} {x:.16g} {y:.16g} {z:.16g}")
22 lines.append("$EndNodes")
23 lines.append("$Elements")
24 lines.append(str(len(quads)))
25 for i, q in enumerate(quads, start=1):
26 a, b, c, d = (v + 1 for v in q) # 0-based -> 1-based
27 lines.append(f"{i} 3 2 1 1 {a} {b} {c} {d}")
28 lines.append("$EndElements")
29 with open(path, "w") as f:
30 f.write("\n".join(lines) + "\n")
33def cubed_sphere(m):
34 """Cube [-1,1]^3 with m-by-m quads per face, projected onto the sphere."""
35 points = []
36 index = {}
38 def vertex(x, y, z):
39 # Normalize onto the unit sphere; dedup shared face-boundary vertices.
40 r = math.sqrt(x * x + y * y + z * z)
41 p = (x / r, y / r, z / r)
42 key = tuple(round(v, 12) for v in p)
43 if key not in index:
44 index[key] = len(points)
45 points.append(p)
46 return index[key]
48 # Each face: origin corner + two axis vectors spanning the face.
49 faces = [
50 ((-1, -1, 1), (1, 0, 0), (0, 1, 0)), # +z
51 ((-1, 1, -1), (1, 0, 0), (0, -1, 0)), # -z
52 ((-1, -1, -1), (0, 1, 0), (0, 0, 1)), # -x
53 ((1, -1, -1), (0, 0, 1), (0, 1, 0)), # +x
54 ((-1, -1, -1), (0, 0, 1), (1, 0, 0)), # -y
55 ((-1, 1, -1), (1, 0, 0), (0, 0, 1)), # +y
56 ]
57 quads = []
58 for origin, du, dv in faces:
59 for i in range(m):
60 for j in range(m):
61 corners = []
62 for di, dj in ((0, 0), (1, 0), (1, 1), (0, 1)):
63 s = 2 * (i + di) / m
64 t = 2 * (j + dj) / m
65 corners.append(
66 vertex(
67 origin[0] + s * du[0] + t * dv[0],
68 origin[1] + s * du[1] + t * dv[1],
69 origin[2] + s * du[2] + t * dv[2],
70 )
71 )
72 quads.append(tuple(corners))
73 return points, quads
76def torus(nu, nv, R=1.0, r=0.4):
77 points = []
78 for i in range(nu):
79 a = 2 * math.pi * i / nu
80 for j in range(nv):
81 b = 2 * math.pi * j / nv
82 points.append(
83 (
84 (R + r * math.cos(b)) * math.cos(a),
85 (R + r * math.cos(b)) * math.sin(a),
86 r * math.sin(b),
87 )
88 )
89 quads = []
90 for i in range(nu):
91 for j in range(nv):
92 i2 = (i + 1) % nu
93 j2 = (j + 1) % nv
94 quads.append((i * nv + j, i2 * nv + j, i2 * nv + j2, i * nv + j2))
95 return points, quads
98def main():
99 os.makedirs(OUT_DIR, exist_ok=True)
100 pts, quads = cubed_sphere(6)
101 write_msh(os.path.join(OUT_DIR, "sphere.msh"), pts, quads)
102 print(f"sphere.msh: {len(pts)} nodes, {len(quads)} quads")
103 pts, quads = torus(24, 12)
104 write_msh(os.path.join(OUT_DIR, "torus.msh"), pts, quads)
105 print(f"torus.msh: {len(pts)} nodes, {len(quads)} quads")
108if __name__ == "__main__":
109 main()