[michakraus]

A parameter square whose left radial edge is a pole — one point of the physical domain reached from every angle, as on a mapped disk — carries no tensor-product spline space that is even C⁰ there. The geometry map collapses the whole circle s = 0 to a point, and nothing constrains a tensor-product basis's θ-dependence on it. This is a stronger obstruction than the loss of C¹ one expects from a polar singularity, and no refinement removes it.

PolarSplineBasis replaces the first two radial rows — 2·Nθ functions — with three, the pole triangle, spanning the constants and the two linear functions of a pseudo-Cartesian chart at the pole. The space is C⁰ and C¹ there by construction, not by refinement. Construction after Toshniwal, Speleers, Hiemstra & Hughes, CMAME 316 (2017) 1005–1061, doi:10.1016/j.cma.2016.11.009 — the reference Zoni & Güçlü's disk solver is built on, read from their §3.3 rather than from an abstract.

What it looks like from outside

A pole couples the two axes, so it is not a boundary condition and is not reachable through BSplineBasis(mesh, p, bc). The tensor-product layer's independent-per-axis promise is untouched, because a polar space is not a tensor product.

The polar space is the parent conjugated by a sparse recombination_matrix, exactly as a RecombinedBSplineBasis is one dimension down — so there is no second quadrature rule, and every assembly is R' · (the parent's) · R.

Two things do not carry over, both consequences of a pole function reaching around the whole angular axis:

• The index set is not a product. Coefficients are a vector, not an array, and evaluate_all returns the indices of the nonzero block rather than its first index, because the block is not contiguous. parent_coefficients is the bridge back to the parent's array shape.
• There is no KroneckerMass. The mass matrix is assembled and factorised — a sparse Cholesky, with only a 3×3 dense block — instead of one one-dimensional solve per axis.

Six new exports: PolarSplineBasis, PolarSplineQuadrature, pole, pole_triangle, pseudo_cartesian, parent_coefficients.

The chart, which is the one thing worth arguing about

C(θ) = Σⱼ cos θⱼ Mⱼ and S(θ) = Σⱼ sin θⱼ Mⱼ are splines, not the trigonometric functions: cos θ is not in the angular spline space and no choice of coefficients puts it there. The chart is therefore built from them, x̃ = (s−a)C(θ), ỹ = (s−a)S(θ), and in it the C¹ property is exact rather than approximate — the space contains 1, x̃ and ỹ to round-off. Its unit circle sits some 2.5 % inside the true one at 16 angular cells, closing at the angular order. A geometry map that is itself represented in this space carries the smoothness to the physical domain.

Measurements

Four scripts under scripts/. The three that check a property exit non-zero when a check or a control fails. polar_mass_cost.jl is a timing report: it has no threshold to fail against and always exits zero.

Continuity, 32 angles, scripts/polar_continuity.jl:

	polar	control	control is
C⁰ spread, worst single basis function	4.996e-16	2.320e+00	the plain TensorProductBasis
C¹, relative residual of one gradient fitted over 32 angles	4.444e-16	1.022e+00	a space imposing C⁰ and nothing more
C⁰ spread of that same C¹ control	—	4.441e-16	it must still pass C⁰

C⁰ is measured on every basis function, not only a random combination: a spread hiding in one function and cancelling in a combination would pass. The C¹ gradient is fitted by least squares against the chart directions, and the radial derivative is taken analytically — a difference quotient would measure the limit and the round-off together.

Dropping one of the three pole functions is not a C¹ control. It leaves a subspace of the polar space, and a subspace of a C¹ space is C¹: measured at 5.329e-15, the same as the full space. It is a must-fail control for completeness, and is used as one — the relative L² error of the projection of ỹ on the pole cells goes from 5.442e-15 to 5.952e-02.

Approximation order, cubic on both axes, L² projection of a smooth function of the disk:

ns × nθ	N	‖e‖ in ds dθ	order	‖e‖ in s ds dθ	order
8 × 16	147	1.1213e-03	—	1.0593e-03	—
16 × 32	547	4.0082e-05	4.806	3.7544e-05	4.818
32 × 64	2115	2.1442e-06	4.224	2.0022e-06	4.229
64 × 128	8323	1.2865e-07	4.059	1.2002e-07	4.060

Full order p+1. The error is not concentrated at the pole — it is smaller there: per-radial-cell RMS residual 5.37e-10 and 4.28e-10 in the two pole cells against 1.34e-07 for the worst cell elsewhere.

Partition of unity and the mapped measure. max |Σₖ Ψₖ − 1| = 1.554e-15 over 176 points including the pole, checked through both evaluate_all and the sum over every index. Least value of a pole function 0.000e+00 — non-negative, which is what the vertex radius buys. 𝟙ᵀ𝕄𝟙 = 2π to 6.217e-15. With the Jacobian of a mapped disk through weighted_matrix, 𝟙ᵀ𝕄_J𝟙 = 114.776991 against an independently computed 114.777, relative difference 8.058e-08; the analytic Jacobian is itself checked against a central difference of the map to 3.171e-09.

What the pole costs, 64×128, N = 8323: assembly 41 ms and Cholesky 16 ms, both once per space; one mass solve 0.74 ms against 0.136 ms for the KroneckerMass solve at the same mesh. A factor of 5.4, under a millisecond, and not where the time of a nonlinear flow goes.

Pre-PR verification

Fixed during the pass

• evaluate_all presizes its index and value vectors. 1184 → 640 B and 167 → 125 ns per call; a deposition loop over 10 000 points goes 1152 → 608 B per point. Output bit-identical.
• The mass matrix is symmetrised in the constructor, as SplineQuadrature already does, so issymmetric(mass_matrix(q)) is true and matches the docstring and both sibling quadratures. nnz unchanged; largest entry moves by 6.9e-18. The one-off space build goes 106 → 117 ms.
• Two assertions of the form @test issymmetric(Symmetric(X)) could not fail, Symmetric(X) being symmetric for every X. The mass matrix now asserts exact symmetry; the stiffness matrix, which is not symmetrised, asserts symmetry to round-off.
• The Ns ≥ 3 constructor guard was unreachable — a clamped basis has ncells + p functions and both are already bounded below — so it is now the invariant as a comment plus a test that measures it.
• The CHANGELOG's Cholesky timing did not reproduce: 115 ms was one loaded run, three cold runs give 16.0 / 16.7 / 16.0 ms, and the solve ratio is 5.4 rather than 5.6. Every other figure in that line reproduced to the quoted digit.

Checked and clean

• Type piracy: Aqua.Piracy.hunt 0 hits; Aqua.test_all(…; persistent_tasks = false) 10/10.
• Inference: PolarSplineQuadrature is concrete in all four parameters and all seven fields; basis_values, mixed_matrix, weighted_matrix, evaluate, evaluate_all, parent_coefficients, pseudo_cartesian, l2_projection and both constructors all return concrete types, with no Any local, Union split or box under code_warntype.
• No duplication of what the dependencies or the sibling modules provide; recombination_matrix is extended rather than reimplemented.
• Comments and docstrings are present tense and describe the code as it is; no historical notes outside CHANGELOG.md.
• 64 895 tests pass, the 391 new polar ones included. Docs build and doctests clean. JuliaFormatter, fatou lint and the NFC check clean.

Pre-existing, reported not fixed

• evaluate(::TensorProductBasis, i::Integer, x, d) allocates 368 B per call.
• evaluate(B, k, x) for one of the three pole functions costs 11 280 B / 4.6 µs against 368 B / 0.21 µs for any other, being O(Nθ) parent evaluations. The docstring states this; reducing it to O(p²) is an API change and belongs in its own PR.
• Aqua's persistent_tasks check was red on every test job here until Pin Aqua to 0.8.16 #16 landed. Aqua v0.8.17 walks each project's deps and weakdeps and raises when Base.locate_package cannot find one, and an uninstalled weak dependency — ChainRulesCore behind AbstractFFTs, reached through FFTW — is the normal case. Pin Aqua to 0.8.16 #16 pinned Aqua to =0.8.16; this branch is rebased onto it and green.

🤖 Generated with Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[michakraus]

CI here is red on every test job for a reason that is not this branch: Aqua v0.8.17, released 2026-09-17 at 16:46 UTC, a few hours before this PR was opened.

Aqua#394 added a manifest walk that iterates each project's deps and weakdeps and raises when Base.locate_package cannot find one. ChainRulesCore is a weak dependency of AbstractFFTs, which this package reaches through FFTW, and it is not loaded in the test process — so persistent_tasks fails on a healthy environment. Aqua.test_all stops at its first failing check, so it took the rest of runtests.jl with it.

main was last green on 2026-09-15, before that release, and this branch changes no dependency — git diff origin/main...HEAD -- Project.toml is empty. Aqua.test_all(SimpleSplines; persistent_tasks = false) passes 10/10 here and the rest of the suite is 64 895 passing.

#16 pins Aqua to exactly 0.8.16 and is green. Once it lands I will rebase this branch onto it and rerun the checks.

[michakraus]

#16 is green on all ten checks — every platform, every Julia version, min through nightly. That confirms Aqua v0.8.17 was the sole cause of the red here: the only difference between the two branches on the main baseline is the [compat] bound.

This branch stays red until #16 lands. Two ways to get there, whichever you prefer:

• Merge Pin Aqua to 0.8.16 #16 first, then I rebase this branch onto the new main and rerun the checks.
• Retarget this PR's base to compat/pin-aqua-0-8-16 and merge the pin into it, which makes it green now and turns the two into a chain. GitHub retargets this one to main automatically when Pin Aqua to 0.8.16 #16 merges.

I have not done either, since both start with merging something.

[michakraus]

Review: approve

The construction is derivable from its own docstring, every quoted figure reproduces in a cold process, and the required checks are green. Seven findings, none of them a bug.

CI

All 7 required contexts pass: Julia 1 and Julia min on ubuntu/macOS/windows, plus Doctests - ubuntu-latest. Julia nightly, Julia pre and Documentation also pass and are not required. The persistent_tasks fault is gone — the branch sits on the Aqua =0.8.16 pin, and a local run resolved 0.8.16 and passed that check in 7.6 s.

The mathematics

Checked by hand against the docstring. Barycentric ℓ_k for vertices v_k = r(cos ψ_k, sin ψ_k) centred on the pole has gradient 2/(3r). The constraint at src/polar.jl:184 needs gradient N₂'(a)/3. Equating gives r = 2/N₂'(a), which is src/polar.jl:181. Partition of unity follows from Σ_k (1 + cos(ψ_k − θ_j))/3 = 1 for the three-fold ψ. The code is consistent with the derivation it states.

Every measurement reproduces

Cold process at the PR head:

claim	in the PR	measured
mass solve, 64×128, N = 8323	0.74 ms	0.736 ms
KroneckerMass solve, same mesh	0.136 ms	0.136 ms
Cholesky	16 ms	16.05 ms
assembly	41 ms	38.55 ms
approximation order, finest pair	4.059 / 4.060	4.059 / 4.060
pole-cell RMS vs worst bulk cell	5.37e-10, 4.28e-10 vs 1.34e-07	identical
C⁰ spread, worst basis function	4.996e-16	4.996e-16
C⁰ tensor-product control	2.320e+00	2.320e+00
C¹ residual / control	4.444e-16 / 1.022	identical
Ψ₃-dropped control, pole cells	5.952e-02	5.952e-02
𝟙ᵀ𝕄𝟙 = 2π	6.217e-15	6.217e-15
mapped 𝟙ᵀ𝕄_J𝟙 vs 114.777	8.058e-08	8.058e-08

Assembly is the one figure that moved, and it is a timing on a loaded machine. Partition of unity measured 5.55e-16 over an independent 176-point grid.

The controls are real. In scripts/polar_continuity.jl:253 each control's failure is a conjunct of the exit condition, so a control that quietly starts passing turns the script red. The script is also honest about the one control that is not a continuity control: dropping Ψ₃ leaves a subspace of a C¹ space, so it is used as a completeness control, and lines 20-22 say so.

Findings

#	severity	file:line	finding
1	quality	src/polar.jl:294-301	pseudo_cartesian rebuilds cos.(θ) and sin.(θ) over all Nθ Greville abscissae on every call, so it is O(Nθ) in time and allocation per point: 3584 B / 1167 ns at 64×128 cubic, against 368 B / 125 ns for an ordinary evaluate. Both test/polar_tests.jl:273 and docs/src/theory/polar.md call it once per quadrature point, ~65 000 calls at that mesh. The two vectors are fixed for the basis and belong in the struct.
2	quality	scripts/polar_mass_cost.jl	The PR body says "four scripts under scripts/, each exiting non-zero when a check or a control fails". This one has no assertion, no pass flag and no exit — it prints a table. A timing script needs no threshold; the sentence describing it is what is wrong.
3	nit	src/polar.jl:402-406	evaluate(B, û, X::AbstractVector) is the vector-of-points method, so a single point passed as [0.3, 1.0] rather than (0.3, 1.0) fails with a bare BoundsError. The PR body files this under "pre-existing", but the method is new in this diff — TensorProductBasis has four evaluate methods and none of them takes a vector of points, and the 1-D form cannot hit this because a 1-D point is a scalar.
4	nit	src/polar.jl:673, src/polar.jl:762	[f(pt) for pt in Iterators.product(x...)][:] allocates the grid array and then copies it. vec is free, and weighted_matrix is on the per-Newton-iteration path its own docstring names.
5	nit	src/polar.jl:494	The signature line reads nq = quadrature_order.(degree(B)); the code at line 544 is map(quadrature_order, degree(B)). Equivalent, but the rendered signature does not match the source.
6	nit	CHANGELOG.md:96-97	Two consecutive blank lines before ## [0.2.0].
7	nit	src/polar.jl:326-332	Base.show(::PolarSplineBasis) has no test, and neither do the meshwidth, mesh and domain forwards at lines 314-316.

Checked and clean

• Exported API. All six exports carry docstrings and all six appear in docs/src/library.md:143-150. docs/make.jl:14 sets checkdocs = :exports and the Documentation job is green, so no exported name is undocumented. Tests in test/polar_tests.jl cover each one, and six jldoctest blocks run under the green Doctests job. Both new .bib entries are cited from docs/src/theory/polar.md:152-153.
• Inference. evaluate → Float64, evaluate_all → Tuple{Vector{Int64},Vector{Float64}}, pseudo_cartesian → Tuple{Float64,Float64}, parent_coefficients → Matrix{Float64}, l2_projection → Vector{Float64}. Both structs and all their fields are concrete.
• Aqua, run per check so a first failure could not hide the rest: ambiguities, unbound_args, undefined_exports, project_extras, stale_deps, deps_compat, piracies, persistent_tasks all pass on 0.8.16.
• No piracy by inspection either. Every Base extension in src/polar.jl dispatches on PolarSplineBasis or PolarSplineQuadrature. Base.length/eachindex without Base.iterate matches src/tensorproduct.jl:101,105.
• CHANGELOG describes the end state, names all six exports, and every figure in it reproduced. It does not narrate the commit sequence.
• JuliaFormatter and the NFC check are clean on all changed files.

Not checked

• The "64 895 tests pass" count. Nine green CI jobs establish that the suite is green; the count itself is not worth a cold run to confirm.
• The reference value 114.777 itself. The comparison runs and agrees to 8.1e-08; whether the reference is right is a separate question.
• Whether the O(Nθ) pole-function evaluation can be reduced. At Nθ = 131 one such evaluation is 90 128 B / 30 µs against 368 B / 125 ns for any other, so the cost is real and grows — but the docstring states it and deferring it is right.

One non-blocking note outside this diff: seven non-public names are imported (@simplify and QMul2 from ContinuumArrays; basis, degree, nnodes, nodes, order from GeometricBase) and two accessed (FFTW.ESTIMATE, FFTW.UNALIGNED). No hit lies in a hunk this branch touched.