# Cubic-Hermite FEM for a Fourth-Order BVP (Euler–Bernoulli Beam) — Task Specification

You are working on a Linux coding task about `C¹` finite element methods for a
fourth-order two-point boundary value problem.

## Task directory

- `base`

## Visible input files

- Problem specification: `base/input/problem_spec.md`
- Python runtime manifest: `base/input/runtime_env/pyproject.toml`
- Python lockfile: `base/input/runtime_env/uv.lock`
- Canonical Python wrapper: `base/software/python.sh`

## Recommended workflow

1. Change into `base`.
2. Read `base/input/problem_spec.md` in full before coding.
3. Use `"base/software/python.sh" your_script.py` or
   `"base/software/python.sh" -c "..."` to run Python with the staged
   NumPy/SciPy environment.
4. Implement your solver and save it to `base/output/fem_solver.py`.
5. Run the solver and write every required artifact under `base/output`.

## Required outputs under `base/output`

- `fem_solver.py`
- `reference_arrays.npz`
- `results.json`

---

## Background

Implement a **C¹ cubic-Hermite finite element method** for the clamped
Euler–Bernoulli beam (a fourth-order two-point boundary value problem):

$$
u''''(x) = f(x), \quad 0 < x < 1, \qquad u(0)=u'(0)=u(1)=u'(1)=0 .
$$

The variational (weak) form is

$$
(u'',\,v'') = (f,\,v) \quad \text{for all admissible } v,
$$

so the bilinear form pairs **second** derivatives — a conforming method needs
globally `C¹` elements, not the usual `C⁰` Lagrange elements.

## Discretization

Uniform mesh `x_j = j h`, `h = 1/n`, `j = 0,…,n`. On each element use the cubic
Hermite basis with two degrees of freedom per node — a **value** DOF and a
**slope** (first-derivative) DOF. On the reference interval `ξ ∈ [0,1]`:

- value functions `H₁(ξ)=2ξ³−3ξ²+1`, `H₂(ξ)=−2ξ³+3ξ²`;
- slope functions `H₁*(ξ)=ξ³−2ξ²+ξ`, `H₂*(ξ)=ξ³−ξ²`, with the corresponding
  global slope basis scaled by the element length so that each slope DOF equals
  the true derivative `u'(x_j)`.

The clamped boundary conditions remove both DOFs at `x=0` and `x=1`, leaving
`2(n−1)` unknowns.

## Manufactured test problem

$$
u_e(x) = e^{x}\,x^2(1-x)^2, \qquad
f(x) = u_e''''(x) = e^{x}\,(x^4 + 14x^3 + 49x^2 + 32x - 12).
$$

## What to compute

### (A) Assembly and solve

Derive the element stiffness matrix for `(u'',v'')` and the load vector,
assemble the global `2(n−1)` system under the clamped BCs, and solve it. No FEM
library — assemble and solve yourself.

Integrate the load with a rule that is **at least 3-point Gauss–Legendre per
element** (or an exact integral). A cruder rule perturbs the discrete solution
at the nodes and destroys the nodal superconvergence that part (B) asks you to
exhibit; two-point Gauss–Legendre is *not* accurate enough here.

### (B) Two error measures

For `n = 2^p`, `p = 5,6,…,11`, compute the `L∞` error of the finite-element
solution in **two** ways and the corresponding orders of convergence
`β = log₂(e_{n} / e_{2n})`:

  - **nodal:** `e_n = max_j |u_e(x_j) − u_h(x_j)|` over the mesh nodes;
  - **dense:** `e*_n = max over a 10×-finer point set` `|u_e(y) − u_h(y)|`,
    `y_i = i/(10n)`, `i = 0,…,10n`, where `u_h(y)` is the finite-element
    solution evaluated on the element that contains `y`.

Tabulate `(n, h, e_n, β)` and `(n, h, e*_n, β*)`. Determine the asymptotic order
of the method, and report what each measure does as `n` grows — they do **not**
behave the same way; explain what you observe.

### (C) Artifact

At `n = 256`, save the dense reconstruction of `u_h` over the `10·256 + 1 = 2561`
points `y_i = i/2560`, together with the exact solution and the pointwise error.

## Required solver interface

`fem_solver.py` must be importable **without side effects** (guard any demo code
behind `if __name__ == "__main__":`) and must expose exactly this entry point:

```python
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)
    """
```

The grader imports your module in an isolated subprocess and calls `solve()` to
regenerate the study, then checks that the regenerated errors — and the `n = 256`
dense field itself — match what you submitted. Arrays your code cannot reproduce
will not be accepted, so do not hand-edit the outputs.

Both returned grids must be the prescribed ones: `x_nodes[j] = j/n` and
`x_dense[i] = i/(10n)`, to within `1e-15` absolute (i.e. exactly, up to
double-precision representation).

## Required outputs (write to `base/output`)

| File | Format | Contents |
|------|--------|----------|
| `fem_solver.py` | Python | Your implementation, exposing `solve(n)` as above. |
| `reference_arrays.npz` | NPZ | The two convergence studies and the `n = 256` artifact — see the array schema below. |
| `results.json` | JSON | `nodal_study`, `dense_study`, `asymptotic_order`, `artifact`, and the metadata below. Must be consistent with the arrays. |

### Scheme fidelity

The grader carries its **own** cubic-Hermite implementation and compares your
dense field against it pointwise at **every mesh** `n = 32 … 2048` (tolerance
`1e-9`). Doing genuine finite elements on some meshes and fabricating the rest
fails. Your returned `u_nodes` must also be your own `u_dense` sampled at the
nodes (they come from the same degrees of freedom, so they agree exactly);
substituting exact values at the nodes to improve the apparent superconvergence
fails.
Implementations that differ only in a legitimate detail — say 5-point instead of
3-point Gauss — agree to about `5e-12`, so this rejects nothing valid.

Field agreement alone is not proof of assembly for this special problem: the
exact cubic-Hermite interpolant can agree with the Ritz field to round-off. The
grader therefore also observes a linear solve of the published `2(n−1)` free-DOF
system at every mesh. At `n = 32` it replaces intercepted solve results with
non-finite values and requires the returned field to change, become non-finite,
or raise (`atol=1e-12` for the unchanged-field test). A discarded dummy solve
therefore fails. A hand-written linear solver that cannot be instrumented is
flagged for manual source review (and scores `0.0` under `--strict`).

### Mandatory regeneration and value types

The grader re-runs your solver and must be able to reproduce **every** mesh —
nothing is best-effort. The reference solver completes all seven meshes in well
under a second, so there is no legitimate reason a conforming submission cannot
be re-run.

| regenerated | if it cannot be reproduced |
|---|---|
| `n = 32, 64, 128, 256` (includes the artifact mesh) | 0.0 |
| `n = 512, 1024, 2048` | 0.0 |
| the `n = 256` dense field, compared **pointwise** | Tier 3 fails → 0.5 |

Re-running the same code in the same pinned environment reproduces the dense
field bitwise, so a pointwise deviation above `1e-11` costs Tier 3 even when the
maximum error still looks correct. Do not hand-edit or re-use a field from
another run.

Every mesh size `n`, in both `reference_arrays.npz` and `results.json`, must be an
**exact integer** (`32`, not `32.0000001` or `32.25`). Non-integer values are
rejected outright rather than silently truncated.

### Array schema (`reference_arrays.npz`)

| key | shape | contents |
|---|---|---|
| `n_values` | `(7,)` | `2^p`, `p = 5..11`, ascending: `32 … 2048` |
| `nodal_errors` | `(7,)` | `e_n` for each mesh, same order as `n_values` |
| `dense_errors` | `(7,)` | `e*_n` for each mesh, same order as `n_values` |
| `nodal_eoc` | `(7,)` | `log₂(e_{n/2}/e_n)`; entry 0 has no predecessor — use `NaN` |
| `dense_eoc` | `(7,)` | `log₂(e*_{n/2}/e*_n)`; entry 0 is `NaN` |
| `x_dense` | `(2561,)` | the artifact grid, `x_dense[i] = i/2560` |
| `uh_dense` | `(2561,)` | the `n = 256` FE field on that grid |
| `uexact_dense` | `(2561,)` | `u_e` on the same grid |
| `error_dense` | `(2561,)` | `uh_dense − uexact_dense` |

`x_dense` must be the exact grid `i/2560` (to `1e-15` absolute), starting at
exactly `0.0` and ending at exactly `1.0`.

`max|uh_dense − uexact_dense|` must equal `dense_errors[3]`, the `n = 256` row of
the dense study: the artifact and the convergence table must come from the same
run.

### `results.json` schema

| key | contents |
|---|---|
| `problem`, `exact_solution`, `f` | short strings describing the BVP and the manufactured data |
| `dofs_per_node` | `["value", "slope"]` |
| `nodal_study.table` | 7 rows of `n`, `h`, `e_n`, `beta` (`beta` is `null` for `n = 32`) |
| `dense_study.table` | 7 rows of `n`, `h`, `e_star_n`, `beta` (`beta` is `null` for `n = 32`) |
| `dense_study.observed_order_clean_regime` | the order fitted over the pre-round-off meshes `n = 32…256` |
| `asymptotic_order` | the asymptotic order of the method in the dense norm |
| `artifact` | `n` (256), `points` (2561), `dense_max_error` |
| `notes` | your explanation of how the two error measures differ as `n` grows |

`h` must equal `1/n`, and every `beta` must equal the order recomputed from the
errors you submitted.

## Rules

- **Python + NumPy + SciPy only.** Do **not** use a FEM/PDE package (FEniCS,
  scikit-fem, sfepy, FiPy, …) — build the element matrices, the global assembly,
  and the solve yourself (`scipy.sparse`, `scipy.linalg`). The grader enforces
  this with an import allowlist over `fem_solver.py`.
- Use the manufactured `f` and the clamped boundary conditions exactly as given,
  the uniform mesh `h = 1/n`, the mesh sequence `n = 2^p`, `p = 5..11`, and the
  dense grid `y_i = i/(10n)`.
- All saved values must be truthful and internally consistent — the grader
  recomputes the exact solution and every convergence order from your arrays,
  and **re-runs your solver** to regenerate the study. If you cannot complete the
  `n = 256` artifact, still write honest study outputs rather than fabricating it.
- Do not modify files outside `base/output`. (Running the staged
  `python.sh` wrapper materialises its virtual environment in a temporary
  location **outside** the task directory — that is expected and permitted, and
  is not a modification of task files. The staged inputs are left untouched, so
  they may be mounted read-only.)