1"""Quadrilateral patch surfaces in the fmm3dbie convention, without the library.
3Each patch is the image of the reference square [-1, 1]^2. It carries
4(norder+1)^2 tensor Gauss-Legendre nodes, with the u index varying fastest,
5and an unnormalized Legendre basis P_j(u) P_i(v) with the same ordering.
6Node data (srcvals) is xyz, dx/du, dx/dv, unit normal, and the geometry is
7also stored as basis coefficients (srccoefs) so that it can be evaluated at
8any point of the reference square. These are the same conventions as
9fmm3dbie's iptype 11, and were checked against the library.
10"""
12import numpy as np
15def gauss_legendre(n):
16 x, w = np.polynomial.legendre.leggauss(n)
17 return x, w
20def legendre_vander(x, norder):
21 """P_0(x), ..., P_norder(x) as an array of shape (norder+1, len(x))."""
22 x = np.asarray(x, dtype=float)
23 P = np.zeros((norder + 1, x.size))
24 P[0] = 1.0
25 if norder >= 1:
26 P[1] = x
27 for n in range(1, norder):
28 P[n + 1] = ((2 * n + 1) * x * P[n] - n * P[n - 1]) / (n + 1)
29 return P
32def disc_nodes(norder):
33 """Reference nodes (2, npols) and weights (npols) of the tensor GL rule."""
34 x, w = gauss_legendre(norder + 1)
35 v, u = np.meshgrid(x, x, indexing="ij") # u fastest
36 wv, wu = np.meshgrid(w, w, indexing="ij")
37 return np.vstack([u.ravel(), v.ravel()]), (wu * wv).ravel()
40def basis_2d(uv, norder):
41 """Basis values (npols, npts): index n = i*(norder+1) + j is P_j(u) P_i(v)."""
42 Pu = legendre_vander(uv[0], norder)
43 Pv = legendre_vander(uv[1], norder)
44 return (Pv[:, None, :] * Pu[None, :, :]).reshape((norder + 1) ** 2, -1)
47def vals_to_coefs_matrix(norder):
48 """U with coefs = U @ vals for values at the tensor GL nodes."""
49 x, w = gauss_legendre(norder + 1)
50 P = legendre_vander(x, norder)
51 U1 = (2 * np.arange(norder + 1)[:, None] + 1) / 2 * P * w[None, :]
52 return np.kron(U1, U1) # (i, j) x (ii, jj) with i, ii slow
55def oversample_matrix(norder, nover):
56 """M with values at the order-nover nodes = M @ values at the order-norder nodes."""
57 uv_o, _ = disc_nodes(nover)
58 return basis_2d(uv_o, norder).T @ vals_to_coefs_matrix(norder)
61class Surface:
62 """exact: optional dict(params=(npatches, 13), a=, p=) giving the superellipsoid face
63 chart of every patch (n, e1, e2, s0, s1, t0, t1); when present the integrals use
64 the exact geometry instead of the polynomial interpolant."""
66 def __init__(self, norder, srcvals, exact=None):
67 self.exact = exact
68 self.norder = norder
69 self.iptype = 11
70 self.srcvals = np.ascontiguousarray(srcvals, dtype=float) # (12, npts)
71 self.npts = self.srcvals.shape[1]
72 self.npols = (norder + 1) ** 2
73 assert self.npts % self.npols == 0
74 self.npatches = self.npts // self.npols
75 self.uvs, self.wref = disc_nodes(norder)
76 self.umat = vals_to_coefs_matrix(norder)
77 vals = self.srcvals[:9].reshape(9, self.npatches, self.npols)
78 self.srccoefs = np.einsum("nj,cpj->cpn", self.umat, vals).reshape(9, self.npts)
79 ru, rv = self.srcvals[3:6], self.srcvals[6:9]
80 self.jac = np.linalg.norm(np.cross(ru.T, rv.T).T, axis=0)
81 self.wts = self.jac * np.tile(self.wref, self.npatches)
82 # centroid and radius of each patch from a 5x5 sample grid (fmm3dbie uses the
83 # four corners, which underestimates the extent of strongly curved patches)
84 g5 = np.linspace(-1, 1, 5)
85 corners = np.array(np.meshgrid(g5, g5, indexing="ij")).reshape(2, -1)
86 Pc = basis_2d(corners, norder) # (npols, 25)
87 xc = np.einsum("cpn,nk->cpk", self.srccoefs[:3].reshape(3, self.npatches, self.npols), Pc)
88 self.centroids = xc.mean(axis=2) # (3, npatches)
89 self.rads = np.linalg.norm(xc - self.centroids[:, :, None], axis=0).max(axis=1)
90 self.patch_id = np.repeat(np.arange(self.npatches), self.npols)
91 self.uv_nodes = np.tile(self.uvs, self.npatches) # (2, npts)
93 def eval_geometry(self, patch, uv):
94 """xyz, r_u, r_v at reference points uv (2, n) of one patch: (9, n)."""
95 if self.exact is not None:
96 e = self.exact
97 return superellipsoid_chart(e["params"][patch], e["a"], e["p"], uv)
98 P = basis_2d(uv, self.norder)
99 c = self.srccoefs[:, patch * self.npols:(patch + 1) * self.npols]
100 return c @ P
102 def oversampled(self, nover):
103 """Nodes of the order-nover rule on every patch: xyz, normal, weights."""
104 if nover == self.norder:
105 return self.srcvals[0:3], self.srcvals[9:12], self.wts
106 uv_o, w_o = disc_nodes(nover)
107 if self.exact is not None:
108 g = np.hstack([self.eval_geometry(j, uv_o) for j in range(self.npatches)])
109 else:
110 P = basis_2d(uv_o, self.norder) # (npols, npols_o)
111 g = np.einsum("cpn,nk->cpk", self.srccoefs.reshape(9, self.npatches, self.npols), P)
112 g = g.reshape(9, -1)
113 nrm = np.cross(g[3:6].T, g[6:9].T).T
114 jac = np.linalg.norm(nrm, axis=0)
115 nrm /= jac
116 return g[0:3], nrm, jac * np.tile(w_o, self.npatches)
119# ---------------------------------------------------------------------------
120# Rounded cube (superellipsoid), the same construction as ../rounded_cube.py
122_FACES = [
123 ((1, 0, 0), (0, 1, 0), (0, 0, 1)),
124 ((-1, 0, 0), (0, 0, 1), (0, 1, 0)),
125 ((0, 1, 0), (0, 0, 1), (1, 0, 0)),
126 ((0, -1, 0), (1, 0, 0), (0, 0, 1)),
127 ((0, 0, 1), (1, 0, 0), (0, 1, 0)),
128 ((0, 0, -1), (0, 1, 0), (1, 0, 0)),
129]
132def graded_breaks(m, grading):
133 t = np.linspace(-1, 1, m + 1)
134 return (1 - grading) * t + grading * np.sin(0.5 * np.pi * t)
137def superellipsoid_chart(par, a, p, uv):
138 """xyz, r_u, r_v (9, n) of the chart par = (n, e1, e2, s0, s1, t0, t1) at uv (2, n)."""
139 n, e1, e2 = par[0:3], par[3:6], par[6:9]
140 s0, s1, t0, t1 = par[9:13]
141 s = s0 + 0.5 * (uv[0] + 1) * (s1 - s0)
142 t = t0 + 0.5 * (uv[1] + 1) * (t1 - t0)
143 c = n[:, None] + s * e1[:, None] + t * e2[:, None]
144 N = np.sum(c**p, axis=0) ** (1.0 / p)
145 gradN = (c / N) ** (p - 1)
147 def deriv(dc, scale):
148 dN = np.sum(gradN * dc[:, None], axis=0)
149 return a * scale * (dc[:, None] / N - c * dN / N**2)
151 return np.vstack([a * c / N, deriv(e1, 0.5 * (s1 - s0)), deriv(e2, 0.5 * (t1 - t0))])
154def rounded_cube(a=1.0, p=8, m=6, norder=8, grading=0.5, exact=True):
155 """Superellipsoid |x|^p+|y|^p+|z|^p = a^p, m x m patches per face.
157 With exact=True the integrals use the analytic parametrization; otherwise the
158 polynomial interpolant of the node data, as fmm3dbie does."""
159 uvs, _ = disc_nodes(norder)
160 br = graded_breaks(m, grading)
161 blocks, params = [], []
162 for n, e1, e2 in _FACES:
163 for i in range(m):
164 for j in range(m):
165 par = np.array(list(n) + list(e1) + list(e2) + [br[i], br[i + 1], br[j], br[j + 1]],
166 dtype=float)
167 params.append(par)
168 g = superellipsoid_chart(par, a, p, uvs)
169 nrm = np.cross(g[3:6].T, g[6:9].T).T
170 nrm /= np.linalg.norm(nrm, axis=0)
171 blocks.append(np.vstack([g, nrm]))
172 ex = dict(params=np.array(params), a=a, p=p) if exact else None
173 return Surface(norder, np.hstack(blocks), exact=ex)
176def superellipsoid_volume(a, p):
177 from math import gamma
178 return 8 * a**3 * gamma(1 + 1 / p) ** 3 / gamma(1 + 3 / p)
181def inside(x, a, p, margin=0.0):
182 return np.sum(np.abs(x / a) ** p, axis=0) < (1 - margin) ** p