[BUG] - extract_region's _build_vertex_map calls _get_kdtree() on CoordinateSystem, which has no such method

[lmoresi]

Symptom

Mesh.extract_region(...) raises AttributeError on the first call,
because the post-construction step _build_vertex_map calls
self.X._get_kdtree() — and mesh.X is a CoordinateSystem object,
not a MeshVariable.

import underworld3 as uw
full = uw.meshing.AnnulusInternalBoundary(
    radiusOuter=1.0, radiusInner=0.5, radiusInternal=0.75, cellSize=0.2,
)
rock = full.extract_region("Inner")   # <-- AttributeError

File ".../discretisation_mesh.py", line 1270, in extract_region
    sub_mesh._build_vertex_map()
File ".../discretisation_mesh.py", line 1287, in _build_vertex_map
    tree = self.X._get_kdtree()
File ".../coordinates.py", line 2034, in __getattr__
    raise AttributeError(...)
AttributeError: 'CoordinateSystem' object has no attribute '_get_kdtree'

Root cause

PR #182 ("Consolidate and unify cached spatial indexing (KDTree)") unified
_get_kdtree() on MeshVariable / Swarm. The corresponding rewrite of
Mesh._build_vertex_map reads:

def _build_vertex_map(self):
    ...
    tree = self.X._get_kdtree()
    dists, indices = tree.query(self.parent.X.coords_nd, sqr_dists=False)
    ...

But mesh.X is self._CoordinateSystem (a CoordinateSystem instance), not a
MeshVariable. CoordinateSystem has no _get_kdtree, no coords_nd, and
its __getattr__ delegates to the underlying sympy Matrix, which doesn't
have these either — so the call fails.

extract_region is unreachable on development as a result. The MPI test
tests/parallel/test_0770_submesh_extract_mpi.py is skipped by default,
so this regression wasn't caught at PR merge time.

Suggested fix (small)

_build_vertex_map only needs a kdtree of vertex coordinates and the parent's
vertex coordinates as a query set — both available directly off the
CoordinateSystem's .coords array (with a .nondimensional() step to match
the original coords_nd semantics, or by reading mesh._coords directly).
Inline construction in the method body removes the dependence on a
non-existent CoordinateSystem method:

def _build_vertex_map(self):
    if hasattr(self, "_vertex_map") and self._vertex_map is not None:
        return self._vertex_map
    sub_coords    = self._coords    # underlying non-dimensional vertex coords
    parent_coords = self.parent._coords
    tree = uw.kdtree.KDTree(np.array(sub_coords))
    dists, indices = tree.query(np.array(parent_coords), sqr_dists=False)
    matched = dists < 1.0e-10
    self._vertex_map = (indices[matched], np.where(matched)[0])
    return self._vertex_map

A similar issue likely affects _build_dof_map (line 1482) — sub_var._get_kdtree()
does work because sub_var is a MeshVariable, but parent_var.coords_nd is fine.
That one is OK.

Why this surfaced now

Discovered while implementing the third submesh flavour (extract_surface,
codimension-1 boundary submesh via DMPlexCreateSubmesh) under
docs/examples/submesh_investigation/. The new prototype follows the
extract_region shape and tripped over this immediately.

Workaround

In investigation code: build the vertex map inline. The PR for the surface
flavour will do that and not depend on this method until the upstream fix
lands.

Underworld development team with AI support from Claude Code

[lmoresi]

Follow-up — there's actually a second bug that surfaces once the
`AttributeError` is worked around: `_build_dof_map` (line 1482) fails
the same way with `IndexError: too many indices for array`.

The root cause is in `uw.kdtree.KDTree.query`
(`src/underworld3/ckdtree.pyx:307`): when `k=1` it reshapes the
indices array to 1-D but leaves the distances array 2-D:

```python
i, d = self.find_closest_n_points(k, coords_contiguous)
if k == 1:
i = i.reshape(-1) # i becomes (n,)
# but d is still (n, 1) !
if sqr_dists:
return d, i
else:
distance_actual = numpy.sqrt(d) # also (n, 1)
...
return distance_actual, i
```

Reproducer:

```python
import numpy as np, underworld3 as uw
mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0,0), maxCoords=(1,1), cellSize=0.3)
v = uw.discretisation.MeshVariable('V', mesh, 1, degree=1)
tree = v._get_kdtree()
d, i = tree.query(np.asarray(v.coords_nd), sqr_dists=False)
print(d.shape, i.shape) # -> (30, 1), (30,)
m = d < 1e-10 # -> shape (30, 1), 2-D
print(i[m]) # -> IndexError: too many indices
```

Both `_build_vertex_map` and `_build_dof_map` write
`tree.query(...); matched = dists < 1e-10; sub_rows = indices[matched]`
— that pattern is only safe if both arrays have the same rank.

Suggested fix

In `KDTree.query`, also flatten `d` to match `i` when `k=1`:

```python
if k == 1:
i = i.reshape(-1)
d = d.reshape(-1)
```

This restores symmetry with the docstring ("Shape is (n_coords,) for
k=1"), which the index array currently satisfies but the distance array
doesn't.

Underworld development team with AI support from Claude Code