1"""Generate the bundled sample meshes (public/samples/*.msh).
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 3-node triangle (type 2) or 4-node quadrangle (type 3) elements — the
6layout surfacemesh.import reads.
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, cells):
18 n, m = len(points), len(cells)
19 etype = 2 if len(cells[0]) == 3 else 3
20 lines = ["$MeshFormat", "4.1 0 8", "$EndMeshFormat"]
21 lines.append("$Nodes")
22 lines.append(f"1 {n} 1 {n}") # numEntityBlocks numNodes minTag maxTag
23 lines.append(f"2 1 0 {n}") # entityDim entityTag parametric numNodes
24 lines.extend(str(i) for i in range(1, n + 1))
25 lines.extend(f"{x:.16g} {y:.16g} {z:.16g}" for x, y, z in points)
26 lines.append("$EndNodes")
27 lines.append("$Elements")
28 lines.append(f"1 {m} 1 {m}") # numEntityBlocks numElements minTag maxTag
29 lines.append(f"2 1 {etype} {m}") # entityDim entityTag elementType numElements
30 for i, cell in enumerate(cells, start=1):
31 lines.append(f"{i} " + " ".join(str(v + 1) for v in cell)) # 0- -> 1-based
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 icosphere(subdiv):
81 """Icosahedron subdivided `subdiv` times, projected onto the unit sphere."""
82 phi = (1 + math.sqrt(5)) / 2
83 norm = math.sqrt(1 + phi * phi)
84 points = [
85 (x / norm, y / norm, z / norm)
86 for x, y, z in (
87 (-1, phi, 0), (1, phi, 0), (-1, -phi, 0), (1, -phi, 0),
88 (0, -1, phi), (0, 1, phi), (0, -1, -phi), (0, 1, -phi),
89 (phi, 0, -1), (phi, 0, 1), (-phi, 0, -1), (-phi, 0, 1),
90 )
91 ]
92 tris = [
93 (0, 11, 5), (0, 5, 1), (0, 1, 7), (0, 7, 10), (0, 10, 11),
94 (1, 5, 9), (5, 11, 4), (11, 10, 2), (10, 7, 6), (7, 1, 8),
95 (3, 9, 4), (3, 4, 2), (3, 2, 6), (3, 6, 8), (3, 8, 9),
96 (4, 9, 5), (2, 4, 11), (6, 2, 10), (8, 6, 7), (9, 8, 1),
97 ]
98 midpoints = {}
100 def midpoint(a, b):
101 key = (a, b) if a < b else (b, a)
102 if key not in midpoints:
103 x = points[a][0] + points[b][0]
104 y = points[a][1] + points[b][1]
105 z = points[a][2] + points[b][2]
106 r = math.sqrt(x * x + y * y + z * z)
107 midpoints[key] = len(points)
108 points.append((x / r, y / r, z / r))
109 return midpoints[key]
111 for _ in range(subdiv):
112 split = []
113 for a, b, c in tris:
114 ab, bc, ca = midpoint(a, b), midpoint(b, c), midpoint(c, a)
115 split += [(a, ab, ca), (ab, b, bc), (ca, bc, c), (ab, bc, ca)]
116 tris = split
117 return points, tris
120def torus(nu, nv, R=1.0, r=0.4):
121 points = []
122 for i in range(nu):
123 a = 2 * math.pi * i / nu
124 for j in range(nv):
125 b = 2 * math.pi * j / nv
126 points.append(
127 (
128 (R + r * math.cos(b)) * math.cos(a),
129 (R + r * math.cos(b)) * math.sin(a),
130 r * math.sin(b),
131 )
132 )
133 quads = []
134 for i in range(nu):
135 for j in range(nv):
136 i2 = (i + 1) % nu
137 j2 = (j + 1) % nv
138 quads.append((i * nv + j, i2 * nv + j, i2 * nv + j2, i * nv + j2))
139 return points, quads
142def main():
143 os.makedirs(OUT_DIR, exist_ok=True)
144 pts, quads = cubed_sphere(6)
145 write_msh(os.path.join(OUT_DIR, "sphere.msh"), pts, quads)
146 print(f"sphere.msh: {len(pts)} nodes, {len(quads)} quads")
147 pts, tris = icosphere(2)
148 write_msh(os.path.join(OUT_DIR, "sphere-tri.msh"), pts, tris)
149 print(f"sphere-tri.msh: {len(pts)} nodes, {len(tris)} triangles")
150 pts, quads = torus(24, 12)
151 write_msh(os.path.join(OUT_DIR, "torus.msh"), pts, quads)
152 print(f"torus.msh: {len(pts)} nodes, {len(quads)} quads")
155if __name__ == "__main__":
156 main()