Minimal Surfaces & Catenaries

1import numpy as np
2
3def solve_minimal_surface_dirichlet(vertices, fixed_boundary_indices, L_cot, max_iter=30):
4    """
5    Solves Plateau's problem for minimal surfaces (zero mean curvature H = 0).
6    Minimizes Dirichlet energy E(V) = 0.5 * Tr(V^T * L_cot * V) by iteratively
7    relaxing interior vertices to the cotangent-weighted average of their 1-ring neighbors.
8    """
9    V = np.array(vertices, dtype=float, copy=True)
10    fixed_set = set(fixed_boundary_indices)
11    free_ids = [i for i in range(len(V)) if i not in fixed_set]
12
13    for _ in range(max_iter):
14        for i in free_ids:
15            # Off-diagonal entries of L_cot hold negative cotangent weights -w_ij
16            row = L_cot[i].copy()
17            diag = row[i]
18            row[i] = 0.0
19            if abs(diag) > 1e-12:
20                V[i] = -np.dot(row, V) / diag
21    return V

1import math
2
3def catenary_curve(x_span=4.0, sag=1.5, num_nodes=25):
4    """
5    Computes exact hyperbolic cosine catenary curve y(x) = a * (cosh(x / a) - 1)
6    by solving a * (cosh(span / (2a)) - 1) = sag via Newton-Raphson iteration.
7    Inverted (-y), this produces Hooke & Gaudi's pure compression arch (M = 0).
8    """
9    half_l = 0.5 * x_span
10    # Initial guess for catenary parameter a = H / w from parabolic approximation
11    a = (half_l ** 2) / (2.0 * max(sag, 1e-5))
12    for _ in range(15):
13        u = half_l / a
14        f_val = a * (math.cosh(u) - 1.0) - sag
15        df_da = math.cosh(u) - 1.0 - u * math.sinh(u)
16        a -= f_val / (df_da - 1e-12)
17
18    y_anchor = a * (math.cosh(half_l / a) - 1.0)
19    coords = []
20    for i in range(num_nodes):
21        x = -half_l + (i / (num_nodes - 1)) * x_span
22        y_hang = a * (math.cosh(x / a) - 1.0) - y_anchor
23        coords.append((x, y_hang, -y_hang))  # (x, hanging_y, inverted_arch_y)
24    return coords