codexpast-exam-of-the-geometry-of-surfaces-course-of-the-university-of-oxford2023-dihedral-free-action.py
past-exam-of-the-geometry-of-surfaces-course-of-the-university-of-oxford/2023-dihedral-free-action.py
"""Original drilled-torus figure; Python 3.14.4, NumPy 2.3.5, Matplotlib 3.10.7.
Uses only the dependencies in the wiki's root pyproject.toml. Run from the
desired output directory; the PNG and mesh-check JSON are written there.
The mesh is a marching-tetrahedra approximation to a smooth maximum of a
solid-torus defining function and three cylindrical bore defining functions.
"""
import json
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
N = 3
R = 2.6
TUBE_RADIUS = 0.8
BORE_RADIUS = 0.35
SMOOTHNESS = 60.0
AXIS_ANGLES = np.arange(N) * np.pi / N
def field(x, y, z):
"""Negative inside the smoothed solid; positive on every rotation axis."""
value = (np.hypot(x, y) - R) ** 2 + z * z - TUBE_RADIUS**2
value = SMOOTHNESS * value
for theta in AXIS_ANGLES:
transverse = -x * np.sin(theta) + y * np.cos(theta)
bore = BORE_RADIUS**2 - transverse**2 - z**2
value = np.logaddexp(value, SMOOTHNESS * bore)
return value / SMOOTHNESS
def extract_mesh():
# A fixed diagonal in every cube makes adjacent tetrahedra conforming.
nx, ny, nz = 137, 137, 49
xs = np.linspace(-3.65, 3.65, nx)
ys = np.linspace(-3.65, 3.65, ny)
zs = np.linspace(-1.0, 1.0, nz)
xx, yy, zz = np.meshgrid(xs, ys, zs, indexing="ij")
xyz = np.column_stack((xx.ravel(), yy.ravel(), zz.ravel()))
values = field(xx, yy, zz).ravel()
ii, jj, kk = np.meshgrid(
np.arange(nx - 1), np.arange(ny - 1), np.arange(nz - 1), indexing="ij"
)
bases = ((ii * ny + jj) * nz + kk).ravel()
offsets = np.array([0, ny * nz, nz, ny * nz + nz,
1, ny * nz + 1, nz + 1, ny * nz + nz + 1])
tetrahedra = [(0, 1, 3, 7), (0, 3, 2, 7), (0, 2, 6, 7),
(0, 6, 4, 7), (0, 4, 5, 7), (0, 5, 1, 7)]
all_triangle_edges = []
for tet in tetrahedra:
ids = bases[:, None] + offsets[list(tet)]
inside = values[ids] < 0
patterns = inside @ np.array([1, 2, 4, 8])
for pattern in range(1, 15):
selected = ids[patterns == pattern]
if not len(selected):
continue
inn = [i for i in range(4) if pattern & (1 << i)]
out = [i for i in range(4) if not pattern & (1 << i)]
if len(inn) in (1, 3):
if len(inn) == 1:
local_edges = [(inn[0], j) for j in out]
else:
local_edges = [(out[0], j) for j in inn]
all_triangle_edges.append(selected[:, local_edges])
else:
a, b = inn
c, d = out
# Intersections a-c, a-d, b-d, b-c form the quadrilateral.
for local_edges in [[(a, c), (a, d), (b, d)],
[(a, c), (b, d), (b, c)]]:
all_triangle_edges.append(selected[:, local_edges])
triangle_edges = np.concatenate(all_triangle_edges)
triangle_edges.sort(axis=2)
edges, inverse = np.unique(triangle_edges.reshape(-1, 2), axis=0,
return_inverse=True)
weights = values[edges[:, 0]] / (values[edges[:, 0]] - values[edges[:, 1]])
vertices = (xyz[edges[:, 0]] * (1 - weights[:, None])
+ xyz[edges[:, 1]] * weights[:, None])
faces = inverse.reshape(-1, 3)
return vertices, faces
def mesh_checks(vertices, faces):
mesh_edges = np.concatenate((faces[:, [0, 1]], faces[:, [1, 2]],
faces[:, [2, 0]]))
mesh_edges.sort(axis=1)
unique_edges, counts = np.unique(mesh_edges, axis=0, return_counts=True)
parents = np.arange(len(vertices))
def root(a):
while parents[a] != a:
parents[a] = parents[parents[a]]
a = parents[a]
return a
for a, b in unique_edges:
parents[root(a)] = root(b)
components = len({root(i) for i in range(len(vertices))})
chi = len(vertices) - len(unique_edges) + len(faces)
# Check invariance analytically at numerical samples, not just on the mesh.
rng = np.random.default_rng(2023)
p = rng.uniform(-3.4, 3.4, (1000, 3))
angle = 2 * np.pi / N
beta = np.array([[np.cos(angle), np.sin(angle), 0],
[-np.sin(angle), np.cos(angle), 0], [0, 0, 1]])
alpha = np.diag([1, -1, -1])
errors = []
for matrix in (alpha, beta):
q = p @ matrix.T
errors.append(float(np.max(np.abs(field(*p.T) - field(*q.T)))))
t = np.linspace(-4, 4, 2001)
axis_minima = [float(np.min(field(0 * t, 0 * t, t)))]
for angle in AXIS_ANGLES:
axis_minima.append(float(np.min(field(t * np.cos(angle),
t * np.sin(angle), 0 * t))))
checks = {"vertices": len(vertices), "edges": len(unique_edges),
"faces": len(faces), "euler_characteristic": chi,
"genus": (2 - chi) // 2, "components": components,
"all_edges_have_two_incident_faces": bool(np.all(counts == 2)),
"generator_invariance_errors": errors,
"minimum_field_on_rotation_axes": axis_minima}
assert components == 1 and chi == -12 and np.all(counts == 2), checks
assert max(errors) < 1e-12 and min(axis_minima) > 0, checks
return checks
def make_figure(vertices, faces):
plt.rcParams.update({"font.size": 10, "figure.facecolor": "white"})
fig = plt.figure(figsize=(10.4, 5.1), layout="constrained")
ax = fig.add_subplot(1, 2, 1, projection="3d")
triangles = vertices[faces]
normals = np.cross(triangles[:, 1] - triangles[:, 0],
triangles[:, 2] - triangles[:, 0])
normals /= np.linalg.norm(normals, axis=1)[:, None]
light = np.array([0.3, -0.5, 0.8])
light /= np.linalg.norm(light)
shade = 0.42 + 0.58 * np.abs(normals @ light)
colors = np.column_stack((0.28 * shade, 0.68 * shade, 0.85 * shade,
np.ones(len(shade))))
mesh = Poly3DCollection(triangles, facecolors=colors, edgecolors="none",
linewidth=0, rasterized=True)
ax.add_collection3d(mesh)
colors_axes = ["#b4462e", "#7a4bab", "#277341"]
for theta, color in zip(AXIS_ANGLES, colors_axes):
t = np.linspace(-3.7, 3.7, 2)
ax.plot(t * np.cos(theta), t * np.sin(theta), 0 * t,
color=color, linestyle="--", linewidth=1.0)
ax.plot([0, 0], [0, 0], [-1.6, 1.8], color="#333333",
linestyle="--", linewidth=1.2)
ax.text(0, 0, 1.92, "vertical rotation axis", ha="center", fontsize=9)
ax.set(xlim=(-3.7, 3.7), ylim=(-3.7, 3.7), zlim=(-1.6, 2.0))
ax.set_box_aspect((1, 1, 0.5))
ax.view_init(elev=27, azim=-76)
ax.set_axis_off()
ax.set_title("Smooth drilled torus, n = 3\ngenus = 1 + 6 = 7", pad=1)
plan = fig.add_subplot(1, 2, 2)
q = np.linspace(-3.7, 3.7, 701)
xx, yy = np.meshgrid(q, q)
ff = field(xx, yy, 0 * xx)
plan.contourf(xx, yy, ff, levels=[-100, 0], colors=["#55a1be"])
plan.contour(xx, yy, ff, levels=[0], colors=["#1d556b"], linewidths=0.8)
for j, (theta, color) in enumerate(zip(AXIS_ANGLES, colors_axes)):
plan.plot(q * np.cos(theta), q * np.sin(theta), linestyle="--",
color=color, linewidth=1.1, label=f"half-turn axis {j+1}")
plan.plot(0, 0, "o", color="#333333", markersize=4)
plan.text(0.12, 0.14, "vertical axis", fontsize=9)
plan.annotate("radial bore", xy=(R, 0), xytext=(3.65, -1.25),
arrowprops={"arrowstyle": "->", "color": "#333333"},
ha="right", fontsize=9)
plan.set_aspect("equal")
plan.set(xlim=(-3.85, 3.85), ylim=(-3.85, 3.85))
plan.axis("off")
plan.set_title("Horizontal section at z = 0\nall rotation axes avoid the surface")
plan.legend(loc="lower center", fontsize=8, frameon=False,
bbox_to_anchor=(0.5, -0.04), ncol=1)
fig.savefig("2023-dihedral-free-action.png", dpi=145, facecolor="white")
plt.close(fig)
if __name__ == "__main__":
vertices, faces = extract_mesh()
checks = mesh_checks(vertices, faces)
make_figure(vertices, faces)
Path("2023-dihedral-free-action-checks.json").write_text(
json.dumps(checks, indent=2) + "\n"
)
print(json.dumps(checks))