1"""Triangular mesh of the projected disk and the per-element geometry derived from it.
3The mesher needs only numpy. Nodes are placed on concentric rings with spacing equal to the
4target element size; ring k carries round(2 pi r_k / dr) nodes, which is close to 6k, so the
5triangles are nearly equilateral. Alternate rings are rotated by half a step to avoid radially
6aligned edges. Consecutive rings are stitched by an angular merge: walking both rings with
7increasing angle, the next triangle uses whichever ring's next node comes first, which yields a
8valid triangulation of the annulus with (n_inner + n_outer) triangles. The centre is a single node
9fanned to the first ring. See docs/model.md, section "Spatial discretization".
10"""
11from __future__ import annotations
13from dataclasses import dataclass, field
15import numpy as np
19class Mesh:
20 """Nodes ``xy`` (N, 2), counter-clockwise triangles ``tri`` (E, 3) and the rim node indices.
22 Derived quantities, computed once: element areas, the constant gradient operators Gx, Gy of
23 the linear (CST) interpolation, and the centroids.
24 """
25 xy: np.ndarray
26 tri: np.ndarray
27 boundary: np.ndarray
28 area: np.ndarray = field(init=False, repr=False)
29 Gx: np.ndarray = field(init=False, repr=False)
30 Gy: np.ndarray = field(init=False, repr=False)
31 centroid: np.ndarray = field(init=False, repr=False)
33 def __post_init__(self):
34 self.xy = np.asarray(self.xy, dtype=np.float64)
35 self.tri = np.asarray(self.tri, dtype=np.int64)
36 x, y = self.xy[self.tri, 0], self.xy[self.tri, 1]
37 x21, x31 = x[:, 1] - x[:, 0], x[:, 2] - x[:, 0]
38 y21, y31 = y[:, 1] - y[:, 0], y[:, 2] - y[:, 0]
39 self.area = 0.5 * (x21 * y31 - x31 * y21)
40 if not (self.area > 0).all():
41 raise ValueError("mesh has non-positive element areas (orientation must be counter-clockwise)")
42 # d(phi)/dx = sum_i Gx[i] phi_i for the linear interpolant on the triangle
43 self.Gx = np.stack([y[:, 1] - y[:, 2], y[:, 2] - y[:, 0], y[:, 0] - y[:, 1]], 1) / (2 * self.area[:, None])
44 self.Gy = np.stack([x[:, 2] - x[:, 1], x[:, 0] - x[:, 2], x[:, 1] - x[:, 0]], 1) / (2 * self.area[:, None])
45 self.centroid = np.stack([x.mean(1), y.mean(1)], 1)
48 def n_nodes(self) -> int:
49 return len(self.xy)
52 def n_elements(self) -> int:
53 return len(self.tri)
55 def edge_lengths(self) -> np.ndarray:
56 """All three edge lengths of every element, shape (E, 3)."""
57 p = self.xy[self.tri]
58 return np.linalg.norm(p - np.roll(p, -1, axis=1), axis=2)
60 def locate(self, px: float, py: float):
61 """Containing element's nodes (3,) and barycentric weights (3,) of a point of the domain."""
62 x, y = self.xy[self.tri, 0], self.xy[self.tri, 1]
63 a2 = 2 * self.area
64 l1 = ((y[:, 1] - y[:, 2]) * (px - x[:, 2]) + (x[:, 2] - x[:, 1]) * (py - y[:, 2])) / a2
65 l2 = ((y[:, 2] - y[:, 0]) * (px - x[:, 2]) + (x[:, 0] - x[:, 2]) * (py - y[:, 2])) / a2
66 l3 = 1 - l1 - l2
67 inside = (l1 >= -1e-9) & (l2 >= -1e-9) & (l3 >= -1e-9)
68 if not inside.any():
69 raise ValueError(f"point ({px:.4f}, {py:.4f}) lies outside the mesh")
70 e = int(np.argmax(inside))
71 return self.tri[e].copy(), np.array([l1[e], l2[e], l3[e]])
73 def interpolate(self, values: np.ndarray, px: float, py: float) -> float:
74 nodes, w = self.locate(px, py)
75 return float(values[nodes] @ w)
78def _stitch(a_idx, a_ang, b_idx, b_ang):
79 """Triangulate the annulus between an inner ring (a) and an outer ring (b).
81 Both rings are given as node indices and angles. Angles are unwrapped so that every outer
82 angle lies in [a_ang[0], a_ang[0] + 2 pi); then the two rings are merged in increasing angle.
83 """
84 na, nb = len(a_idx), len(b_idx)
85 b_ang = (b_ang - a_ang[0]) % (2 * np.pi) + a_ang[0]
86 order = np.argsort(b_ang)
87 b_idx, b_ang = b_idx[order], b_ang[order]
88 ta = np.append(a_ang, a_ang[0] + 2 * np.pi)
89 tb = np.append(b_ang, b_ang[0] + 2 * np.pi)
90 out = []
91 i = j = 0
92 while i < na or j < nb:
93 next_a = ta[i + 1] if i < na else np.inf
94 next_b = tb[j + 1] if j < nb else np.inf
95 if next_a <= next_b:
96 out.append((a_idx[i % na], b_idx[j % nb], a_idx[(i + 1) % na]))
97 i += 1
98 else:
99 out.append((a_idx[i % na], b_idx[j % nb], b_idx[(j + 1) % nb]))
100 j += 1
101 return out
104def _smooth(xy: np.ndarray, tri: np.ndarray, fixed: np.ndarray, n_iter: int) -> np.ndarray:
105 """Laplacian smoothing: move each free node to the mean of its neighbours."""
106 N = len(xy)
107 edges = np.concatenate([tri[:, [0, 1]], tri[:, [1, 2]], tri[:, [2, 0]]])
108 src = np.concatenate([edges[:, 0], edges[:, 1]])
109 dst = np.concatenate([edges[:, 1], edges[:, 0]])
110 deg = np.bincount(src, minlength=N)
111 free = np.ones(N, bool)
112 free[fixed] = False
113 xy = xy.copy()
114 for _ in range(n_iter):
115 acc = np.zeros_like(xy)
116 np.add.at(acc, src, xy[dst])
117 mean = acc / deg[:, None]
118 xy[free] = mean[free]
119 return xy
122def disk_mesh(radius: float, element_size: float, smoothing: int = 1) -> Mesh:
123 """Near-equilateral triangulation of the disk of the given radius. The rim nodes stay on
124 the circle; interior nodes are relaxed by ``smoothing`` Laplacian passes. One pass removes
125 the worst stitching triangles; more passes improve the mean shape but shrink the smallest
126 elements, which lowers the stable time step (see docs/model.md, section 4.1)."""
127 n_rings = max(2, int(round(radius / element_size)))
128 dr = radius / n_rings
129 nodes = [(0.0, 0.0)]
130 prev_idx, prev_ang = np.array([0]), np.array([0.0])
131 tris: list[tuple[int, int, int]] = []
132 for k in range(1, n_rings + 1):
133 r = k * dr
134 n_k = max(6, int(round(2 * np.pi * r / dr)))
135 ang = 2 * np.pi * (np.arange(n_k) + 0.5 * (k % 2)) / n_k
136 idx = len(nodes) + np.arange(n_k)
137 nodes.extend(zip(r * np.cos(ang), r * np.sin(ang)))
138 if k == 1:
139 tris.extend((0, idx[j], idx[(j + 1) % n_k]) for j in range(n_k))
140 else:
141 tris.extend(_stitch(prev_idx, prev_ang, idx, ang))
142 prev_idx, prev_ang = idx, ang
143 xy = np.array(nodes, dtype=np.float64)
144 tri = np.array(tris, dtype=np.int64)
145 if smoothing > 0:
146 xy = _smooth(xy, tri, prev_idx, smoothing)
147 x, y = xy[tri, 0], xy[tri, 1]
148 signed = (x[:, 1] - x[:, 0]) * (y[:, 2] - y[:, 0]) - (x[:, 2] - x[:, 0]) * (y[:, 1] - y[:, 0])
149 tri[signed < 0] = tri[signed < 0][:, [0, 2, 1]]
150 return Mesh(xy=xy, tri=tri, boundary=prev_idx)