1"""Rounded cube (superellipsoid) surface built from quadrilateral patches.
3The surface is |x|^p + |y|^p + |z|^p = a^p with p even. Each of the six cube
4faces (s, t) in [-1, 1]^2 is mapped radially onto the surface,
5x = a c / ||c||_p with c = n + s e1 + t e2. Larger p gives sharper edges.
6Patch breakpoints in s and t are graded toward the face edges, where the
7curvature is concentrated.
8"""
10import numpy as np
12import fmm3dbie_ctypes as fb
14# (n, e1, e2) with e1 x e2 = n, so that d/ds x d/dt points outward
15_FACES = [
16 ((1, 0, 0), (0, 1, 0), (0, 0, 1)),
17 ((-1, 0, 0), (0, 0, 1), (0, 1, 0)),
18 ((0, 1, 0), (0, 0, 1), (1, 0, 0)),
19 ((0, -1, 0), (1, 0, 0), (0, 0, 1)),
20 ((0, 0, 1), (1, 0, 0), (0, 1, 0)),
21 ((0, 0, -1), (0, 1, 0), (1, 0, 0)),
22]
25def graded_breaks(m, grading):
26 """m+1 breakpoints on [-1, 1]; grading in [0, 1] blends uniform and sine-clustered."""
27 t = np.linspace(-1, 1, m + 1)
28 return (1 - grading) * t + grading * np.sin(0.5 * np.pi * t)
31def rounded_cube(a=1.0, p=8, m=6, norder=8, grading=0.5):
32 """Return an fmm3dbie Surface for the superellipsoid of half-width a.
34 m x m patches per face, each discretized with (norder+1)^2 Gauss-Legendre nodes.
35 """
36 iptype = 11
37 uvs = fb.get_disc_nodes(norder, iptype) # (2, npols) on [-1, 1]^2
38 br = graded_breaks(m, grading)
39 blocks = []
40 for n, e1, e2 in _FACES:
41 n, e1, e2 = (np.array(v, dtype=float) for v in (n, e1, e2))
42 for i in range(m):
43 for j in range(m):
44 s0, s1 = br[i], br[i + 1]
45 t0, t1 = br[j], br[j + 1]
46 s = s0 + 0.5 * (uvs[0] + 1) * (s1 - s0)
47 t = t0 + 0.5 * (uvs[1] + 1) * (t1 - t0)
48 c = n[:, None] + s * e1[:, None] + t * e2[:, None] # (3, npols)
49 N = np.sum(c**p, axis=0) ** (1.0 / p)
50 gradN = c ** (p - 1) / N ** (p - 1) # dN/dc
51 xyz = a * c / N
53 def deriv(dc, scale):
54 dN = np.sum(gradN * dc[:, None], axis=0)
55 return a * scale * (dc[:, None] / N - c * dN / N**2)
57 du = deriv(e1, 0.5 * (s1 - s0))
58 dv = deriv(e2, 0.5 * (t1 - t0))
59 nrm = np.cross(du.T, dv.T).T
60 nrm /= np.linalg.norm(nrm, axis=0)
61 blocks.append(np.vstack([xyz, du, dv, nrm]))
62 srcvals = np.hstack(blocks)
63 return fb.Surface(norder, iptype, srcvals)
66def superellipsoid_volume(a, p):
67 """Exact volume of |x|^p+|y|^p+|z|^p <= a^p (used to check the discretization)."""
68 from math import gamma
69 return 8 * a**3 * gamma(1 + 1 / p) ** 3 / gamma(1 + 3 / p)