1"""A cymbal as a thin closed body built from quadrilateral patches (pure NumPy).
3fmm3dbie's integral equations need a closed surface, so the cymbal is not a
4zero-thickness shell but a thin body around an axisymmetric midsurface
5z = F(rho), rho = x^2 + y^2, with thickness T(rho). The midsurface has a raised
6bell at the center and a shallow bow outside it; the thickness decreases from
7the center to the rim, with zero slope at the rim. The top and bottom faces
8are the offset surfaces M +- (T/2) N, with M the midsurface and N its upward
9unit normal, and a half-torus of diameter T closes the body at the rim. Every
10patch lies entirely on one face or on the rim, so every patch is analytic even
11though the curvature jumps where the faces meet the rim.
13The faces are meshed with an O-grid in the (x, y) plane: an equiangular central
14square of m_c x m_c patches, then rings of patches between the square and a
15circle and between successive circles. Rings are graded toward the rim, where
16the density varies on the scale of the thickness. Neighboring rings need not
17have the same number of patches around, so the mesh is not conforming; neither
18fmm3dbie nor the GPU code requires that it be.
20Node derivatives dx/du, dx/dv are computed by the complex step method, which is
21exact to rounding for this analytic parametrization.
22"""
24import numpy as np
26from . import surface
28CSTEP = 1e-30
31class CymbalShape:
32 """Axisymmetric cymbal shape; lengths in meters.
34 radius: outer radius, including the rim.
35 t_rim: thickness at the rim; taper: thickness at the center / t_rim.
36 bell_radius, bell_height, bell_width: the bell is a logistic step in rho,
37 of height bell_height, centered at r = bell_radius, with width about
38 bell_width in r.
39 bow_height: fall of the midsurface from the center to the rim, apart from the bell.
40 """
42 def __init__(self, radius=0.2, t_rim=0.01, taper=2.0, bell_radius=0.06, bell_height=0.025,
43 bell_width=0.01, bow_height=0.015):
44 self.radius, self.t_rim, self.taper = radius, t_rim, taper
45 self.bell_radius, self.bell_height, self.bell_width = bell_radius, bell_height, bell_width
46 self.bow_height = bow_height
47 self.r_mid = radius - t_rim / 2 # radius of the midsurface, where the rim is attached
49 # midsurface height, its rho-derivative, and the thickness, all functions of rho
50 def _bell_arg(self, rho):
51 return (rho - self.bell_radius ** 2) / (2 * self.bell_radius * self.bell_width)
53 def F(self, rho):
54 return (self.bow_height * (1 - rho / self.r_mid ** 2)
55 + self.bell_height / (1 + np.exp(self._bell_arg(rho))))
57 def dF(self, rho):
58 e = np.exp(self._bell_arg(rho))
59 return (-self.bow_height / self.r_mid ** 2
60 - self.bell_height * e / (1 + e) ** 2 / (2 * self.bell_radius * self.bell_width))
62 def T(self, rho):
63 return self.t_rim * (1 + (self.taper - 1) * (1 - rho / self.r_mid ** 2) ** 2)
65 def face_point(self, x, y, side):
66 """Point of the top (side = +1) or bottom (side = -1) face above/below (x, y) of
67 the midsurface. Complex-safe."""
68 rho = x * x + y * y
69 g = 2 * self.dF(rho) # grad F = g (x, y)
70 q = np.sqrt(1 + g * g * rho)
71 s = side * self.T(rho) / 2 / q
72 return np.array([x - s * g * x, y - s * g * y, self.F(rho) + s])
74 def rim_point(self, theta, phi):
75 """Point of the rim at azimuth theta; phi from -pi/2 (bottom face) to pi/2 (top face)."""
76 R = self.r_mid
77 rho = R * R
78 fp = 2 * R * self.dF(rho) # dz/dr of the midsurface at the rim
79 q = np.sqrt(1 + fp * fp)
80 c, s = np.cos(theta), np.sin(theta)
81 h = self.T(rho) / 2
82 a, b = h * np.cos(phi) / q, h * np.sin(phi) / q # along the tangent (1, fp) and normal (-fp, 1)
83 rr = R + a - b * fp
84 return np.array([rr * c, rr * s, self.F(rho) + a * fp + b])
86 def cross_section(self, n=400):
87 """Closed outline of the body in the (r, z) half plane, (2, m), counterclockwise."""
88 r = np.linspace(0, self.r_mid, n)
89 top = self.face_point(r, 0 * r, 1)
90 bot = self.face_point(r, 0 * r, -1)
91 phi = np.linspace(-np.pi / 2, np.pi / 2, n // 4)
92 rim = self.rim_point(0 * phi, phi)
93 rz = np.hstack([bot[[0, 2]], rim[[0, 2]], top[[0, 2]][:, ::-1]])
94 return rz
96 def inside(self, pts, margin=0.0):
97 """True for points (3, n) inside the body or within margin of its surface.
99 Every vertical line meets the cross section in one interval, so we offset the
100 outline outward by margin, split it at its outermost point into a lower and an
101 upper chain, and compare z with both chains at the point's radius."""
102 rz = self.cross_section(4000)
103 if margin > 0:
104 t = np.roll(rz, -1, axis=1) - np.roll(rz, 1, axis=1) # counterclockwise tangent
105 t /= np.linalg.norm(t, axis=0)
106 rz = rz + margin * np.array([t[1], -t[0]]) # outward normal
107 rz[0] = np.maximum(rz[0], 0)
108 i = int(np.argmax(rz[0]))
109 lo, up = rz[:, :i + 1], rz[:, i:][:, ::-1]
110 r = np.hypot(pts[0], pts[1])
111 z_lo = np.interp(r, np.maximum.accumulate(lo[0]), lo[1])
112 z_up = np.interp(r, np.maximum.accumulate(up[0]), up[1])
113 return (r <= rz[0, i]) & (pts[2] > z_lo) & (pts[2] < z_up)
115 def volume(self, n=20000):
116 """Volume of the body by Pappus' theorem on the cross section (a check)."""
117 rz = self.cross_section(n)
118 r, z = rz
119 r1, z1 = np.roll(r, -1), np.roll(z, -1)
120 # int r dA over the cross section, by Green's theorem: (1/2) oint r^2 dz
121 return 2 * np.pi * np.sum((r * r + r * r1 + r1 * r1) / 6 * (z1 - z))
124def ring_layout(shape, h, grade=1.5, c_frac=0.12):
125 """Radii and patch counts of the O-grid for a target patch size h.
127 Returns (m_c, c, radii, counts): the central square has half-width c and
128 m_c x m_c patches (m_c >= 2); ring i lies between radii[i] and radii[i + 1]
129 (radii[0] = c means the square) and has counts[i] patches around (a multiple
130 of 4). The ring around the square reaches r = 2c, so that it is 0.59c wide at
131 the square's corners. Radial widths are at most h; they grow geometrically
132 from t_rim at the rim by the factor grade, and they are at most twice the
133 bell width across the bell.
134 """
135 R = shape.r_mid
136 c = c_frac * R
137 m_c = max(2, int(np.ceil(2 * c / h)))
138 first = 2 * c # outer radius of the ring around the square
139 b0 = shape.bell_radius - 2.5 * shape.bell_width
140 b1 = shape.bell_radius + 2.5 * shape.bell_width
142 def wmax(edge):
143 w = min(h, shape.t_rim + (grade - 1) * (R - edge))
144 if edge > b0 and edge - w < b1:
145 w = min(w, 2 * shape.bell_width)
146 return w
148 radii = [R]
149 while True:
150 w = wmax(radii[-1])
151 if radii[-1] - w < first + 0.5 * w:
152 break
153 radii.append(radii[-1] - w)
154 radii = np.array([c, first] + radii[::-1])
155 counts = []
156 for i in range(len(radii) - 1):
157 width = radii[i + 1] - radii[i]
158 arc = min(h, 2.5 * width)
159 n = 4 * max(m_c, int(np.ceil(2 * np.pi * radii[i + 1] / (4 * arc))))
160 counts.append(n)
161 counts[0] = 4 * m_c # the first ring matches the square
162 # rings never have fewer patches around than the ring inside them
163 counts = list(np.maximum.accumulate(counts))
164 return m_c, c, radii, counts
167def _face_charts(m_c, c, radii, counts):
168 """Planar charts (u, v) in [-1, 1]^2 -> (x, y) for all face patches, as closures.
169 Each chart is oriented so that d/du x d/dv points up."""
170 charts = []
171 q = np.pi / 4
172 # central square, equiangular so that its edges match the first ring
173 br = np.linspace(-1, 1, m_c + 1)
174 for i in range(m_c):
175 for j in range(m_c):
176 def ch(u, v, a0=br[i], a1=br[i + 1], b0=br[j], b1=br[j + 1]):
177 s = a0 + 0.5 * (u + 1) * (a1 - a0)
178 t = b0 + 0.5 * (v + 1) * (b1 - b0)
179 return c * np.tan(q * s), c * np.tan(q * t)
180 charts.append(ch)
181 # rings; u is radial, v is the angle (counterclockwise)
182 for ir in range(len(radii) - 1):
183 r0, r1, n = radii[ir], radii[ir + 1], counts[ir]
184 dth = 2 * np.pi / n
185 for k in range(n):
186 th0 = -q + k * dth # patch angles start at the square corner (-pi/4)
187 if ir == 0:
188 def ch(u, v, th0=th0, r1=r1, dth=dth):
189 th = th0 + 0.5 * (v + 1) * dth
190 sec = np.floor((th.real + q) / (2 * q)) # sector index 0..3
191 loc = th - sec * 2 * q # angle relative to the sector's center
192 rin = c / np.cos(loc)
193 a = 0.5 * (u + 1)
194 rad = (1 - a) * rin + a * r1
195 return rad * np.cos(th), rad * np.sin(th)
196 else:
197 def ch(u, v, th0=th0, r0=r0, r1=r1, dth=dth):
198 th = th0 + 0.5 * (v + 1) * dth
199 rad = r0 + 0.5 * (u + 1) * (r1 - r0)
200 return rad * np.cos(th), rad * np.sin(th)
201 charts.append(ch)
202 return charts
205def _node_data(fun, uvs):
206 """xyz, r_u, r_v, unit normal (12, n) of the map fun(u, v) -> (3, n), with complex step derivatives."""
207 u, v = uvs[0].astype(complex), uvs[1].astype(complex)
208 x = np.real(fun(u, v))
209 du = np.imag(fun(u + 1j * CSTEP, v)) / CSTEP
210 dv = np.imag(fun(u, v + 1j * CSTEP)) / CSTEP
211 nrm = np.cross(du.T, dv.T).T
212 nrm /= np.linalg.norm(nrm, axis=0)
213 return np.vstack([x, du, dv, nrm])
216def patch_maps(shape, h, grade=1.5):
217 """All patch maps (u, v) -> xyz (3, n) with outward orientation, and the ring layout."""
218 m_c, c, radii, counts = ring_layout(shape, h, grade)
219 charts = _face_charts(m_c, c, radii, counts)
220 maps = []
221 for ch in charts: # top faces
222 maps.append(lambda u, v, ch=ch: shape.face_point(*ch(u, v), 1))
223 for ch in charts: # bottom faces: swapping u and v flips the normal
224 maps.append(lambda u, v, ch=ch: shape.face_point(*ch(v, u), -1))
225 nrim = counts[-1]
226 for k in range(nrim): # rim: u is the azimuth, v the angle over the rim
227 th0 = -np.pi / 4 + 2 * np.pi * k / nrim
228 maps.append(lambda u, v, th0=th0: shape.rim_point(th0 + 0.5 * (u + 1) * 2 * np.pi / nrim,
229 0.5 * np.pi * v))
230 return maps, dict(m_c=m_c, c=c, radii=radii, counts=counts, nrim=nrim)
233def cymbal_srcvals(shape, h, norder, grade=1.5):
234 """Node data (12, npts) in the fmm3dbie iptype 11 convention, and the layout."""
235 uvs, _ = surface.disc_nodes(norder)
236 maps, layout = patch_maps(shape, h, grade)
237 return np.hstack([_node_data(f, uvs) for f in maps]), layout
240def cymbal(shape, h, norder, grade=1.5):
241 """A surface.Surface for the cymbal, with the interpolated geometry inside the
242 integrals, as in fmm3dbie."""
243 srcvals, layout = cymbal_srcvals(shape, h, norder, grade)
244 s = surface.Surface(norder, srcvals)
245 s.layout = layout
246 return s