"""Gold reference: cubic-Hermite (C^1) finite element method for the Euler-Bernoulli beam (a 4th-order two-point BVP). u''''(x) = f(x), 0 < x < 1, u(0) = u'(0) = u(1) = u'(1) = 0 (clamped both ends). Weak form: (u'', v'') = (f, v) for all admissible v. Manufactured test problem: u_e(x) = e^x x^2 (1-x)^2, f(x) = u_e''''(x) = e^x (x^4 + 14 x^3 + 49 x^2 + 32 x - 12). Discretization: uniform mesh x_j = j h, h = 1/n. Two degrees of freedom per node (value and slope). The slope shape functions are scaled by h so that the slope DOF equals the actual derivative u'(x_j). The clamped boundary conditions remove both DOFs at x = 0 and x = 1, leaving 2(n-1) unknowns. Two error measures at L-infinity: * nodal error e_n = max_j |u_e(x_j) - u_h(x_j)| (nodally superconvergent), * dense error e*_n = max over a 10x-finer point set, where u_h between nodes is reconstructed from the LOCAL element basis -- which REQUIRES the slope DOFs and slope basis functions, not only the values. Dependencies: Python + NumPy + SciPy only. """ import argparse import json import os import shutil import sys import numpy as np from scipy.sparse import lil_matrix, csr_matrix from scipy.sparse.linalg import spsolve # n = 2^p for p = 5..11, as in the original problem P_RANGE = list(range(5, 12)) DENSE = 10 # sub-intervals per element for e*_n (10n points) # 3-point Gauss-Legendre on [0, 1] (5th order) GL_XI = np.array([(1 - np.sqrt(3 / 5)) / 2, 0.5, (1 + np.sqrt(3 / 5)) / 2]) GL_W = np.array([5 / 18, 8 / 18, 5 / 18]) def u_exact(x): return np.exp(x) * (x * (1 - x))**2 def f_rhs(x): return np.exp(x) * (x**4 + 14 * x**3 + 49 * x**2 + 32 * x - 12) def shape(xi, h): """Cubic Hermite shape functions [v_L, theta_L, v_R, theta_R] on ref [0,1]. Slope functions carry the factor h so the slope DOF is the true derivative. """ H1 = 2 * xi**3 - 3 * xi**2 + 1 G1 = h * (xi**3 - 2 * xi**2 + xi) H2 = -2 * xi**3 + 3 * xi**2 G2 = h * (xi**3 - xi**2) return np.array([H1, G1, H2, G2]) def solve_fem(n): """Assemble and solve. Returns (h, U) with U the full 2(n+1) DOF vector.""" h = 1.0 / n ndof = 2 * (n + 1) Ke = (1.0 / h**3) * np.array([ [12, 6 * h, -12, 6 * h], [6 * h, 4 * h * h, -6 * h, 2 * h * h], [-12, -6 * h, 12, -6 * h], [6 * h, 2 * h * h, -6 * h, 4 * h * h], ]) K = lil_matrix((ndof, ndof)) b = np.zeros(ndof) for e in range(n): xL = e * h dofs = [2 * e, 2 * e + 1, 2 * e + 2, 2 * e + 3] be = np.zeros(4) for xi, w in zip(GL_XI, GL_W): be += w * f_rhs(xL + h * xi) * shape(xi, h) be *= h for a in range(4): b[dofs[a]] += be[a] for c in range(4): K[dofs[a], dofs[c]] += Ke[a, c] free = np.arange(2, ndof - 2) # clamp both ends U = np.zeros(ndof) U[free] = spsolve(csr_matrix(K)[free][:, free], b[free]) return h, U def nodal_error(n, U, h): nodes = np.arange(n + 1) * h return float(np.max(np.abs(U[0::2] - u_exact(nodes)))) def dense_reconstruct(n, U, h, dense=DENSE): """Reconstruct u_h on a 10x-finer grid using the local element basis.""" xs, uh = [], [] for e in range(n): xL = e * h coef = U[[2 * e, 2 * e + 1, 2 * e + 2, 2 * e + 3]] for j in range(dense + (1 if e == n - 1 else 0)): xi = j / dense xs.append(xL + h * xi) uh.append(float(coef @ shape(xi, h))) return np.array(xs), np.array(uh) def dense_error(n, U, h, dense=DENSE): x, uh = dense_reconstruct(n, U, h, dense) return float(np.max(np.abs(uh - u_exact(x)))) # --------------------------------------------------------------------------- # Required benchmark entry point (see problem_spec.md, "Required solver # interface"). The grader imports this module in an isolated subprocess and # calls solve() to regenerate the whole study, so it must be side-effect free. # --------------------------------------------------------------------------- def solve(n): """Assemble, solve, and reconstruct the cubic-Hermite FE solution. Parameters ---------- n : int Number of uniform elements; mesh size h = 1/n. Returns ------- (x_nodes, u_nodes, x_dense, u_dense) x_nodes : (n+1,) mesh nodes, x_nodes[j] = j/n u_nodes : (n+1,) FE solution at the nodes (the value DOFs) x_dense : (10n+1,) dense evaluation points, x_dense[i] = i/(10n) u_dense : (10n+1,) FE field at those points, reconstructed from the local element basis using ALL four element DOFs (values AND slopes) """ n = int(n) h, U = solve_fem(n) x_nodes = np.arange(n + 1) * h u_nodes = np.asarray(U[0::2], dtype=float).copy() x_dense, u_dense = dense_reconstruct(n, U, h) return x_nodes, u_nodes, x_dense, u_dense def run(): rows = [] for p in P_RANGE: n = 2**p h, U = solve_fem(n) rows.append((n, nodal_error(n, U, h), dense_error(n, U, h))) return rows def build_outputs(): """Regenerate every required array and the matching JSON sidecar.""" rows = run() n_values = np.asarray([r[0] for r in rows], dtype=np.int64) nodal_errors = np.asarray([r[1] for r in rows], dtype=float) dense_errors = np.asarray([r[2] for r in rows], dtype=float) nodal_eoc = np.r_[np.nan, np.log2(nodal_errors[:-1] / nodal_errors[1:])] dense_eoc = np.r_[np.nan, np.log2(dense_errors[:-1] / dense_errors[1:])] artifact_index = int(np.where(n_values == 256)[0][0]) x_nodes, u_nodes, x_dense, u_dense = solve(256) del x_nodes, u_nodes uexact_dense = u_exact(x_dense) error_dense = u_dense - uexact_dense arrays = { "n_values": n_values, "nodal_errors": nodal_errors, "dense_errors": dense_errors, "nodal_eoc": nodal_eoc, "dense_eoc": dense_eoc, "x_dense": x_dense, "uh_dense": u_dense, "uexact_dense": uexact_dense, "error_dense": error_dense, } def table(errors, eoc, error_key): return [ { "n": int(n), "h": float(1.0 / n), error_key: float(errors[i]), "beta": None if i == 0 else float(eoc[i]), } for i, n in enumerate(n_values) ] clean_fit = float(-np.polyfit(np.arange(4), np.log2(dense_errors[:4]), 1)[0]) results = { "problem": "Cubic-Hermite FEM for Euler-Bernoulli beam u''''=f, clamped BCs", "exact_solution": "exp(x)*x^2*(1-x)^2", "f": "exp(x)*(x^4+14x^3+49x^2+32x-12)", "dofs_per_node": ["value", "slope"], "nodal_study": { "description": "L-inf error at the nodes; nodally superconvergent", "table": table(nodal_errors, nodal_eoc, "e_n"), }, "dense_study": { "description": "L-inf error on a 10x-finer grid via local element basis; order 4 then round-off", "table": table(dense_errors, dense_eoc, "e_star_n"), "observed_order_clean_regime": round(clean_fit, 2), }, "asymptotic_order": 4, "notes": "nodal values superconvergent (~round-off); dense error is the genuine O(h^4) field error and requires the slope DOFs in the reconstruction; round-off (ill-conditioning ~ n^4) dominates for n >= 512.", "artifact": { "n": 256, "points": int(x_dense.size), "dense_max_error": float(dense_errors[artifact_index]), }, } return arrays, results def emit(out_dir): """Write fem_solver.py, reference_arrays.npz and results.json to out_dir.""" out_dir = os.path.abspath(out_dir) os.makedirs(out_dir, exist_ok=True) arrays, results = build_outputs() np.savez(os.path.join(out_dir, "reference_arrays.npz"), **arrays) with open(os.path.join(out_dir, "results.json"), "w") as fh: json.dump(results, fh, indent=2) fh.write("\n") src = os.path.abspath(__file__) dst = os.path.join(out_dir, "fem_solver.py") if os.path.realpath(src) != os.path.realpath(dst): shutil.copyfile(src, dst) print(f"wrote {dst}") print(f"wrote {os.path.join(out_dir, 'reference_arrays.npz')} ({len(arrays)} arrays)") print(f"wrote {os.path.join(out_dir, 'results.json')}") if __name__ == "__main__": ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) ap.add_argument("--emit", metavar="DIR", default=None, help="write fem_solver.py, reference_arrays.npz and " "results.json into DIR") args = ap.parse_args() if args.emit: emit(args.emit) sys.exit(0) print(f"{'n':>6} {'nodal e_n':>14} {'beta':>7} {'dense e*_n':>14} {'beta*':>7}") pn = pd = None for n, en, ed in run(): bn = np.log2(pn / en) if pn else float("nan") bd = np.log2(pd / ed) if pd else float("nan") print(f"{n:>6} {en:>14.5e} {bn:>7.3f} {ed:>14.5e} {bd:>7.3f}") pn, pd = en, ed