diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py
index cd25fc26..90f912bc 100644
--- a/src/underworld3/swarm.py
+++ b/src/underworld3/swarm.py
@@ -2558,6 +2558,22 @@ def _update_proxy_variables(self):
                 # n, d, b = kd_swarm.find_closest_point(self._meshLevelSetVars[0].coords)
                 d, n = kd_swarm.query(self._meshLevelSetVars[0].coords, k=1, sqr_dists=False)
 
+                # Which (particle, node) pairs are valid:
+                # - node is at same distance as the nearest node
+                # - node is within radius_s
+                is_nearest = np.isclose(n_distance, n_distance[:, [0]])
+                within_radius = n_distance < self.radius_s
+                valid = is_nearest & within_radius
+
+                # IDW weights (only for valid pairs; others zeroed)
+                weights = 1.0 / (n_distance + 1e-16)
+                weights[~valid] = 0.0
+
+                # Total weight per node — material independent, compute once
+                n_mesh_nodes = self._meshLevelSetVars[0].data.shape[0]
+                w = np.zeros(n_mesh_nodes)
+                np.add.at(w, n_indices, weights)
+
             for ii in range(self.indices):
                 meshVar = self._meshLevelSetVars[ii]
 
@@ -2571,23 +2587,15 @@ def _update_proxy_variables(self):
                 final_values = np.array(meshVar.data[:, 0], copy=True)
 
                 if not starved:
+                    # Material presence mask: (n_particles, 1) for broadcasting
+                    mat_mask = (self.data.flatten() == ii).astype(weights.dtype)[:, None]
+
+                    # Weighted sum per node
                     node_values = np.zeros(final_values.shape[0])
-                    w = np.zeros(final_values.shape[0])
-
-                    for i in range(self.swarm.local_size):
-                        tem = np.isclose(n_distance[i, :], n_distance[i, 0])
-                        dist = n_distance[i, tem]
-                        indices = n_indices[i, tem]
-                        tem = dist < self.radius_s
-                        dist = dist[tem]
-                        indices = indices[tem]
-                        for j, ind in enumerate(indices):
-                            node_values[ind] += (
-                                np.isclose(self.data[i], ii) / (1.0e-16 + dist[j])
-                            )[0]
-                            w[ind] += 1.0 / (1.0e-16 + dist[j])
-
-                    node_values[np.where(w > 0.0)[0]] /= w[np.where(w > 0.0)[0]]
+                    np.add.at(node_values, n_indices, weights * mat_mask)
+
+                    # Normalize
+                    node_values[w > 0] /= w[w > 0]
                     final_values = node_values
 
                     # if there is no material found,
@@ -2602,49 +2610,67 @@ def _update_proxy_variables(self):
                 # current values, i.e. the proxy is left unchanged there)
                 meshVar.data[:, 0] = final_values
         elif self.update_type == 1:
-            # NOTE: this branch performs data-dependent MeshVariable writes
-            # outside any access context, which is not parallel-safe
-            # independently of the starved-rank issue (pre-existing).
-            # The guard here only prevents the empty-rank KDTree crash.
-            if starved:
-                return
-            kd = uw.kdtree.KDTree(self.swarm._particle_coordinates.data)
-            n_distance, n_indices = kd.query(
-                self._meshLevelSetVars[0].coords, k=self.nnn, sqr_dists=False
-            )
+            if not starved:
+                kd = uw.kdtree.KDTree(self.swarm._particle_coordinates.data)
+                n_distance, n_indices = kd.query(
+                    self._meshLevelSetVars[0].coords, k=self.nnn, sqr_dists=False
+                )
+
+                # IDW weights and validity mask for all (node, particle) pairs
+                valid = n_distance < self.radius_s
+                a = 1.0 / (n_distance + 1e-16)
+                a[~valid] = 0.0
+
+                # Total weight per node (material-independent)
+                w = a.sum(axis=1)
+
+                # Handle boundary nodes: restrict to nnn_bc particles
+                if hasattr(self, 'ind_bc') and self.ind_bc is not None:
+                    bc_idx = np.array(list(self.ind_bc))
+                    bc_idx = bc_idx[bc_idx < a.shape[0]]
+                    if len(bc_idx) > 0:
+                        valid_bc = n_distance[bc_idx, :self.nnn_bc] < self.radius_s
+                        a_bc = 1.0 / (n_distance[bc_idx, :self.nnn_bc] + 1e-16)
+                        a_bc[~valid_bc] = 0.0
+                        w_bc = a_bc.sum(axis=1)
+
+                        a[bc_idx] = 0.0
+                        a[bc_idx, :self.nnn_bc] = a_bc
+                        w[bc_idx] = w_bc
 
             for ii in range(self.indices):
                 meshVar = self._meshLevelSetVars[ii]
-                node_values = np.zeros((meshVar.data.shape[0],))
-                w = np.zeros((meshVar.data.shape[0],))
-                for i in range(meshVar.data.shape[0]):
-                    if i not in self.ind_bc:
-                        ind = np.where(n_distance[i, :] < self.radius_s)
-                        a = 1.0 / (n_distance[i, ind] + 1.0e-16)
-                        w[i] = np.sum(a)
-                        b = np.isclose(self.data[n_indices[i, ind]], ii)
-                        node_values[i] = np.sum(np.dot(a, b))
-                        if ind[0].size == 0:
-                            w[i] = 0
-                    else:
-                        ind = np.where(n_distance[i, : self.nnn_bc] < self.radius_s)
-                        a = 1.0 / (n_distance[i, : self.nnn_bc][ind] + 1.0e-16)
-                        w[i] = np.sum(a)
-                        b = np.isclose(self.data[n_indices[i, : self.nnn_bc][ind]], ii)
-                        node_values[i] = np.sum(np.dot(a, b))
-                        if ind[0].size == 0:
-                            w[i] = 0
-
-                node_values[np.where(w > 0.0)[0]] /= w[np.where(w > 0.0)[0]]
-                meshVar.data[:, 0] = node_values[...]
-
-                # if there is no material found,
-                # impose a near-neighbour hunt for a valid material and set that one
-                ind_w0 = np.where(w == 0.0)[0]
-                if len(ind_w0) > 0:
-                    ind_ = np.where(self.data[n_indices[ind_w0]] == ii)[0]
-                    if len(ind_) > 0:
-                        meshVar.data[ind_w0[ind_]] = 1.0
+
+                # MeshVariable reads/writes perform collective ghost
+                # synchronisation, so every rank must execute exactly the
+                # same read-then-write sequence per level set (same pattern
+                # as update_type=0 above). Starved ranks read their current
+                # proxy values and write them back unchanged; populated ranks
+                # compute and write new values.
+                final_values = np.array(meshVar.data[:, 0], copy=True)
+
+                if not starved:
+                    # Material presence at each (node, particle) pair
+                    # self.data has shape (n_particles, 1); flatten to (n_particles,)
+                    # so fancy indexing yields (n_nodes, nnn) — matching `a`.
+                    mat_present = (self.data.flatten()[n_indices] == ii).astype(a.dtype)
+
+                    # Weighted sum per node, then normalize
+                    node_values = (a * mat_present).sum(axis=1)
+                    node_values[w > 0] /= w[w > 0]
+                    final_values = node_values
+
+                    # if there is no material found,
+                    # impose a near-neighbour hunt for a valid material and set that one
+                    ind_w0 = np.where(w == 0.0)[0]
+                    if len(ind_w0) > 0:
+                        ind_ = np.where(self.data[n_indices[ind_w0]] == ii)[0]
+                        if len(ind_) > 0:
+                            final_values[ind_w0[ind_]] = 1.0
+
+                # single symmetric write (starved ranks write back their
+                # current values, i.e. the proxy is left unchanged there)
+                meshVar.data[:, 0] = final_values
         return
 
 
diff --git a/tests/benchmark_index_swarm_vectorized.py b/tests/benchmark_index_swarm_vectorized.py
new file mode 100644
index 00000000..33ba1901
--- /dev/null
+++ b/tests/benchmark_index_swarm_vectorized.py
@@ -0,0 +1,95 @@
+"""
+Benchmark the vectorized IndexSwarmVariable._update_proxy_variables
+against the original loop-based reference.
+
+Uses a larger mesh to make the speedup measurable.
+Run with: python tests/benchmark_index_swarm_vectorized.py
+"""
+import numpy as np
+import underworld3 as uw
+import time
+
+from test_0115_index_swarm_vectorized import (
+    _reference_update_type0,
+    _reference_update_type1,
+)
+
+
+def benchmark(update_type=1, mesh_res=(40, 40), fill_param=5, nnn=10, radius=0.1):
+    uw.reset_default_model()
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=mesh_res,
+        minCoords=(0.0, 0.0),
+        maxCoords=(10.0, 10.0),
+    )
+    swarm = uw.swarm.Swarm(mesh)
+    material = uw.swarm.IndexSwarmVariable(
+        "M_bench", swarm, indices=4, proxy_degree=1, proxy_continuous=True,
+        update_type=update_type, npoints=nnn, radius=radius,
+    )
+    swarm.populate(fill_param=fill_param)
+    with swarm.access():
+        # 4-layer material assignment
+        y = swarm.data[:, 1]
+        material.data[:, 0] = np.where(
+            y < 2.5, 0,
+            np.where(y < 5.0, 1,
+                     np.where(y < 7.5, 2, 3))
+        ).astype(int)
+
+    print(f"\n=== update_type={update_type}, mesh={mesh_res}, fill={fill_param}, "
+          f"nnn={nnn}, radius={radius} ===")
+    print(f"  Particles (local): {swarm.local_size}")
+    print(f"  Mesh nodes: {material._meshLevelSetVars[0].data.shape[0]}")
+    print(f"  Materials: {material.indices}")
+
+    # Benchmark vectorized (in-class)
+    material._proxy_stale = True
+    t0 = time.perf_counter()
+    material._update_proxy_if_stale()
+    t_vec = time.perf_counter() - t0
+    print(f"  Vectorized: {t_vec*1000:.2f} ms")
+
+    # Read vectorized result
+    vec = np.zeros((material._meshLevelSetVars[0].data.shape[0], material.indices))
+    for ii in range(material.indices):
+        vec[:, ii] = material._meshLevelSetVars[ii].data[:, 0]
+
+    # Benchmark reference (original loop)
+    t0 = time.perf_counter()
+    if update_type == 0:
+        ref = _reference_update_type0(swarm, material, nnn, radius)
+    else:
+        ref = _reference_update_type1(swarm, material, nnn, radius)
+    t_ref = time.perf_counter() - t0
+    print(f"  Reference:  {t_ref*1000:.2f} ms")
+
+    speedup = t_ref / t_vec if t_vec > 0 else float('inf')
+    print(f"  Speedup:    {speedup:.1f}x")
+
+    # Verify correctness
+    max_diff = np.max(np.abs(ref - vec))
+    print(f"  Max abs diff: {max_diff:.2e}")
+    assert np.allclose(ref, vec, atol=1e-10, rtol=1e-7), (
+        f"Vectorized output does not match reference! Max diff: {max_diff}"
+    )
+    print(f"  ✓ Output matches reference")
+    return speedup
+
+
+if __name__ == "__main__":
+    print("Benchmarking vectorized IndexSwarmVariable._update_proxy_variables")
+    print("=" * 70)
+
+    speedups = []
+    for res in [(10, 10), (20, 20), (40, 40), (60, 60)]:
+        for ut in [0, 1]:
+            s = benchmark(update_type=ut, mesh_res=res, nnn=10, radius=0.5)
+            speedups.append((ut, res, s))
+
+    print("\n" + "=" * 70)
+    print("Summary:")
+    print(f"{'update_type':<15} {'mesh':<15} {'speedup':<10}")
+    print("-" * 40)
+    for ut, res, s in speedups:
+        print(f"{ut:<15} {str(res):<15} {s:.1f}x")
diff --git a/tests/parallel/test_0767_index_swarm_advection_mpi.py b/tests/parallel/test_0767_index_swarm_advection_mpi.py
new file mode 100644
index 00000000..541208b8
--- /dev/null
+++ b/tests/parallel/test_0767_index_swarm_advection_mpi.py
@@ -0,0 +1,276 @@
+"""
+Parallel regression tests for IndexSwarmVariable proxy refresh after advection.
+
+Run under MPI, e.g.::
+
+    mpirun -np 4 python -m pytest tests/parallel/test_0767_index_swarm_advection_mpi.py
+    mpirun -np 2 python -m pytest tests/parallel/test_0767_index_swarm_advection_mpi.py
+
+Covers:
+
+- SWARM-07 / BF-06 follow-up: ``IndexSwarmVariable._update_proxy_variables()``
+  with ``update_type=1`` previously returned early on starved ranks (<= 1 local
+  particle) *without* participating in the collective MeshVariable read/write
+  sequence, causing an MPI deadlock when the populated ranks' ghost sync waited
+  for all ranks. The fix restructures the ``update_type=1`` branch to follow
+  the same collective pattern as ``update_type=0`` (and the base
+  ``SwarmVariable._rbf_to_meshVar``): starved ranks still read and write their
+  level-set proxy variables — they just leave the values unchanged.
+
+- Advection across rank boundaries followed by a projection solve that reads
+  the refreshed proxy, exercising the full advection → migrate → proxy-refresh
+  → solve path under MPI.
+"""
+
+import pytest
+import numpy as np
+import underworld3 as uw
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(120)]
+
+
+# ==============================================================================
+# Helpers
+# ==============================================================================
+
+SENTINEL = -123.0
+
+
+def _box(cell=0.25):
+    """Simple unit-square mesh for all tests."""
+    return uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cell,
+    )
+
+
+def _sentinel_level_sets(mat):
+    """Write SENTINEL into every level-set proxy variable (all ranks)."""
+    for lv in mat._meshLevelSetVars:
+        lv.data[:, 0] = SENTINEL
+
+
+def _check_starved_kept_sentinel(mat, starved):
+    """Assert starved-rank proxies still hold SENTINEL; populated ranks don't."""
+    for lv in mat._meshLevelSetVars:
+        vals = np.asarray(lv.data[:, 0])
+        if starved:
+            frac = float(np.mean(np.isclose(vals, SENTINEL)))
+            assert frac > 0.5, (
+                f"rank {uw.mpi.rank}: starved-rank proxy was overwritten "
+                f"(only {frac:.0%} of nodes kept the sentinel)"
+            )
+        else:
+            assert not np.allclose(vals, SENTINEL), (
+                f"rank {uw.mpi.rank}: populated-rank proxy was NOT refreshed "
+                f"(all nodes still hold the sentinel)"
+            )
+
+
+# ==============================================================================
+# Test 1: update_type=1 starved-rank proxy refresh (THE HANG REPRODUCER)
+# ==============================================================================
+
+@pytest.mark.mpi(min_size=4)
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_update_type1_starved_proxy_refresh():
+    """IndexSwarmVariable(update_type=1): starved ranks must participate in the
+    collective mesh-variable read/write, not return early. Before the fix this
+    deadlocked the populated ranks' ghost sync."""
+    mesh = _box()
+    swarm = uw.swarm.Swarm(mesh)
+    mat = uw.swarm.IndexSwarmVariable(
+        "M_starved_t1", swarm, indices=2, proxy_degree=1,
+        update_type=1, npoints=5, radius=0.5,
+    )
+
+    # Cluster particles near the origin — on np≥4 some ranks get ≤1 particle
+    xs = np.linspace(0.02, 0.08, 3)
+    cluster = np.column_stack([
+        np.repeat(xs, len(xs)),
+        np.tile(xs, len(xs)),
+    ])
+    swarm.add_particles_with_coordinates(cluster)
+
+    mat.data[...] = 0
+
+    starved = swarm.local_size <= 1
+    n_starved = uw.mpi.comm.allreduce(int(starved), op=uw.MPI.SUM)
+    assert n_starved >= 1, "test premise: at least one rank must be starved"
+    assert n_starved < uw.mpi.size, "test premise: one rank holds the cluster"
+
+    _sentinel_level_sets(mat)
+
+    # THIS is the code path that deadlocked before the fix:
+    # starved ranks return from _update_proxy_variables() without touching
+    # meshVar.data, while populated ranks do — causing a collective hang.
+    mat._proxy_stale = True
+    mat._update_proxy_if_stale()
+
+    _check_starved_kept_sentinel(mat, starved)
+
+    del swarm, mesh
+
+
+# ==============================================================================
+# Test 2: update_type=0 starved-rank proxy refresh (verification — must pass)
+# ==============================================================================
+
+@pytest.mark.mpi(min_size=4)
+@pytest.mark.level_1
+@pytest.mark.tier_a
+def test_update_type0_starved_proxy_refresh():
+    """IndexSwarmVariable(update_type=0): starved-rank collective read/write
+    was already correct; this test verifies nothing is broken."""
+    mesh = _box()
+    swarm = uw.swarm.Swarm(mesh)
+    mat = uw.swarm.IndexSwarmVariable(
+        "M_starved_t0", swarm, indices=2, proxy_degree=1,
+        update_type=0, npoints=5, radius=0.5,
+    )
+
+    xs = np.linspace(0.02, 0.08, 3)
+    cluster = np.column_stack([
+        np.repeat(xs, len(xs)),
+        np.tile(xs, len(xs)),
+    ])
+    swarm.add_particles_with_coordinates(cluster)
+    mat.data[...] = 0
+
+    starved = swarm.local_size <= 1
+    n_starved = uw.mpi.comm.allreduce(int(starved), op=uw.MPI.SUM)
+    assert n_starved >= 1
+    assert n_starved < uw.mpi.size
+
+    _sentinel_level_sets(mat)
+
+    mat._proxy_stale = True
+    mat._update_proxy_if_stale()
+
+    _check_starved_kept_sentinel(mat, starved)
+
+    del swarm, mesh
+
+
+# ==============================================================================
+# Test 3: advection across rank boundaries, then proxy solve
+# ==============================================================================
+
+@pytest.mark.mpi(min_size=2)
+@pytest.mark.level_1
+@pytest.mark.tier_a
+@pytest.mark.parametrize("update_type", [0, 1])
+def test_advection_then_proxy_solve(update_type):
+    """Advect a material blob across rank boundaries, then run a projection
+    solve that reads the proxy via .sym. Verifies the full advection→migrate
+    →proxy-refresh→solve path doesn't hang and produces correct values.
+
+    Uses a rightward-sheared velocity (v_x = y, v_y = 0) so every particle
+    at y > 0 moves right while staying inside [0, 1] — no domain-exit
+    deletion. Particles close to y=0 barely move; those near y=1 sweep
+    across partition boundaries.
+    """
+    import sympy
+    mesh = _box(cell=1.0 / 16.0)
+
+    # Shear flow: v_x = y, v_y = 0.  dt = 0.6 moves particles at y=1
+    # rightward by 0.6, staying in [0, 1] since the max x is 1.0 + 0.6
+    # but the box is [0, 1], so some may leave at the far right —
+    # but the shear means most particles remain in the box.
+    x, y = mesh.X
+    v_fn = sympy.Matrix([y, 0.0])
+
+    swarm = uw.swarm.Swarm(mesh)
+    mat = uw.swarm.IndexSwarmVariable(
+        f"M_shear_t{update_type}", swarm, indices=2, proxy_degree=1,
+        update_type=update_type, npoints=5, radius=0.3,
+    )
+    swarm.populate(fill_param=4)
+
+    # Tag left-half particles (x < 0.3) as material 1
+    mat.data[...] = 0
+    pc = swarm._particle_coordinates.data
+    mat.data[pc[:, 0] < 0.3, 0] = 1
+
+    def mat1_count():
+        return int((mat.data[:, 0] > 0.5).sum())
+
+    count_before = mat1_count()
+
+    # Advection: shear moves the left-half blob rightward; high-y particles
+    # move fast, low-y particles move slow. The blob should shift right
+    # across the partition boundary.
+    swarm.advection(v_fn, 0.5, order=2, step_limit=True)
+
+    # At least some particles must survive (the ones at low y stayed in box)
+    count_after = mat1_count()
+    assert count_after > 0, "all material-1 particles were deleted"
+
+    # Projection solve: reads mat.sym via the lazy path. Must NOT hang.
+    proj_var = uw.discretisation.MeshVariable(
+        "proj_0767", mesh, 1, degree=1,
+    )
+    proj = uw.systems.Projection(mesh, proj_var)
+    proj.uw_function = mat.createMask([0.0, 1.0])  # material 1 fraction
+    proj.petsc_options.delValue("ksp_monitor")
+    proj.solve()
+
+    # The solve must have produced some material-1 values. Not a tight
+    # correctness check — the test is primarily about "no hang" — but
+    # a sanity that the proxy was refreshed at all.
+    frac1 = np.asarray(proj_var.data[:, 0])
+    mean1 = float(frac1.mean())
+    assert mean1 > 0.01, (
+        f"Material-1 fraction implausibly low after advection "
+        f"(mean={mean1:.6f})"
+    )
+
+    del proj_var, proj, swarm, mesh
+
+
+# ==============================================================================
+# Test 4: advection that creates starved ranks, then proxy refresh
+# ==============================================================================
+
+@pytest.mark.mpi(min_size=4)
+@pytest.mark.level_1
+@pytest.mark.tier_a
+@pytest.mark.parametrize("update_type", [0, 1])
+def test_advection_creates_starved_ranks(update_type):
+    """Advection that pushes all particles into one side of the domain,
+    concentrating them on fewer ranks and leaving others starved. Then
+    triggers a proxy refresh via the solve path — must not hang."""
+    import sympy
+    mesh = _box(cell=1.0 / 8.0)
+
+    # Shear flow: v_x = y, v_y = 0.  High-y particles move far right;
+    # low-y particles stay near the left. With a dense uniform population,
+    # after migration the left-side ranks may become starved.
+    x, y = mesh.X
+    v_fn = sympy.Matrix([y, 0.0])
+
+    swarm = uw.swarm.Swarm(mesh)
+    mat = uw.swarm.IndexSwarmVariable(
+        f"M_starve_t{update_type}", swarm, indices=2, proxy_degree=1,
+        update_type=update_type, npoints=5, radius=0.3,
+    )
+    # Dense uniform population
+    swarm.populate(fill_param=6)
+
+    mat.data[...] = 0
+
+    # Advection: large dt pushes most particles to the right side.
+    # The migration concentrates them on right-side ranks.
+    swarm.advection(v_fn, 1.0, order=2, step_limit=True)
+
+    # Sentinel and force proxy refresh (exercises the collective path).
+    # This is what _sync_before_assembly does before a solve.
+    _sentinel_level_sets(mat)
+    mat._proxy_stale = True
+    mat._update_proxy_if_stale()
+
+    # No hang = success. On starved ranks the sentinel is preserved.
+    starved = swarm.local_size <= 1
+    _check_starved_kept_sentinel(mat, starved)
+
+    del swarm, mesh
diff --git a/tests/test_0115_index_swarm_vectorized.py b/tests/test_0115_index_swarm_vectorized.py
new file mode 100644
index 00000000..bf88c29a
--- /dev/null
+++ b/tests/test_0115_index_swarm_vectorized.py
@@ -0,0 +1,223 @@
+"""
+Verify the vectorized IndexSwarmVariable._update_proxy_variables against
+a reference implementation that mirrors the original (slow) per-particle loop.
+
+We construct a small mesh, populate a swarm, assign materials based on
+a horizontal interface, then compare the fractions produced by the
+in-class (vectorized) method against an independent reference computed
+by walking the particles one-by-one (the original algorithm).
+
+Both update_type=0 and update_type=1 are tested.
+"""
+import numpy as np
+import underworld3 as uw
+import pytest
+
+# All tests in this module are quick core tests
+pytestmark = pytest.mark.level_1
+
+
+def _reference_update_type0(swarm, material, nnn, radius_s):
+    """
+    Independent reimplementation of the ORIGINAL update_type=0 algorithm:
+    for each particle, find its nearest mesh node(s) (np.isclose),
+    accumulate 1/dist onto those nodes, normalize per node.
+    Returns an (n_nodes, n_indices) array of fractions.
+    """
+    meshVar0 = material._meshLevelSetVars[0]
+    kd_nodes = meshVar0._get_kdtree()
+    n_distance, n_indices = kd_nodes.query(
+        swarm._particle_coordinates.data, k=nnn, sqr_dists=False
+    )
+    kd_swarm = swarm._get_kdtree()
+    _, nearest_particle_per_node = kd_swarm.query(
+        meshVar0.coords, k=1, sqr_dists=False
+    )
+
+    n_nodes = meshVar0.data.shape[0]
+    n_indices_mat = material.indices
+    out = np.zeros((n_nodes, n_indices_mat))
+
+    for ii in range(n_indices_mat):
+        node_values = np.zeros(n_nodes)
+        w = np.zeros(n_nodes)
+        for i in range(swarm.local_size):
+            tem = np.isclose(n_distance[i, :], n_distance[i, 0])
+            dist = n_distance[i, tem]
+            indices = n_indices[i, tem]
+            tem = dist < radius_s
+            dist = dist[tem]
+            indices = indices[tem]
+            for j, ind in enumerate(indices):
+                node_values[ind] += (
+                    np.isclose(material.data[i], ii) / (1.0e-16 + dist[j])
+                )[0]
+                w[ind] += 1.0 / (1.0e-16 + dist[j])
+        node_values[w > 0] /= w[w > 0]
+        # fallback for uncovered nodes
+        ind_w0 = np.where(w == 0.0)[0]
+        if len(ind_w0) > 0:
+            ind_ = np.where(material.data[nearest_particle_per_node[ind_w0]] == ii)[0]
+            if len(ind_) > 0:
+                node_values[ind_w0[ind_]] = 1.0
+        out[:, ii] = node_values
+    return out
+
+
+def _reference_update_type1(swarm, material, nnn, radius_s, nnn_bc=None, ind_bc=None):
+    """
+    Independent reimplementation of the ORIGINAL update_type=1 algorithm:
+    for each mesh node, find its nnn nearest particles within radius_s,
+    compute IDW-weighted material fractions.
+    """
+    meshVar0 = material._meshLevelSetVars[0]
+    kd = uw.kdtree.KDTree(swarm._particle_coordinates.data)
+    n_distance, n_indices = kd.query(meshVar0.coords, k=nnn, sqr_dists=False)
+
+    n_nodes = meshVar0.data.shape[0]
+    n_indices_mat = material.indices
+    out = np.zeros((n_nodes, n_indices_mat))
+
+    bc_set = set(ind_bc) if ind_bc is not None else set()
+
+    for ii in range(n_indices_mat):
+        node_values = np.zeros(n_nodes)
+        w = np.zeros(n_nodes)
+        for i in range(n_nodes):
+            if i not in bc_set:
+                ind = np.where(n_distance[i, :] < radius_s)
+                a = 1.0 / (n_distance[i, ind] + 1.0e-16)
+                w[i] = np.sum(a)
+                b = np.isclose(material.data[n_indices[i, ind]], ii)
+                node_values[i] = np.sum(np.dot(a, b))
+                if ind[0].size == 0:
+                    w[i] = 0
+            else:
+                ind = np.where(n_distance[i, :nnn_bc] < radius_s)
+                a = 1.0 / (n_distance[i, :nnn_bc][ind] + 1.0e-16)
+                w[i] = np.sum(a)
+                b = np.isclose(material.data[n_indices[i, :nnn_bc][ind]], ii)
+                node_values[i] = np.sum(np.dot(a, b))
+                if ind[0].size == 0:
+                    w[i] = 0
+        node_values[w > 0] /= w[w > 0]
+        out[:, ii] = node_values
+    return out
+
+
+def _build_mesh_and_swarm(fill_param=5, nnn=5, radius=0.5, update_type=0,
+                           nnn_bc=None, ind_bc=None, interface_y=0.0,
+                           mesh_res=(4, 4)):
+    """Build a small mesh with a 2-material interface at y=interface_y."""
+    uw.reset_default_model()
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=mesh_res,
+        minCoords=(0.0, 0.0),
+        maxCoords=(1.0, 1.0),
+    )
+    swarm = uw.swarm.Swarm(mesh)
+    material = uw.swarm.IndexSwarmVariable(
+        "M_test", swarm, indices=2, proxy_degree=1, proxy_continuous=True,
+        update_type=update_type, npoints=nnn, radius=radius,
+        npoints_bc=nnn_bc if nnn_bc is not None else 2,
+        ind_bc=ind_bc,
+    )
+    swarm.populate(fill_param=fill_param)
+
+    # Assign materials based on y-coordinate relative to interface
+    with swarm.access():
+        material.data[:, 0] = np.where(
+            swarm.data[:, 1] <= interface_y, 0, 1
+        ).astype(int)
+    return mesh, swarm, material
+
+
+@pytest.mark.parametrize("update_type", [0, 1])
+@pytest.mark.parametrize("radius", [0.1, 0.3, 0.5, 1.0])
+@pytest.mark.parametrize("nnn", [3, 5, 10])
+def test_vectorized_matches_reference(update_type, radius, nnn):
+    """Compare the vectorized in-class update against the reference loop."""
+    mesh, swarm, material = _build_mesh_and_swarm(
+        fill_param=5, nnn=nnn, radius=radius, update_type=update_type,
+        interface_y=0.5, mesh_res=(4, 4),
+    )
+
+    # Compute reference using the slow loop
+    if update_type == 0:
+        ref = _reference_update_type0(swarm, material, nnn, radius)
+    else:
+        ref = _reference_update_type1(swarm, material, nnn, radius)
+
+    # Trigger the vectorized in-class update
+    material._proxy_stale = True
+    material._update_proxy_if_stale()
+
+    # Read the in-class fractions
+    vec = np.zeros_like(ref)
+    for ii in range(material.indices):
+        vec[:, ii] = material._meshLevelSetVars[ii].data[:, 0]
+
+    # Compare — allow small tolerance for floating-point summation order
+    assert np.allclose(ref, vec, atol=1e-10, rtol=1e-7), (
+        f"Mismatch for update_type={update_type}, radius={radius}, nnn={nnn}\n"
+        f"Max abs diff: {np.max(np.abs(ref - vec))}\n"
+        f"ref sum per material: {ref.sum(axis=0)}\n"
+        f"vec sum per material: {vec.sum(axis=0)}"
+    )
+
+
+def test_fractions_sum_to_one():
+    """For nodes with at least one particle nearby, fractions should sum to ~1."""
+    mesh, swarm, material = _build_mesh_and_swarm(
+        fill_param=5, nnn=10, radius=0.5, update_type=1,
+        interface_y=0.5, mesh_res=(6, 6),
+    )
+    material._proxy_stale = True
+    material._update_proxy_if_stale()
+
+    total = np.zeros(material._meshLevelSetVars[0].data.shape[0])
+    for ii in range(material.indices):
+        total += material._meshLevelSetVars[ii].data[:, 0]
+    # Where there's at least some weight, the sum should be ~1.0
+    # (allowing for the fallback behavior that may overshoot)
+    assert np.all(total >= -1e-10), "Fractions should be non-negative"
+    # Most interior nodes should sum to 1.0
+    assert np.allclose(total[total > 0.5], 1.0, atol=1e-6), (
+        f"Fractions don't sum to 1 where they should. "
+        f"Min/Max of total: {total.min():.6f} / {total.max():.6f}"
+    )
+
+
+def test_interface_location():
+    """The material interface should be near y=0.5."""
+    mesh, swarm, material = _build_mesh_and_swarm(
+        fill_param=5, nnn=10, radius=0.2, update_type=1,
+        interface_y=0.5, mesh_res=(8, 8),
+    )
+    material._proxy_stale = True
+    material._update_proxy_if_stale()
+
+    # Material 0 should be 1.0 below the interface, 0.0 above
+    coords = material._meshLevelSetVars[0].coords
+    frac0 = material._meshLevelSetVars[0].data[:, 0]
+
+    below = coords[:, 1] < 0.4
+    above = coords[:, 1] > 0.6
+    if below.any():
+        assert np.mean(frac0[below]) > 0.9, (
+            f"Material 0 fraction below interface too low: {np.mean(frac0[below]):.3f}"
+        )
+    if above.any():
+        assert np.mean(frac0[above]) < 0.1, (
+            f"Material 0 fraction above interface too high: {np.mean(frac0[above]):.3f}"
+        )
+
+
+if __name__ == "__main__":
+    # Run a quick smoke test
+    print("Running smoke test...")
+    test_vectorized_matches_reference(0, 0.5, 5)
+    test_vectorized_matches_reference(1, 0.5, 5)
+    test_fractions_sum_to_one()
+    test_interface_location()
+    print("All smoke tests passed.")