MOLE Original Curvilinear GI13 reads from wrong plane

`grad3DCurv` needs the logical eta- and zeta-derivatives evaluated at x-face
positions, but they are computed on y- and z-faces. `GI13` moves them between
face sets. Its index map is wrong.

Feeding `GI13` a field whose values encode their own indices, source y-faces
`(m, n+1, o)` to target x-faces `(m+1, n, o)`:

```
out(1,1,2) <- src(1,4,1) *** wrong zeta-plane
out(2,1,2) <- src(2,4,1) *** wrong zeta-plane
out(3,1,2) <- src(3,4,1) *** wrong zeta-plane
out(4,1,2) <- src(4,4,1) *** wrong zeta-plane
15 of 45 output entries read from the wrong zeta-plane
```

No floating point is involved — this is purely which source entry each output
entry is copied from.

**Cause.** `GI13` builds `kron(speye(n*o), I1)`, walking the source in `n*o`
blocks, but the source has `(n+1)*o` blocks. The zero-column pad is appended at
the end rather than one eta-plane being skipped per zeta-level, so the map
drifts by one plane for every zeta-level. The corrupted fraction grows with `o`.

**Second, smaller defect.** Even where the index is right, `GI13` is a
rectangular identity — nearest-neighbour injection, with the last value
duplicated at the end. Compare the 2-D analogue, which performs the same kind of
staggered shift:

| | `GI2` (2-D) | `GI13` (3-D) |
|---|---|---|
| interior | 4-point average, weights `.25` | `speye(m+1,m)` — nearest neighbour |
| ends | one-sided `.5 .25 -.25` | `I1(end,end)=1` — duplicate last value |

This is why `grad2DCurv` converges and `grad3DCurv` does not.

### Consequence

`f = x^2+y^2+z^2` on a curvilinear grid, distortion amplitude 0.10, `k = 2`:

```
STOCK GI13
n grad rms order L=D*G rms order
13 5.8640e-01 - 4.1488e+00 -
21 5.6031e-01 0.09 4.6321e+00 -0.22
33 5.4836e-01 0.05 5.0521e+00 -0.18
49 5.4146e-01 0.03 5.4626e+00 -0.19
```

The error does not decrease under refinement. It is not boundary contamination
either — trimming 8 layers instead of 2 makes it worse, because the sinusoidal
grid's most distorted region is interior.

Isolating the two halves of the composition on the same grid:

| n | `div3DCurv` alone (`F=(x,y,z)`, div=3) | `grad3DCurv` alone |
|---|---|---|
| 13 | 0.0113 | 0.586 |
| 21 | 0.0047 | 0.560 |
| 33 | 0.0020 | 0.548 |

Divergence converges at second order. The gradient does not converge at all.
On a purely *linear* field `f = x+2y+3z`, where the answer is a constant, at
`n = 21`:

| grid distortion | rms err | max err |
|---|---|---|
| 0.00 (Cartesian) | 0.00000 | 0.00000 |
| 0.05 | 0.0893 | 0.341 |
| 0.10 | 0.390 | 1.834 |
| 0.20 | 16.4 | 270.9 |

### Patch

All six shifts `GI13` handles are the same object: a tensor product of three 1-D
operators, one per axis, each either identity, node-to-centre, or
centre-to-node.

| type | source -> target | xi | eta | zeta |
|---|---|---|---|---|
| `Gn` | y-face -> x-face | c->n | n->c | I |
| `Gc` | z-face -> x-face | c->n | I | n->c |
| `Ge` | x-face -> y-face | n->c | c->n | I |
| `Gcy` | z-face -> y-face | I | c->n | n->c |
| `Gee` | x-face -> z-face | n->c | I | c->n |
| `Gnn` | y-face -> z-face | I | n->c | c->n |

with node-to-centre a two-point average, centre-to-node a two-point average in
the interior plus linear extrapolation at the two ends, assembled as
`kron(Az, kron(Ay, Ax))` to match the xi-fastest flattening. This fixes both
defects and replaces the hand-rolled `kron` chains with something checkable.

`GI13.m` in this directory is a drop-in replacement with the same signature.

```
PATCHED GI13
n grad rms order L=D*G rms order
13 2.6199e-02 - 2.6693e-01 -
21 9.2276e-03 2.04 1.0673e-01 1.79
33 3.4612e-03 2.09 4.4987e-02 1.84
49 1.4859e-03 2.09 2.1077e-02 1.87

Cartesian regression (must stay at roundoff):
n=11 rms 6.534e-14 n=17 rms 0.000e+00 n=25 rms 4.898e-13
```

**Scope.** `GI13` is called from three lines, all in `grad3DCurvLegacy`, which
is what the four-argument `grad3DCurv(k,X,Y,Z)` dispatches to. The
six-argument path uses different machinery and is unaffected.

**Known limitation.** At `k = 4` the patched operator still converges at order
2, because the interpolation stencils are second order. `grad2DCurv` behaves the
same way (`GI2` at `k = 4` also gives order 2.02), so this is pre-existing, not
introduced here. Recovering fourth order would need higher-order stencils in
both `GI2` and `GI13`.

This was originally reported by @manuelvalera

The solution could be this

```matlab
function I = GI13(M, m, n, o, type)
% Interpolates a logical-gradient component from the staggered face set it
% is computed on, to the face set it is needed on, for grad3DCurv.
%
% Parameters:
% M : the operator block to be interpolated (Ge, Gn or Gc)
% m, n, o : number of CELLS along xi, eta, zeta
% type : which shift to perform (see table below)
%
% Face layouts, flattened xi-fastest then eta then zeta:
% x-faces : (m+1) x n x o xi on nodes
% y-faces : m x (n+1) x o eta on nodes
% z-faces : m x n x (o+1) zeta on nodes
%
% Each shift is a tensor product of three 1-D operators, one per axis:
% type source -> target xi eta zeta
% 'Gn' y-face -> x-face c->n n->c I
% 'Gc' z-face -> x-face c->n I n->c
% 'Ge' x-face -> y-face n->c c->n I
% 'Gcy' z-face -> y-face I c->n n->c
% 'Gee' x-face -> z-face n->c I c->n
% 'Gnn' y-face -> z-face I n->c c->n

switch type
case 'Gn', Ax = P(m); Ay = Q(n); Az = speye(o);
case 'Gc', Ax = P(m); Ay = speye(n); Az = Q(o);
case 'Ge', Ax = Q(m); Ay = P(n); Az = speye(o);
case 'Gcy', Ax = speye(m); Ay = P(n); Az = Q(o);
case 'Gee', Ax = Q(m); Ay = speye(n); Az = P(o);
case 'Gnn', Ax = speye(m); Ay = Q(n); Az = P(o);
otherwise, error('GI13:BadType', 'unknown type "%s"', type);
end

I = kron(Az, kron(Ay, Ax)) * M;
end

function A = Q(N)
% node -> center, N+1 values on nodes -> N values on centers. Midpoint average.
A = spdiags(0.5*ones(N, 2), [0 1], N, N+1);
end

function A = P(N)
% center -> node, N values on centers -> N+1 values on nodes.
% Midpoint average in the interior, linear extrapolation at the two ends.
A = spdiags(0.5*ones(N+1, 2), [-1 0], N+1, N);
A(1, 1) = 1.5; A(1, 2) = -0.5; % extrapolate to node 0
A(N+1, N) = 1.5; A(N+1, N-1) = -0.5; % extrapolate to node N
end
```