9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 1"""Green's identity S_k[du/dn] - D_k[u] = u/2 at boundary nodes, brute force over all patches,
2with three models of the patch geometry used inside the integrals:
3 interp: xyz, r_u, r_v each interpolated from node values (fmm3dbie's model)
4 consistent: xyz interpolated, r_u and r_v the derivatives of that interpolant
5 exact: the analytic superellipsoid parametrization
6The density (u and du/dn) is always the polynomial interpolant of node values."""
7import os
8import sys
10import numpy as np
12sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
13sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
15from fmm3dbie_gpu import surface as S # noqa: E402
bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 16from problems.common import green, green_dn # noqa: E402
9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 17from test_quadrature import kernels6 # noqa: E402
19a_, p_, m_, norder = 1.0, 8, 4, 6
20surf = S.rounded_cube(a=a_, p=p_, m=m_, norder=norder)
21br = S.graded_breaks(m_, 0.5)
22npols = surf.npols
23k = 10.0; xe = np.array([1.7, 0.6, -0.5])
24u = green(surf.srcvals[0:3], xe, k); un = green_dn(surf.srcvals[0:3], surf.srcvals[9:12], xe, k)
26def exact_geom(patch, uv):
27 f, rem = divmod(patch, m_ * m_); i, j = divmod(rem, m_)
28 n, e1, e2 = (np.array(v, float) for v in S._FACES[f])
29 s0, s1, t0, t1 = br[i], br[i + 1], br[j], br[j + 1]
30 s = s0 + 0.5 * (uv[0] + 1) * (s1 - s0); t = t0 + 0.5 * (uv[1] + 1) * (t1 - t0)
31 c = n[:, None] + s * e1[:, None] + t * e2[:, None]
32 N = np.sum(c ** p_, axis=0) ** (1.0 / p_); gradN = c ** (p_ - 1) / N ** (p_ - 1)
33 def deriv(dc, scale):
34 dN = np.sum(gradN * dc[:, None], axis=0)
35 return a_ * scale * (dc[:, None] / N - c * dN / N ** 2)
36 return np.vstack([a_ * c / N, deriv(e1, 0.5 * (s1 - s0)), deriv(e2, 0.5 * (t1 - t0))])
38def dlegendre_vander(x, n):
39 P = S.legendre_vander(x, n); D = np.zeros_like(P)
40 for k_ in range(1, n + 1): # P_k' = sum over j = k-1, k-3, ... of (2j+1) P_j
41 for j in range(k_ - 1, -1, -2):
42 D[k_] += (2 * j + 1) * P[j]
43 return D
45def consistent_geom(patch, uv):
46 c = surf.srccoefs[0:3, patch * npols:(patch + 1) * npols]
47 Pu, Pv = S.legendre_vander(uv[0], norder), S.legendre_vander(uv[1], norder)
48 Du, Dv = dlegendre_vander(uv[0], norder), dlegendre_vander(uv[1], norder)
49 B = (Pv[:, None] * Pu[None]).reshape(npols, -1)
50 Bu = (Pv[:, None] * Du[None]).reshape(npols, -1)
51 Bv = (Dv[:, None] * Pu[None]).reshape(npols, -1)
52 return np.vstack([c @ B, c @ Bu, c @ Bv])
54models = {"interp": lambda P, uv: surf.eval_geometry(P, uv), "consistent": consistent_geom, "exact": exact_geom}
56def moments(geom, patch, uv, w, xt, nt):
57 g = geom(patch, uv)
58 ns = np.cross(g[3:6].T, g[6:9].T).T; jac = np.linalg.norm(ns, axis=0); ns = ns / jac
59 K = kernels6(xt, nt, g[0:3], ns, k)
60 return (K * (jac * w)[None]) @ S.basis_2d(uv, norder).T
62def self_rule(u0, v0, nx=200, nr=40):
63 xr, wr = np.polynomial.legendre.leggauss(nr); xr, wr = 0.5 * (xr + 1), 0.5 * wr
64 xx, wx = np.polynomial.legendre.leggauss(nx)
65 U, W = [], []
66 for e in range(4):
67 h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e]
68 x1, x2 = (-1 - v0, 1 - v0) if e < 2 else (-1 - u0, 1 - u0)
69 for xa, xb in [(x1, 0.0), (0.0, x2)]:
70 x = xa + 0.5 * (xb - xa) * (xx + 1); dx = 0.5 * (xb - xa) * wx
71 R, X = np.meshgrid(xr, x, indexing="ij"); WW = np.outer(wr, dx) * abs(h) * R
72 uv = np.vstack([(u0 + R * h).ravel(), (v0 + R * X).ravel()]) if e < 2 else \
73 np.vstack([(u0 + R * X).ravel(), (v0 + R * h).ravel()])
74 U.append(uv); W.append(WW.ravel())
75 return np.hstack(U), np.hstack(W)
77def leaf_rule(depth, nq=8):
78 gx, gw = np.polynomial.legendre.leggauss(nq); nl = 2 ** depth; h = 1.0 / nl
79 c = -1 + h + 2 * h * np.arange(nl)
80 Uc, Vc, A, B = np.meshgrid(c, c, gx, gx, indexing="ij")
81 uv = np.vstack([(Uc + h * A).ravel(), (Vc + h * B).ravel()])
82 w = (h * h * np.outer(gw, gw)[None, None]).repeat(nl, 0).repeat(nl, 1).ravel()
83 return uv, w
85leaf6, leaf3 = leaf_rule(6), leaf_rule(3)
86P = 13
87for loc in [0, 3, 21, 24]:
88 i = P * npols + loc
89 xt, nt = surf.srcvals[0:3, i], surf.srcvals[9:12, i]
90 out = []
91 for name, geom in models.items():
92 Sv = Dv = 0
93 for j in range(surf.npatches):
94 if j == P:
95 uv, w = self_rule(*surf.uv_nodes[:, i])
96 else:
97 d = np.linalg.norm(xt - surf.centroids[:, j]) / surf.rads[j]
98 uv, w = leaf6 if d < 3 else leaf3
99 Wt = moments(geom, j, uv, w, xt, nt) @ surf.umat
100 sl = slice(j * npols, (j + 1) * npols)
101 Sv += Wt[4] @ un[sl]; Dv += Wt[5] @ u[sl]
102 out.append(f"{name} {abs(Sv - Dv - u[i] / 2) / abs(u[i]):.1e}")
103 print(f"node position {loc} (u,v)=({surf.uv_nodes[0, i]:+.2f},{surf.uv_nodes[1, i]:+.2f}): relative residual " + " ".join(out), flush=True)