From a781acccd47cc66745c7046d27eaec42b84e66a9 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Sun, 1 Nov 2026 03:13:04 +0000
Subject: [PATCH 01/12] Contour fields evaluated where they are read

A space can mark itself evaluated, read through its view's value and
corner_value alone with no stored values to trace: contour's subdivision then
computes it at each lattice point, and its exact crossings find the field's
roots through value. A space marked continuous shares crossings between cells
as H1 does. Together they let contour cut a field computed on demand, an
expression of other fields say, without storing it.
---
 changelog.d/evaluated-fields.md                          | 8 ++++++++
 packages/tack-vis/src/tack/dataset/algorithms/contour.py | 9 ++++++---
 .../tack-vis/src/tack/dataset/algorithms/subdivide.py    | 5 ++++-
 3 files changed, 18 insertions(+), 4 deletions(-)
 create mode 100644 changelog.d/evaluated-fields.md

diff --git a/changelog.d/evaluated-fields.md b/changelog.d/evaluated-fields.md
new file mode 100644
index 000000000..4b768c318
--- /dev/null
+++ b/changelog.d/evaluated-fields.md
@@ -0,0 +1,8 @@
+### Added
+
+- Fields evaluated where they are read: a space whose `evaluated` attribute is
+  true is read through its view's `value` and `corner_value` alone, and
+  `contour`'s subdivision computes it at each lattice point; one whose
+  `continuous` attribute is true shares its crossings between cells, as `H1`
+  does, so its surface is watertight. An expression of continuous fields is
+  such a field.
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/contour.py b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
index fb53c1ee2..d848c8607 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/contour.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
@@ -458,8 +458,8 @@ def contour(data, field, value, merge_points=True, fields=None, dimension=None,
     is placed on an edge to a missing value.
 
     With ``merge_points``, triangles share their points: by global edge when
-    the field and the geometry are both continuous, so the surface is
-    watertight; otherwise only within each cell, so a DG field's surface
+    the field and the geometry are both continuous (an ``H1`` field, or a space
+    whose ``continuous`` attribute says so), so the surface is watertight; otherwise only within each cell, so a DG field's surface
     keeps its jumps and an ``L2`` geometry's its gaps. ``H1`` fields carry
     onto the surface's points, interpolated along the same edges; cell
     fields (``Constant``, values on cells) onto its triangles, from the cell
@@ -517,7 +517,10 @@ def contour(data, field, value, merge_points=True, fields=None, dimension=None,
     space = field.space
     if not space.interpolated or space.on not in ("cells", "points"):
         raise TypeError(f"contour needs a field with a basis on the cells, not {space!r}")
-    per_cell = 0 if isinstance(space, H1) and isinstance(data.geometry.space, H1) else 1
+    # A space that declares itself continuous (an expression of continuous fields
+    # evaluated where it is read, say) shares its crossings between cells as H1 does.
+    continuous = isinstance(space, H1) or getattr(space, "continuous", False)
+    per_cell = 0 if continuous and isinstance(data.geometry.space, H1) else 1
     groups = [g for g in data.launch_groups("cells", [field])
               if g.count and g.shape.CONTOUR_TRIANGLES]
     n_cells = data.num_cells
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py b/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
index 7d846f491..ce2106b9b 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
@@ -1026,7 +1026,10 @@ def _refined_family(data, flags, level, fields, elsewhere, cut, dimension, cell_
     linear = UnstructuredTopology(types, offsets, conn, num_points=points, id_dtype=idt)
     point_fields = {}
     for name, field in fields.items():
-        if not isinstance(field.space, (H1, L2)) or field.space.varies:
+        # A space marked evaluated has no stored values to trace: it is computed where it
+        # is read, here at each lattice point by its value.
+        evaluated = getattr(field.space, "evaluated", False)
+        if not (isinstance(field.space, (H1, L2)) or evaluated) or field.space.varies:
             continue
         out = tack.empty((points,), _element(field.values))
         if count:

From 377dca430bdb49537577bc8dc86ac1486269dc41 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Tue, 3 Nov 2026 20:15:16 +0000
Subject: [PATCH 02/12] One resolver for field arguments; filters carry given
 gradients and rebind evaluated fields

Text that names no field goes to a registered resolver, so field expressions
work wherever a field does. A space evaluated where it is read gives its own
default subdivision, is made again on a filter's output, and is written as its
values.
---
 changelog.d/field-resolver.md                 | 15 +++++++
 .../tack-vis/src/tack/dataset/__init__.py     |  3 +-
 .../src/tack/dataset/algorithms/_common.py    | 30 ++++++++++++--
 .../src/tack/dataset/algorithms/contour.py    |  3 +-
 .../src/tack/dataset/algorithms/subdivide.py  | 12 +++---
 packages/tack-vis/src/tack/dataset/carry.py   | 21 +++++++++-
 packages/tack-vis/src/tack/dataset/dataset.py | 15 +++++++
 .../tack-vis/src/tack/interop/vtk/__init__.py |  4 ++
 .../tack-vis/tests/test_field_resolver.py     | 40 +++++++++++++++++++
 .../tack-vis/tests/test_given_gradients.py    | 15 +++++++
 10 files changed, 146 insertions(+), 12 deletions(-)
 create mode 100644 changelog.d/field-resolver.md
 create mode 100644 packages/tack-vis/tests/test_field_resolver.py

diff --git a/changelog.d/field-resolver.md b/changelog.d/field-resolver.md
new file mode 100644
index 000000000..ad6572eeb
--- /dev/null
+++ b/changelog.d/field-resolver.md
@@ -0,0 +1,15 @@
+### Added
+
+- `register_field_resolver(resolve)`: every algorithm's field argument may be text
+  that names no field, which `resolve(data, text)` turns into a `Field`. Field
+  expressions register one, so `contour(data, "mag(U)", 1.0)` contours an
+  expression.
+- A space evaluated where it is read may give `default_subdivision`, the lattice
+  steps contour, clip and threshold cut its cells into by default.
+
+### Changed
+
+- Filters carry a field's given gradient with it, and make a field evaluated where
+  it is read again on their output (its space's `rebind`) instead of carrying
+  its placeholder values.
+- `dataset_to_vtk` writes such a field's values (its space's `stored`).
diff --git a/packages/tack-vis/src/tack/dataset/__init__.py b/packages/tack-vis/src/tack/dataset/__init__.py
index a71a6782c..6c81e746b 100644
--- a/packages/tack-vis/src/tack/dataset/__init__.py
+++ b/packages/tack-vis/src/tack/dataset/__init__.py
@@ -109,7 +109,7 @@
     voxelize,
     warp,
 )
-from tack.dataset.algorithms._common import edge_crossing
+from tack.dataset.algorithms._common import edge_crossing, register_field_resolver
 from tack.dataset.algorithms.gradients import world_gradient, world_hessian, world_jet
 from tack.dataset.amr import AMRDataSet
 from tack.dataset.carry import Interpolate, Pieces, Same, Take, carry
@@ -282,6 +282,7 @@
     "probe",
     "project",
     "rectilinear_grid",
+    "register_field_resolver",
     "remove_ghosts",
     "shapes",
     "shrink",
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_common.py b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
index 612c2b058..4db25cfdf 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -136,13 +136,37 @@ def _same_space(*fields):
     return space, n
 
 
+_RESOLVERS = []
+
+
+def register_field_resolver(resolve):
+    """Lets every algorithm's field argument be text that names no field: ``resolve(data,
+    text)`` gives a ``Field`` on ``data`` for it, or ``None`` when the text is not its kind
+    (the field expressions of ``tack_expr`` register one). The last registered is asked
+    first."""
+    if resolve not in _RESOLVERS:
+        _RESOLVERS.append(resolve)
+
+
+def _resolved(data, text):
+    if not isinstance(text, str) or text in data.fields:
+        return data.fields[text]
+    for resolve in reversed(_RESOLVERS):
+        field = resolve(data, text)
+        if field is not None:
+            return field
+    return data.fields[text]
+
+
 def _field(data, field):
-    """A field by name or object, checked against the dataset before launch."""
+    """A field by name, object or registered resolver (``register_field_resolver``),
+    checked against the dataset before launch."""
     from tack.dataset.dataset import Field, is_field_key
 
-    field = data.fields[field] if is_field_key(field) else field
+    field = _resolved(data, field) if is_field_key(field) else field
     if not isinstance(field, Field):
-        raise TypeError("a field argument is a tack.dataset.Field or its name")
+        raise TypeError("a field argument is a tack.dataset.Field, its name, or text a "
+                        "registered resolver reads")
     if field.space.topology is not data.topology:
         raise ValueError("a field argument belongs to another topology")
     return field
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/contour.py b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
index d848c8607..e5bd2ff31 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/contour.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
@@ -465,7 +465,8 @@ def contour(data, field, value, merge_points=True, fields=None, dimension=None,
     fields (``Constant``, values on cells) onto its triangles, from the cell
     each came from.
     """
-    value = _in_unit_of(_field(data, field), value, "contour's value")
+    field = _field(data, field)
+    value = _in_unit_of(field, value, "contour's value")
     from tack.dataset.ghosts import visible
 
     data, field = visible(data, field)
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py b/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
index ce2106b9b..d78528ea2 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
@@ -862,7 +862,8 @@ def _order(field):
     space = field.space
     if isinstance(space, (H1, L2)) and not space.varies:
         return space.order
-    return 1
+    # A space evaluated where it is read (an expression's) says how finely to sample it.
+    return getattr(space, "default_subdivision", 1)
 
 
 def default_level(data, *fields):
@@ -1373,10 +1374,11 @@ def subdivides(data, field):
     """Whether cutting ``field`` on ``data`` takes subdivision: a curved geometry or a
     field above order 1, on tensor cells (segments, quads, pixels, hexahedra,
     voxels)."""
-    nonlinear = _order(data.geometry) > 1 or (isinstance(field.space, (H1, L2))
-                                              and not field.space.varies
-                                              and (field.space.order > 1
-                                                   or getattr(field.space, "tensor", False)))
+    space = field.space
+    nonlinear = (_order(data.geometry) > 1
+                 or (isinstance(space, (H1, L2)) and not space.varies
+                     and (space.order > 1 or getattr(space, "tensor", False)))
+                 or getattr(space, "default_subdivision", 1) > 1)
     groups = [g for g in data.topology.groups() if g.count]
     dims = {g.shape.DIMENSION for g in groups}
     return nonlinear and len(dims) <= 1 and all(bases.has_bases(g.shape) for g in groups)
diff --git a/packages/tack-vis/src/tack/dataset/carry.py b/packages/tack-vis/src/tack/dataset/carry.py
index f1b2c5c84..4b700e8fe 100644
--- a/packages/tack-vis/src/tack/dataset/carry.py
+++ b/packages/tack-vis/src/tack/dataset/carry.py
@@ -499,9 +499,19 @@ def carry(data, out, *, points=None, cells=None, faces=None, faces_to_cells=None
         if geometry is None:
             raise ValueError("the geometry cannot be carried by these maps; pass it")
     carried = {}
+    rebound = []
     for name in selected(data, fields):
-        field = _carry_field(data.fields[name], data, out, points, cells, faces,
+        source = data.fields[name]
+        if getattr(source.space, "evaluated", False) and hasattr(source.space, "rebind"):
+            rebound.append((name, source))
+            continue
+        field = _carry_field(source, data, out, points, cells, faces,
                              faces_to_cells, edges, edges_to_cells)
+        if field is not None and getattr(source, "gradient", None) is not None:
+            gradient = _carry_field(source.gradient, data, out, points, cells, faces,
+                                    faces_to_cells, edges, edges_to_cells)
+            if gradient is not None:
+                field = Field(field.space, field.values, field.unit, gradient=gradient)
         if field is not None:
             carried[name] = field
     if sets is None:
@@ -518,5 +528,12 @@ def carry(data, out, *, points=None, cells=None, faces=None, faces_to_cells=None
     made = DataSet(out, geometry, fields=carried, sets=sets, **data._frame())
     carry_ids(data, made, points=points, cells=None if faces_to_cells is not None or
               edges_to_cells is not None else cells)
-    return carry_marks(data, made, points=points, cells=cells, faces_to_cells=faces_to_cells,
+    made = carry_marks(data, made, points=points, cells=cells, faces_to_cells=faces_to_cells,
                        edges_to_cells=edges_to_cells)
+    # A field evaluated where it is read has nothing to move: its space makes it again
+    # on the output, from the fields carried there (None where they were not).
+    for name, source in rebound:
+        field = source.space.rebind(source, made)
+        if field is not None:
+            made.fields[name] = field
+    return made
diff --git a/packages/tack-vis/src/tack/dataset/dataset.py b/packages/tack-vis/src/tack/dataset/dataset.py
index cf4bd65b9..858dd6f1f 100644
--- a/packages/tack-vis/src/tack/dataset/dataset.py
+++ b/packages/tack-vis/src/tack/dataset/dataset.py
@@ -733,6 +733,21 @@ def __repr__(self):
         return f"Fields({{{', '.join(f'{k!r}: {f.space!r}' for k, f in self.items())}}})"
 
 
+def stored(data):
+    """``data`` with each field whose space is evaluated where it is read (an
+    expression's) replaced by its values, as the space's ``stored`` gives them: for
+    writers, which take arrays. ``data`` itself when it has none."""
+    entries = data.fields.entries()
+    if not any(callable(getattr(f.space, "stored", None)) for _, f in entries):
+        return data
+    import copy
+
+    out = copy.copy(data)
+    out.fields = [(name, f.space.stored(f, data) if callable(getattr(f.space, "stored", None))
+                   else f) for name, f in entries]
+    return out
+
+
 class DataSet:
     """A topology, a geometry, named fields and named sets. The geometry is a
     ``Field`` in ``H1`` or ``L2`` of 3-vectors -- or positions per point (a host
diff --git a/packages/tack-vis/src/tack/interop/vtk/__init__.py b/packages/tack-vis/src/tack/interop/vtk/__init__.py
index 23cfd5f92..b48605f37 100644
--- a/packages/tack-vis/src/tack/interop/vtk/__init__.py
+++ b/packages/tack-vis/src/tack/interop/vtk/__init__.py
@@ -2197,8 +2197,12 @@ def dataset_to_vtk(data, copy=True, higher_order=False, cell_grid=False):
     ``Constant`` fields and values on cells ``constant``; fields of other spaces
     are left out with a warning. The cells come type by type, in Tack's order
     within each. It needs ``vtkDGCell.SetCellSpec``, newer than VTK 9.7.1."""
+    from tack.dataset import DataSet
+    from tack.dataset.dataset import stored
     from tack.interop.vtk import _lagrange
 
+    if isinstance(data, DataSet):
+        data = stored(data)
     if cell_grid:
         from tack.interop.vtk import _cellgrid
 
diff --git a/packages/tack-vis/tests/test_field_resolver.py b/packages/tack-vis/tests/test_field_resolver.py
new file mode 100644
index 000000000..ef04dc5f0
--- /dev/null
+++ b/packages/tack-vis/tests/test_field_resolver.py
@@ -0,0 +1,40 @@
+"""Field arguments as text a registered resolver reads."""
+
+import numpy as np
+import pytest
+
+import tack
+import tack.dataset as td
+from tack.dataset.algorithms import _common
+
+
+@tack.func
+def _u(x):
+    return x[0] + 2.0 * x[1] - x[2]
+
+
+@pytest.fixture
+def twice():
+    def resolve(data, text):
+        if not text.startswith("twice "):
+            return None
+        u = data.fields[text[len("twice "):]]
+        values = tack.array(2.0 * u.values.numpy())
+        return td.Field(u.space, values, u.unit)
+
+    td.register_field_resolver(resolve)
+    yield resolve
+    _common._RESOLVERS.remove(resolve)
+
+
+def test_algorithms_take_text_a_resolver_reads(backend, twice):
+    data = td.uniform_grid((5, 4, 3), origin=(0.1, 0.2, 0.3), spacing=(0.3, 0.25, 0.4))
+    data.fields["u"] = td.interpolate(data, _u, td.H1(data))
+    doubled = td.Field(data.fields["u"].space, tack.array(2.0 * data.fields["u"].values.numpy()))
+    a = td.contour(data, "twice u", 1.5)
+    b = td.contour(data, doubled, 1.5)
+    assert a.num_points == b.num_points > 0
+    np.testing.assert_allclose(td.gradients(data, "twice u").values.numpy(),
+                               td.gradients(data, doubled).values.numpy())
+    with pytest.raises(KeyError):
+        td.contour(data, "missing", 1.0)
diff --git a/packages/tack-vis/tests/test_given_gradients.py b/packages/tack-vis/tests/test_given_gradients.py
index b01df215c..08b070a7d 100644
--- a/packages/tack-vis/tests/test_given_gradients.py
+++ b/packages/tack-vis/tests/test_given_gradients.py
@@ -112,3 +112,18 @@ def test_point_data_keeps_its_given_gradient(backend):
     given = np.tile(np.float32([4.0, 5.0, 6.0]), (data.num_cells, 1))
     data.fields["u"] = td.Field(td.Values(data, "points"), u.values, gradient=_per_cell(data, given))
     np.testing.assert_allclose(td.gradients(data, "u").values.numpy(), given)
+
+
+def test_a_given_gradient_goes_with_its_field_through_filters(backend):
+    data = _data()
+    given = np.random.default_rng(2).uniform(-1, 1, (data.num_cells, 3)).astype(np.float32)
+    u = td.interpolate(data, _u, td.H1(data))
+    data.fields["u"] = td.Field(u.space, u.values, gradient=_per_cell(data, given))
+    centers = td.cell_centers(data).values.numpy()
+    keep = centers[:, 0] < 0.7
+    data.fields["x"] = td.Field(td.Constant(data), tack.array(centers[:, 0].astype(np.float32)))
+    out = td.threshold(data, "x", -1.0, 0.7)
+    assert out.num_cells == int(keep.sum())
+    assert out.fields["u"].gradient is not None
+    np.testing.assert_allclose(td.gradients(out, "u").values.numpy(), given[keep])
+

From 74076e3b25ed2432dba0b371e87c28e819f6e657 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Tue, 3 Nov 2026 21:23:53 +0000
Subject: [PATCH 03/12] Reference entry for register_field_resolver

---
 docs/vis-guide/reference/data-model.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/vis-guide/reference/data-model.md b/docs/vis-guide/reference/data-model.md
index 4a29b6c13..e7874603d 100644
--- a/docs/vis-guide/reference/data-model.md
+++ b/docs/vis-guide/reference/data-model.md
@@ -205,6 +205,8 @@ do, so pieces sharing an edge compute it alike.
 
 ::: tack.dataset.edge_crossing
 
+::: tack.dataset.register_field_resolver
+
 ::: tack.dataset.keep_marks
 
 ::: tack.dataset.drop_marks

From 49cf7c15ef0d9aac3e193c38ed46b7b39d4c0bcc Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Tue, 3 Nov 2026 19:54:25 +0000
Subject: [PATCH 04/12] Fields given their gradient; world derivatives from the
 field's view

A solver's derivatives are its own: a finite-volume reconstruction, or a
lattice Boltzmann solver's from its populations. Field(..., gradient=) takes
that gradient as a field on the same topology (gradients' layout), and
world_gradient and world_jet ask the field's view, which reads it, or a space's
own derivatives, or else differentiates the basis as before. gradients splits
a vector's given gradient into its rows, so flow_quantities and every kernel on
world_gradient use the solver's derivatives.
---
 changelog.d/given-gradients.md                |  14 ++
 .../src/tack/dataset/algorithms/_common.py    |   4 +-
 .../src/tack/dataset/algorithms/gradients.py  | 148 ++++++++++++------
 packages/tack-vis/src/tack/dataset/dataset.py |  49 +++++-
 .../tack-vis/tests/test_given_gradients.py    | 114 ++++++++++++++
 5 files changed, 278 insertions(+), 51 deletions(-)
 create mode 100644 changelog.d/given-gradients.md
 create mode 100644 packages/tack-vis/tests/test_given_gradients.py

diff --git a/changelog.d/given-gradients.md b/changelog.d/given-gradients.md
new file mode 100644
index 000000000..9e3ba6b8c
--- /dev/null
+++ b/changelog.d/given-gradients.md
@@ -0,0 +1,14 @@
+### Added
+
+- `Field(space, values, unit, gradient=)`: a field's world gradient where
+  something other than its basis gives it, such as a solver's own
+  reconstruction, as a field on the same topology (3 values each for a scalar,
+  3k row by row for a vector of k). `world_gradient`, `world_jet` and
+  `gradients` read it, so `flow_quantities` and every kernel built on them use
+  the solver's derivatives; the Hessian it does not determine is NaN.
+
+### Changed
+
+- `world_gradient` and `world_jet` ask the field's view, whose default
+  (on every field view whose basis gives the derivatives) differentiates its basis as
+  before. A space that computes its values can give its own.
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_common.py b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
index 5af7990a6..612c2b058 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -157,9 +157,9 @@ def _with_basis(data, field):
 
     space = field.space
     if isinstance(space, Values) and space.on == "points":
-        return Field(H1(data), field.values, field.unit)
+        return Field(H1(data), field.values, field.unit, gradient=field.gradient)
     if isinstance(space, Values) and space.on == "cells":
-        return Field(Constant(data), field.values, field.unit)
+        return Field(Constant(data), field.values, field.unit, gradient=field.gradient)
     return field
 
 
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/gradients.py b/packages/tack-vis/src/tack/dataset/algorithms/gradients.py
index 74efcd1b9..523ba137b 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/gradients.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/gradients.py
@@ -17,20 +17,33 @@
 
 @tack.func
 def world_gradient(cells, u, c, pc):
-    """`u`'s gradient in the world at parametric `pc` of cell `c`: its
-    parametric gradient through the geometry's Jacobian. For kernels that need a
-    field's gradient where `gradients` does not put it."""
-    # d(u)/d(x) from d(u)/d(pc), through the geometry's Jacobian J. A solid's is
-    # J^-T d(u)/d(pc). A 2D cell's lies in its tangent plane, spanned by J's
-    # columns a and b: g = alpha a + beta b with g . a and g . b the parametric
-    # derivatives, VTK's and Viskores' in-plane gradient (no area: NaN, as a
-    # solid of no volume). Lines and vertices: zero.
-    g = tack.vec3(0.0, 0.0, 0.0)
-    if cells.DIMENSION == 3:
-        g = cells.geometry_jacobian(c, pc).inverse().transpose() @ u.parametric_gradient(c, pc)
-    elif cells.DIMENSION == 2:
-        g = _in_plane(cells.geometry_jacobian(c, pc), u.parametric_gradient(c, pc))
-    return g
+    """`u`'s gradient in the world at parametric `pc` of cell `c`, as its view gives it:
+    from its basis by default (`_BasisGradient`), the gradient it was given
+    (`Field(..., gradient=)`), or a space's own. For kernels that need a field's
+    gradient where `gradients` does not put it."""
+    return u.world_gradient(cells, c, pc)
+
+
+@tack.func
+def world_hessian(cells, u, c, pc):
+    """`u`'s second derivatives in the world at parametric `pc` of solid cell `c`, a
+    symmetric 3x3 matrix, as its view gives them (`world_jet`). From the basis: by the
+    chain rule through the geometry, the parametric Hessian less the geometry's
+    curvature along the world gradient, taken to world axes on both sides by the
+    inverse Jacobian; analytic on every shape at orders 1 and 2, and at any order on
+    tensor shapes, simplices and wedges; zero on surfaces, lines and vertices. NaN for
+    a field given only its gradient."""
+    _, _, h = world_jet(cells, u, c, pc)
+    return h
+
+
+@tack.func
+def world_jet(cells, u, c, pc):
+    """Scalar `u`'s value, world gradient and world Hessian at parametric `pc` of cell
+    `c`, as its view gives them, from one pass over the cell's values where it can:
+    what a kernel needing several of them reads once. The gradient is
+    `world_gradient`'s, and the Hessian `world_hessian`'s."""
+    return u.world_jet(cells, c, pc)
 
 
 @tack.func
@@ -46,36 +59,57 @@ def _in_plane(m, d):
     return ((bb * d[0] - ab * d[1]) * a + (aa * d[1] - ab * d[0]) * b) / det
 
 
-
-@tack.func
-def world_hessian(cells, u, c, pc):
-    """`u`'s second derivatives in the world at parametric `pc` of solid cell `c`, a
-    symmetric 3x3 matrix: by the chain rule through the geometry, the parametric
-    Hessian less the geometry's curvature along the world gradient, taken to world
-    axes on both sides by the inverse Jacobian. Analytic from the basis on every
-    shape at orders 1 and 2, and at any order on tensor shapes, simplices and
-    wedges. Zero on surfaces, lines and vertices."""
-    _, _, h = world_jet(cells, u, c, pc)
-    return h
-
-
-@tack.func
-def world_jet(cells, u, c, pc):
-    """Scalar `u`'s value, world gradient and world Hessian at parametric `pc` of cell
-    `c`, from the field view's `parametric_jet`, one pass over the cell's values: what a
-    kernel needing several of them reads once. The gradient is `world_gradient`'s on
-    every cell, in the cell's plane on surfaces, and the Hessian `world_hessian`'s, zero
-    on surfaces, lines and vertices."""
-    v, pg, ph = u.parametric_jet(c, pc)
-    g = tack.vec3(0.0, 0.0, 0.0)
-    h = tack.mat3([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0])
-    if cells.DIMENSION == 3:
-        inverse = cells.geometry_jacobian(c, pc).inverse()
-        g = inverse.transpose() @ pg
-        h = inverse.transpose() @ (ph - cells.geometry_curvature(c, pc, g)) @ inverse
-    elif cells.DIMENSION == 2:
-        g = _in_plane(cells.geometry_jacobian(c, pc), pg)
-    return v, g, h
+class _BasisGradient:
+    """A field view's world gradient from its basis: on every field view with a
+    ``parametric_gradient`` and no other source (``_derivative_mixins``)."""
+
+    @tack.func
+    def world_gradient(self, cells, c, pc):
+        # d(u)/d(x) from d(u)/d(pc), through the geometry's Jacobian J. A solid's is
+        # J^-T d(u)/d(pc). A 2D cell's lies in its tangent plane, spanned by J's
+        # columns a and b: g = alpha a + beta b with g . a and g . b the parametric
+        # derivatives, VTK's and Viskores' in-plane gradient (no area: NaN, as a
+        # solid of no volume). Lines and vertices: zero.
+        g = tack.vec3(0.0, 0.0, 0.0)
+        if cells.DIMENSION == 3:
+            g = cells.geometry_jacobian(c, pc).inverse().transpose() @ self.parametric_gradient(c, pc)
+        elif cells.DIMENSION == 2:
+            g = _in_plane(cells.geometry_jacobian(c, pc), self.parametric_gradient(c, pc))
+        return g
+
+
+class _BasisJet:
+    """A field view's value, world gradient and world Hessian from its basis, in one
+    pass: on every field view with a ``parametric_jet``."""
+
+    @tack.func
+    def world_jet(self, cells, c, pc):
+        v, pg, ph = self.parametric_jet(c, pc)
+        g = tack.vec3(0.0, 0.0, 0.0)
+        h = tack.mat3([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0])
+        if cells.DIMENSION == 3:
+            inverse = cells.geometry_jacobian(c, pc).inverse()
+            g = inverse.transpose() @ pg
+            h = inverse.transpose() @ (ph - cells.geometry_curvature(c, pc, g)) @ inverse
+        elif cells.DIMENSION == 2:
+            g = _in_plane(cells.geometry_jacobian(c, pc), pg)
+        return v, g, h
+
+
+class _GivenGradient:
+    """A field given its gradient (``Field(..., gradient=)``): ``world_gradient`` is that
+    field's value at ``pc`` of cell ``c``, and ``world_jet`` pairs it with the field's
+    value and a Hessian of NaN, which a gradient alone does not determine."""
+
+    @tack.func
+    def world_gradient(self, cells, c, pc):
+        return self.given_gradient.value(c, pc)
+
+    @tack.func
+    def world_jet(self, cells, c, pc):
+        n = self.no_hessian
+        return (self.value(c, pc), self.given_gradient.value(c, pc),
+                tack.mat3([n, n, n], [n, n, n], [n, n, n]))
 
 
 @tack.kernel
@@ -275,6 +309,24 @@ def _polyhedral_scalar_gradients(data, field, at):
     return out
 
 
+@tack.kernel
+def _row(values, k, out):
+    for i in range(out.shape[0]):
+        v = values[i]
+        out[i] = tack.vec3(v[3 * k], v[3 * k + 1], v[3 * k + 2])
+
+
+def _gradient_row(gradient, k):
+    """Row ``k`` of a vector field's given gradient: its component ``k``'s."""
+    if gradient is None:
+        return None
+    values = materialize(gradient.values)
+    out = tack.empty((size_of(values),), tack.vec(3, arrays.dtype_of(values)))
+    if out.shape[0]:
+        _row(values, k, out)
+    return Field(gradient.space, out, gradient.unit)
+
+
 @tack.kernel
 def _component(values, k, out):
     for i in range(out.shape[0]):
@@ -352,6 +404,9 @@ def gradients(data, field, at="cells"):
     Viskores adds nothing but still counts the cell). ``flow_quantities`` derives divergence, vorticity and
     the Q-criterion from a vector field's gradient.
 
+    A field given its gradient (``Field(..., gradient=)``) gives that instead of its
+    basis', evaluated at the same points and averaged the same way.
+
     On a polyhedral topology there is no reference element. A cell's gradient
     of point data is then Green-Gauss over its boundary: the faces walked
     outward and fanned from their points' means (the triangles ``cell_geometry``
@@ -362,6 +417,9 @@ def gradients(data, field, at="cells"):
     get zero. At the points, the cells' gradients are averaged over the cells
     around each point."""
     field = _with_basis(data, _field(data, field))
+    if not getattr(data.topology, "reference_cells", True) and field.gradient is not None:
+        raise NotImplementedError("gradients of a field given its gradient on a mesh without "
+                                  "reference elements")
     if (not getattr(data.topology, "reference_cells", True)
             and not isinstance(field.space, (H1, Constant))):
         raise NotImplementedError(f"gradients of a field on {field.space!r} on a mesh "
@@ -391,7 +449,7 @@ def gradients(data, field, at="cells"):
         part = tack.empty((size_of(field.values),), arrays.dtype_of(field.values))
         if part.shape[0]:
             _component(values, k, part)
-        rows = _scalar_gradients(data, Field(field.space, part), at)
+        rows = _scalar_gradients(data, Field(field.space, part, gradient=_gradient_row(field.gradient, k)), at)
         if n:
             _place_rows(rows, k, out)
     return Field(space, out, unit)
diff --git a/packages/tack-vis/src/tack/dataset/dataset.py b/packages/tack-vis/src/tack/dataset/dataset.py
index e005fb4b9..cf4bd65b9 100644
--- a/packages/tack-vis/src/tack/dataset/dataset.py
+++ b/packages/tack-vis/src/tack/dataset/dataset.py
@@ -51,9 +51,16 @@ class Field:
     ``unit`` is the values' physical unit, a string such as ``"m/s"`` kept in
     unyt's canonical spelling (``tack.dataset.units``), or ``None`` for unknown:
     meaning, which kernels never see. Filters that keep a quantity keep its
-    unit; anything else gives ``None`` rather than a guess."""
+    unit; anything else gives ``None`` rather than a guess.
 
-    def __init__(self, space, values, unit=None):
+    ``gradient``, a field on the same topology, is its world gradient where something
+    other than its basis gives it, such as a solver's own reconstruction: 3 values
+    each for a scalar, 3k row by row for a vector of k components (``gradients``'
+    layout). ``world_gradient``, ``world_jet`` and ``gradients`` read it, evaluated
+    at the point asked for as that field's view evaluates it; the Hessian, which a
+    gradient does not determine, is NaN."""
+
+    def __init__(self, space, values, unit=None, gradient=None):
         if not isinstance(space, Space):
             raise TypeError(f"a field's space is a Space on a topology, not {space!r}")
         views.storage(values)                       # refuses what is not an array
@@ -64,6 +71,16 @@ def __init__(self, space, values, unit=None):
         self.values = values
         self._unit = units.canonical(unit)
         _ = self._device                    # refuses values on another device
+        if gradient is not None:
+            if not isinstance(gradient, Field):
+                raise TypeError(f"a field's gradient is a tack.dataset.Field, not {gradient!r}")
+            if gradient.space.topology is not space.topology:
+                raise ValueError("a field's gradient belongs to another topology")
+            width = max(1, arrays.width_of(values) or 0)
+            if (arrays.width_of(gradient.values) or 1) != 3 * width:
+                raise ValueError(f"a field of {width} components has a gradient of {3 * width}, "
+                                 f"not {arrays.width_of(gradient.values) or 1}")
+        self.gradient = gradient
 
     @property
     def unit(self):
@@ -93,7 +110,8 @@ def __repr__(self):
         values = self.values
         stored = repr(values) if isinstance(values, arrays._Implicit) else "stored"
         unit = "" if self._unit is None else f", in {self._unit}"
-        return f"Field({self.space!r}, {arrays.size_of(values)} values, {stored}{unit})"
+        given = "" if self.gradient is None else ", with its gradient"
+        return f"Field({self.space!r}, {arrays.size_of(values)} values, {stored}{unit}{given})"
 
 
 def _is_cells(group):
@@ -130,6 +148,8 @@ def _space_parts(field, group, space_mixin=None, corners_only=False):
         # A geometry's incidence is the domain view's own.
         mixins.extend(space.extra_mixins(group))
     mixins.append(storage)
+    if not space_mixin:
+        mixins.extend(_derivative_mixins(field, group, mixins, attributes))
     # Reading only the points' values -- an order-2 geometry on faces and edges
     # -- needs none of the space's own layout.
     layout = {} if corners_only else space.attributes(group, dtype=_basis_type(field))
@@ -139,6 +159,26 @@ def _space_parts(field, group, space_mixin=None, corners_only=False):
     return mixins, attributes
 
 
+def _derivative_mixins(field, group, mixins, attributes):
+    """Where a field view's ``world_gradient`` and ``world_jet`` come from: the gradient
+    the field was given, a space's own mixin, or else the basis -- the parts of it the
+    view has, since every method of a view resolves when a kernel takes it."""
+    from tack.dataset.algorithms.gradients import _BasisGradient, _BasisJet, _GivenGradient
+
+    given = getattr(field, "gradient", None)
+    if given is not None:
+        attributes.update(given_gradient=given.view(group), no_hessian=float("nan"))
+        return [_GivenGradient]
+    out = []
+    if (any(hasattr(m, "parametric_gradient") for m in mixins)
+            and not any(hasattr(m, "world_gradient") for m in mixins)):
+        out.append(_BasisGradient)
+    if (any(hasattr(m, "parametric_jet") for m in mixins)
+            and not any(hasattr(m, "world_jet") for m in mixins)):
+        out.append(_BasisJet)
+    return out
+
+
 def _basis_type(field):
     """The float type ``field``'s basis is evaluated in: its own (``spaces._float_type``)."""
     from tack.dataset.spaces import _float_type
@@ -778,7 +818,8 @@ def moved(field):
             values = field.values.to_device(target)
             if isinstance(space, L2) and space.layout != field.space.layout:
                 values = _relaid(values, field.space, space)      # "auto" chose another
-            return Field(space, values, field.unit)
+            gradient = None if field.gradient is None else moved(field.gradient)
+            return Field(space, values, field.unit, gradient=gradient)
 
         sets = {name: Set._from_sorted_ids(topology, s.on, s.ids.to_device(target))
                 for name, s in self.sets.items()}
diff --git a/packages/tack-vis/tests/test_given_gradients.py b/packages/tack-vis/tests/test_given_gradients.py
new file mode 100644
index 000000000..b01df215c
--- /dev/null
+++ b/packages/tack-vis/tests/test_given_gradients.py
@@ -0,0 +1,114 @@
+"""A field given its gradient: world_gradient, world_jet and gradients read it instead of
+differentiating the basis."""
+
+import numpy as np
+import pytest
+
+import tack
+import tack.dataset as td
+
+
+@tack.func
+def _u(x):
+    return x[0] * x[0] + 2.0 * x[1] - x[2]
+
+
+@tack.func
+def _v(x):
+    return tack.vec3(x[0] * x[1], x[2], -x[0])
+
+
+@tack.kernel
+def _gradients_at(cells, u, out):
+    for c in cells:
+        out[cells.entity_id(c)] = td.world_gradient(cells, u, c, cells.parametric_center())
+
+
+@tack.kernel
+def _jets_at(cells, u, values, grads, hess):
+    for c in cells:
+        v, g, h = td.world_jet(cells, u, c, cells.parametric_center())
+        e = cells.entity_id(c)
+        values[e] = v
+        grads[e] = g
+        hess[e] = h[0, 0]
+
+
+def _data():
+    return td.uniform_grid((5, 4, 3), origin=(0.1, 0.2, 0.3), spacing=(0.3, 0.25, 0.4))
+
+
+def _vectors(rows):
+    rows = np.asarray(rows, np.float32)
+    out = tack.empty((len(rows),), tack.vec(rows.shape[1], tack.f32))
+    tack.copyto(out, rows)
+    return out
+
+
+def _per_cell(data, rows):
+    return td.Field(td.Constant(data), _vectors(rows))
+
+
+def test_world_gradient_and_gradients_read_a_given_gradient(backend):
+    data = _data()
+    given = np.random.default_rng(1).uniform(-1, 1, (data.num_cells, 3)).astype(np.float32)
+    u = td.interpolate(data, _u, td.H1(data))
+    data.fields["u"] = td.Field(u.space, u.values, gradient=_per_cell(data, given))
+    out = tack.zeros((data.num_cells,), tack.vec3(tack.f32))
+    td.for_each(_gradients_at, data, "cells", data.fields["u"], out)
+    np.testing.assert_allclose(out.numpy(), given)
+    np.testing.assert_allclose(td.gradients(data, "u").values.numpy(), given)
+    plain = td.gradients(data, td.Field(u.space, u.values), at="points").values.numpy()
+    at_points = td.gradients(data, "u", at="points").values.numpy()
+    assert not np.allclose(plain, at_points)
+
+
+def test_world_jet_pairs_the_given_gradient_with_the_value(backend):
+    data = _data()
+    given = np.tile(np.float32([1.0, 2.0, 3.0]), (data.num_cells, 1))
+    u = td.interpolate(data, _u, td.H1(data))
+    field = td.Field(u.space, u.values, gradient=_per_cell(data, given))
+    values = tack.zeros((data.num_cells,), tack.f32)
+    grads = tack.zeros((data.num_cells,), tack.vec3(tack.f32))
+    hess = tack.zeros((data.num_cells,), tack.f32)
+    td.for_each(_jets_at, data, "cells", field, values, grads, hess)
+    centers = td.cell_centers(data).values.numpy()
+    np.testing.assert_allclose(grads.numpy(), given)
+    np.testing.assert_allclose(values.numpy(), _u_numpy(centers), atol=0.05)
+    assert np.isnan(hess.numpy()).all()
+
+
+def _u_numpy(x):
+    return x[:, 0] ** 2 + 2.0 * x[:, 1] - x[:, 2]
+
+
+def test_a_vector_field_gives_its_rows(backend):
+    data = _data()
+    jacobian = np.random.default_rng(2).uniform(-1, 1, (data.num_cells, 9)).astype(np.float32)
+    v = td.interpolate(data, _v, td.H1(data), dtype=tack.vec3(tack.f32))
+    data.fields["v"] = td.Field(v.space, v.values,
+                                gradient=td.Field(td.Constant(data), _vectors(jacobian)))
+    np.testing.assert_allclose(td.gradients(data, "v").values.numpy(), jacobian)
+    q = td.flow_quantities(data, td.gradients(data, "v"))
+    J = jacobian.reshape(-1, 3, 3)
+    np.testing.assert_allclose(q[0].values.numpy(), np.trace(J, axis1=1, axis2=2), rtol=1e-5, atol=1e-6)
+
+
+def test_a_given_gradient_is_checked(backend):
+    data = _data()
+    u = td.interpolate(data, _u, td.H1(data))
+    with pytest.raises(ValueError):
+        td.Field(u.space, u.values, gradient=td.Field(td.Constant(data), tack.zeros((data.num_cells,), tack.f32)))
+    with pytest.raises(ValueError):
+        other = _data()
+        td.Field(u.space, u.values, gradient=_per_cell(other, np.zeros((other.num_cells, 3))))
+    with pytest.raises(TypeError):
+        td.Field(u.space, u.values, gradient=np.zeros(3))
+
+
+def test_point_data_keeps_its_given_gradient(backend):
+    data = _data()
+    u = td.interpolate(data, _u, td.H1(data))
+    given = np.tile(np.float32([4.0, 5.0, 6.0]), (data.num_cells, 1))
+    data.fields["u"] = td.Field(td.Values(data, "points"), u.values, gradient=_per_cell(data, given))
+    np.testing.assert_allclose(td.gradients(data, "u").values.numpy(), given)

From 2a433132fcfa538a74e2804e3214d3d6baeeb6b6 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Wed, 4 Nov 2026 20:43:06 +0000
Subject: [PATCH 05/12] SPDX header for test_given_gradients.py

---
 packages/tack-vis/tests/test_given_gradients.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/packages/tack-vis/tests/test_given_gradients.py b/packages/tack-vis/tests/test_given_gradients.py
index 08b070a7d..4134e5d18 100644
--- a/packages/tack-vis/tests/test_given_gradients.py
+++ b/packages/tack-vis/tests/test_given_gradients.py
@@ -1,3 +1,6 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
 """A field given its gradient: world_gradient, world_jet and gradients read it instead of
 differentiating the basis."""
 

From a273ae699289812402244f0b55551d33460338c8 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Wed, 4 Nov 2026 20:43:06 +0000
Subject: [PATCH 06/12] SPDX header for test_given_gradients.py

---
 packages/tack-vis/tests/test_given_gradients.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/packages/tack-vis/tests/test_given_gradients.py b/packages/tack-vis/tests/test_given_gradients.py
index b01df215c..55c8d2d7a 100644
--- a/packages/tack-vis/tests/test_given_gradients.py
+++ b/packages/tack-vis/tests/test_given_gradients.py
@@ -1,3 +1,6 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
 """A field given its gradient: world_gradient, world_jet and gradients read it instead of
 differentiating the basis."""
 

From fce2b4351fdb61c04415a6093c0d054d5bfd72de Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Wed, 4 Nov 2026 20:43:06 +0000
Subject: [PATCH 07/12] SPDX header for test_field_resolver.py

---
 packages/tack-vis/tests/test_field_resolver.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/packages/tack-vis/tests/test_field_resolver.py b/packages/tack-vis/tests/test_field_resolver.py
index ef04dc5f0..fc34f95b9 100644
--- a/packages/tack-vis/tests/test_field_resolver.py
+++ b/packages/tack-vis/tests/test_field_resolver.py
@@ -1,3 +1,6 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
 """Field arguments as text a registered resolver reads."""
 
 import numpy as np

From 97ae286f4d3b1c1b3e0ebf1bf9473aab12022d67 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Wed, 4 Nov 2026 21:42:45 +0000
Subject: [PATCH 08/12] SPDX header for _vtk_algorithms_mpi.py

---
 packages/tack-vis/tests/_vtk_algorithms_mpi.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/packages/tack-vis/tests/_vtk_algorithms_mpi.py b/packages/tack-vis/tests/_vtk_algorithms_mpi.py
index 6029609d4..60c76f411 100644
--- a/packages/tack-vis/tests/_vtk_algorithms_mpi.py
+++ b/packages/tack-vis/tests/_vtk_algorithms_mpi.py
@@ -1,3 +1,6 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
 """``test_vtk_algorithms``' check across MPI ranks: ``mpiexec -n 3 python
 _vtk_algorithms_mpi.py comm|vtk``. Each rank holds one slab of VTK's wavelet, the
 slabs sharing their planes of points, in a ``vtkPartitionedDataSet`` of one

From d6c4806a2d5a4ec2e06eed38c6f8dc824c80b8b6 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Wed, 4 Nov 2026 21:53:07 +0000
Subject: [PATCH 09/12] statistics, entropy and integrate take fields through
 the field resolver

A field evaluated where it is read was reduced by its placeholder values; now by
its stored values, and integrate reads it through its view.
---
 changelog.d/reductions-take-fields.md         |  7 +++
 .../src/tack/dataset/algorithms/_common.py    |  9 +++
 .../src/tack/dataset/algorithms/entropy.py    |  6 +-
 .../src/tack/dataset/algorithms/integrate.py  | 55 +++++++++++++++----
 .../src/tack/dataset/algorithms/statistics.py |  7 ++-
 packages/tack-vis/src/tack/dataset/ghosts.py  |  8 +--
 .../tack-vis/src/tack/dataset/partitioned.py  | 10 +++-
 .../tack-vis/tests/test_field_resolver.py     | 22 ++++++++
 8 files changed, 100 insertions(+), 24 deletions(-)
 create mode 100644 changelog.d/reductions-take-fields.md

diff --git a/changelog.d/reductions-take-fields.md b/changelog.d/reductions-take-fields.md
new file mode 100644
index 000000000..aa344c705
--- /dev/null
+++ b/changelog.d/reductions-take-fields.md
@@ -0,0 +1,7 @@
+### Changed
+
+- `statistics`, `entropy` and `integrate` take any field argument: a name, a
+  `Field`, or text a registered resolver reads. A field evaluated where it is read
+  is reduced by its stored values (`owned_values` stores it) and integrated
+  through its view by quadrature. `integrate` takes names, text, or a mapping of
+  names to any of these.
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_common.py b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
index 297f43166..2503aecd1 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -175,6 +175,15 @@ def _field(data, field):
     return field
 
 
+def _stored_field(data, field):
+    """``_field``, with a field evaluated where it is read (an expression's) replaced
+    by its values, as its space's ``stored`` gives them: for an algorithm that
+    reduces a field's values rather than reading it through its view."""
+    field = _field(data, field)
+    stored = getattr(field.space, "stored", None)
+    return stored(field, data) if callable(stored) else field
+
+
 def _with_basis(data, field):
     """Values on points or cells read as the space they are: order-1 ``H1`` and
     ``Constant``, so every algorithm that takes a field with a basis takes them
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/entropy.py b/packages/tack-vis/src/tack/dataset/algorithms/entropy.py
index 5f81eac6f..0303209b7 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/entropy.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/entropy.py
@@ -6,7 +6,7 @@
 import numpy as np
 
 from tack.arrays import size_of, width_of
-from tack.dataset.algorithms._common import _field
+from tack.dataset.algorithms._common import _stored_field
 from tack.dataset.algorithms.statistics import _over_pieces
 from tack.dataset.ghosts import owned_values
 from tack.dataset.partitioned import PartitionedDataSet, fields_of
@@ -19,7 +19,7 @@
 def _extrema(data, field):
     """The minimum and maximum of the values of one dataset's scalar field that it
     owns, or ``None`` for no values."""
-    field = _field(data, field)
+    field = _stored_field(data, field)
     if width_of(field.values):
         raise TypeError("entropy takes a scalar field")
     values = owned_values(data, field)
@@ -32,7 +32,7 @@ def _extrema(data, field):
 def _counts(data, field, bins, value_range, weights=None):
     from tack.algorithms import histogram
 
-    field = _field(data, field)
+    field = _stored_field(data, field)
     values = owned_values(data, field)
     if weights is not None:
         from tack.dataset.algorithms.statistics import _weights
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/integrate.py b/packages/tack-vis/src/tack/dataset/algorithms/integrate.py
index fea1ab349..02ba2c5f1 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/integrate.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/integrate.py
@@ -13,7 +13,7 @@
 import tack
 from tack import arrays
 from tack.dataset import shapes, units
-from tack.dataset.algorithms._common import _field, _with_basis
+from tack.dataset.algorithms._common import _with_basis
 from tack.dataset.algorithms.cell_geometry import _signed_twice_area
 from tack.dataset.dataset import Field, is_field_key
 from tack.dataset.spaces import H1, L2, Constant, Values
@@ -206,7 +206,10 @@ def quadrature(shape, n):
 
 
 def _order(space):
-    """The polynomial order a space's values carry along an axis."""
+    """The polynomial order a space's values carry along an axis: for one evaluated
+    where it is read, twice its inputs' (a product of two)."""
+    if getattr(space, "evaluated", False):
+        return 2 * int(getattr(space, "default_subdivision", 1))
     if isinstance(space, (Constant, Values)):
         return 0 if isinstance(space, Constant) or space.on == "cells" else 1
     return max(int(getattr(space, "order", 1)), 0)
@@ -312,17 +315,23 @@ def _no_reference_values(data, field):
 
 
 def _piece(data, names):
-    """One dataset's owned cells' measure and integrals, summed in cell order (f64)."""
+    """One dataset's owned cells' measure and integrals, summed in cell order (f64):
+    ``names`` as ``_names`` gives them."""
+    from tack.dataset.algorithms._common import _field
     from tack.dataset.algorithms.cell_geometry import cell_geometry
     from tack.dataset.ghosts import owned_values, visible
 
+    given = data
     data = visible(data)
     measures = Field(Values(data, "cells"), _top_cells(data, arrays.materialize(
         cell_geometry(data).measure.values)))
     measure = float(np.sum(np.asarray(owned_values(data, measures).numpy(), np.float64)))
     values = {}
-    for name in names:
-        field = data.fields[name]
+    for name, spec in names:
+        field = visible(given, spec)[1] if isinstance(spec, Field) else _field(data, spec)
+        stored = getattr(field.space, "stored", None)
+        if callable(stored) and not getattr(data.topology, "reference_cells", True):
+            field = stored(field, data)
         if not getattr(data.topology, "reference_cells", True):
             space = _with_basis(data, field).space
             if _on_points(space):
@@ -375,29 +384,51 @@ def integrate(data, fields=None):
 
         pieces = list(data)
         names = _names(pieces[0], fields) if pieces else []
+        if any(isinstance(spec, Field) for _, spec in names):
+            raise TypeError("integrate over a PartitionedDataSet takes fields by name or "
+                            "text: a Field belongs to one partition's topology")
         parts = _over_pieces(data, lambda piece: _piece(piece, names))
         measure = math.fsum(m for m, _ in parts)
-        values = {name: np.sum([v[name] for _, v in parts], axis=0) for name in names}
+        values = {name: np.sum([v[name] for _, v in parts], axis=0) for name, _ in names}
         first = pieces[0] if pieces else None
     else:
         names = _names(data, fields)
         measure, values = _piece(data, names)
         first = data
-    out_units = {name: units.derived((first.fields[name].unit, 1),
+    out_units = {name: units.derived((_unit_of(first, spec), 1),
                                      (_measure_unit(first), 1)) if first is not None else None
-                 for name in names}
+                 for name, spec in names}
     return Integrals(measure=measure, measure_unit=_measure_unit(first) if first else None,
                      values={name: (float(v) if np.ndim(v) == 0 else tuple(map(float, v)))
                              for name, v in values.items()},
                      units=out_units)
 
 
+def _unit_of(data, spec):
+    from tack.dataset.algorithms._common import _field
+
+    return spec.unit if isinstance(spec, Field) else _field(data, spec).unit
+
+
 def _names(data, fields):
+    """``[(name, spec)]``: every stored field with a basis or on the cells (and every
+    one evaluated where it is read) by default; or the names, or text a
+    registered resolver reads, each its own name; or a mapping of names to any
+    of these or to ``Field``s."""
+    from tack.dataset.algorithms._common import _field
+
     if fields is None:
-        return [name for name, field in data.fields.items()
+        return [(name, name) for name, field in data.fields.items()
                 if isinstance(field.space, (H1, L2, Constant))
+                or getattr(field.space, "evaluated", False)
                 or (isinstance(field.space, Values) and field.space.on in ("points", "cells"))]
-    names = [fields] if is_field_key(fields) else list(fields)
-    for name in names:
-        _field(data, name)
+    if hasattr(fields, "items"):
+        names = list(fields.items())
+    else:
+        entries = [fields] if is_field_key(fields) or isinstance(fields, Field) else list(fields)
+        if any(isinstance(e, Field) for e in entries):
+            raise TypeError("integrate names each Field it is given: pass {name: field}")
+        names = [(e, e) for e in entries]
+    for _, spec in names:
+        _field(data, spec)
     return names
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/statistics.py b/packages/tack-vis/src/tack/dataset/algorithms/statistics.py
index c51d873a4..0ff5bd271 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/statistics.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/statistics.py
@@ -10,7 +10,7 @@
 
 import tack
 from tack.arrays import size_of, width_of
-from tack.dataset.algorithms._common import _converted, _field, _floating_dtype
+from tack.dataset.algorithms._common import _converted, _field, _floating_dtype, _stored_field
 from tack.dataset.ghosts import owned_values
 from tack.dataset.partitioned import PartitionedDataSet, fields_of
 from tack.runtime.dispatch import on_inputs_device
@@ -107,7 +107,7 @@ def _moments(data, field, weights=None):
     over the values it owns (not its ghosts), the moments about their own mean
     in a second pass, weighted by ``weights`` (``W`` their sum, ``W2`` the sum of
     their squares; both ``n`` without); ``None`` for no values."""
-    field = _field(data, field)
+    field = _stored_field(data, field)
     if width_of(field.values):
         raise TypeError("statistics takes a scalar field")
     source = owned_values(data, field)
@@ -133,12 +133,13 @@ def _weights(data, field, weights):
     from tack.dataset.algorithms.cell_geometry import cell_geometry
     from tack.dataset.algorithms.ghost_cells import _kind
 
+    field = _stored_field(data, field)
     kind = _kind(field)
     if isinstance(weights, str) and weights == "volume":
         if kind != "cells":
             raise TypeError("volume weights are a cell's measure: they take a cell field")
         return cell_geometry(data).measure
-    weights = _field(data, weights)
+    weights = _stored_field(data, weights)
     if _kind(weights) != kind or width_of(weights.values):
         raise TypeError("weights are a scalar field on the same entities as the field "
                         "(points or cells), or \"volume\"")
diff --git a/packages/tack-vis/src/tack/dataset/ghosts.py b/packages/tack-vis/src/tack/dataset/ghosts.py
index 79a6ada4b..4b40b6f9a 100644
--- a/packages/tack-vis/src/tack/dataset/ghosts.py
+++ b/packages/tack-vis/src/tack/dataset/ghosts.py
@@ -223,17 +223,17 @@ def _kept(flags, slots, ids):
 
 @on_inputs_device
 def owned_values(data, field):
-    """The values of ``field`` (a ``Field`` of ``data``, or its name) that ``data`` owns: all of
+    """The values of ``field`` (a ``Field`` of ``data``, its name, or text a registered
+    resolver reads; one evaluated where it is read, by its stored values) that ``data`` owns: all of
     them when it has no marks for where they live, else those of cells and
     points not marked a duplicate or hidden, materialized. Fields on faces or
     edges, and order-2 fields, have no marks of their own yet; on a dataset
     with ghost marks they are refused rather than counted twice."""
+    from tack.dataset.algorithms._common import _stored_field
     from tack.dataset.carry import _take
-    from tack.dataset.dataset import is_field_key
     from tack.dataset.spaces import H1, L2, Constant, Values
 
-    if is_field_key(field):
-        field = data.fields[field]
+    field = _stored_field(data, field)
     space = field.space
     on_points = (isinstance(space, H1) and space.order == 1) or (
         isinstance(space, Values) and space.on == "points")
diff --git a/packages/tack-vis/src/tack/dataset/partitioned.py b/packages/tack-vis/src/tack/dataset/partitioned.py
index 13ab9a67f..c6c793844 100644
--- a/packages/tack-vis/src/tack/dataset/partitioned.py
+++ b/packages/tack-vis/src/tack/dataset/partitioned.py
@@ -78,10 +78,16 @@ def fields_of(data, field, what):
     if not is_field_key(field):
         raise TypeError(f"{what} of a PartitionedDataSet takes a field's name (or name "
                         "and space): a Field belongs to one partition's topology")
+    from tack.dataset.algorithms._common import _field
+
     for i, piece in enumerate(data):
-        if field not in piece.fields:
+        if field in piece.fields:
+            continue
+        try:
+            _field(piece, field)        # text a registered resolver reads
+        except KeyError:
             raise KeyError(f"partition {i} has no field {field!r}: it has "
-                           f"{', '.join(map(repr, piece.fields)) or 'none'}")
+                           f"{', '.join(map(repr, piece.fields)) or 'none'}") from None
     return field
 
 
diff --git a/packages/tack-vis/tests/test_field_resolver.py b/packages/tack-vis/tests/test_field_resolver.py
index fc34f95b9..698f18090 100644
--- a/packages/tack-vis/tests/test_field_resolver.py
+++ b/packages/tack-vis/tests/test_field_resolver.py
@@ -41,3 +41,25 @@ def test_algorithms_take_text_a_resolver_reads(backend, twice):
                                td.gradients(data, doubled).values.numpy())
     with pytest.raises(KeyError):
         td.contour(data, "missing", 1.0)
+
+
+def test_reductions_take_text_a_resolver_reads(backend, twice):
+    data = td.uniform_grid((5, 4, 3), origin=(0.1, 0.2, 0.3), spacing=(0.3, 0.25, 0.4))
+    data.fields["u"] = td.interpolate(data, _u, td.H1(data))
+    data.fields["2u"] = td.Field(data.fields["u"].space, tack.array(2.0 * data.fields["u"].values.numpy()))
+    by_text = td.statistics(data, "twice u")
+    stored = td.statistics(data, "2u")
+    assert (by_text.n, by_text.min, by_text.max) == (stored.n, stored.min, stored.max)
+    np.testing.assert_allclose(by_text.mean, stored.mean)
+    assert td.entropy(data, "twice u") == td.entropy(data, "2u")
+    a = td.integrate(data, ["twice u"])
+    b = td.integrate(data, {"doubled": "2u", "field": data.fields["2u"]})
+    np.testing.assert_allclose(a.values["twice u"], b.values["doubled"])
+    np.testing.assert_allclose(b.values["field"], b.values["doubled"])
+    with pytest.raises(TypeError, match="names each Field"):
+        td.integrate(data, [data.fields["2u"]])
+    pieces = td.PartitionedDataSet([data])
+    np.testing.assert_allclose(td.statistics(pieces, "twice u").mean, stored.mean)
+    np.testing.assert_allclose(td.integrate(pieces, ["twice u"]).values["twice u"],
+                               a.values["twice u"])
+

From 7dcbff3f038557509570d64651aa47af4ceca919 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Wed, 4 Nov 2026 22:44:30 +0000
Subject: [PATCH 10/12] Derived fields: one name, read at the location asked
 for

register_derived_field(name, compute, at) lets a name stand, at points, cells or
faces, for a field a dataset makes from itself where it stores none of that
name there; named_field(data, name, at) gives the stored one there first, then
the derived one, made once and kept while the geometry is the same. Every field
argument reads them. Tack registers normals (a surface's at its points and
cells, every mesh's at its faces), measures (cells and faces) and centroids.
---
 changelog.d/derived-fields.md                 | 10 +++
 docs/vis-guide/reference/data-model.md        |  6 ++
 .../tack-vis/src/tack/dataset/__init__.py     | 12 ++-
 .../src/tack/dataset/algorithms/_common.py    | 77 ++++++++++++++++++-
 .../src/tack/dataset/algorithms/derived.py    | 43 +++++++++++
 .../tack-vis/tests/test_derived_fields.py     | 70 +++++++++++++++++
 6 files changed, 215 insertions(+), 3 deletions(-)
 create mode 100644 changelog.d/derived-fields.md
 create mode 100644 packages/tack-vis/src/tack/dataset/algorithms/derived.py
 create mode 100644 packages/tack-vis/tests/test_derived_fields.py

diff --git a/changelog.d/derived-fields.md b/changelog.d/derived-fields.md
new file mode 100644
index 000000000..23e3a85cb
--- /dev/null
+++ b/changelog.d/derived-fields.md
@@ -0,0 +1,10 @@
+### Added
+
+- Derived fields: a name a dataset stores no field under, at the location it is read,
+  may stand for a field the dataset makes from itself (`register_derived_field(name,
+  compute, at)`), made the first time it is asked for and kept while the geometry
+  is the same. One name is one quantity: `normals` is a surface's normals at its
+  points and cells and every mesh's face normals at its faces; `measures` each
+  cell's volume, area or length and each face's area; `centroids` each cell's.
+  Every field argument reads them (`named_field`), and so do field expressions,
+  so `integrate(surface, ["dot(U, normals)"])` is a flux with no filter first.
diff --git a/docs/vis-guide/reference/data-model.md b/docs/vis-guide/reference/data-model.md
index e7874603d..ac22b254a 100644
--- a/docs/vis-guide/reference/data-model.md
+++ b/docs/vis-guide/reference/data-model.md
@@ -207,6 +207,12 @@ do, so pieces sharing an edge compute it alike.
 
 ::: tack.dataset.register_field_resolver
 
+::: tack.dataset.register_derived_field
+
+::: tack.dataset.named_field
+
+::: tack.dataset.derived_names
+
 ::: tack.dataset.keep_marks
 
 ::: tack.dataset.drop_marks
diff --git a/packages/tack-vis/src/tack/dataset/__init__.py b/packages/tack-vis/src/tack/dataset/__init__.py
index 9b5b93156..87c433f93 100644
--- a/packages/tack-vis/src/tack/dataset/__init__.py
+++ b/packages/tack-vis/src/tack/dataset/__init__.py
@@ -14,6 +14,7 @@
 one module each, all re-exported here.
 """
 
+import tack.dataset.algorithms.derived
 from tack.dataset import (
     algorithms,
     implicit,
@@ -112,7 +113,13 @@
     voxelize,
     warp,
 )
-from tack.dataset.algorithms._common import edge_crossing, register_field_resolver
+from tack.dataset.algorithms._common import (
+    derived_names,
+    edge_crossing,
+    named_field,
+    register_derived_field,
+    register_field_resolver,
+)
 from tack.dataset.algorithms.gradients import world_gradient, world_hessian, world_jet
 from tack.dataset.amr import AMRDataSet
 from tack.dataset.carry import Interpolate, Pieces, Same, Take, carry
@@ -232,6 +239,7 @@
     "cross",
     "cylindrical",
     "deposit",
+    "derived_names",
     "divergence",
     "dot",
     "drop_marks",
@@ -272,6 +280,7 @@
     "mask_points",
     "merge_pieces",
     "merge_points",
+    "named_field",
     "net_flux",
     "orient",
     "outliers",
@@ -285,6 +294,7 @@
     "probe",
     "project",
     "rectilinear_grid",
+    "register_derived_field",
     "register_field_resolver",
     "remove_ghosts",
     "shapes",
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_common.py b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
index 2503aecd1..a921ccf62 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -140,6 +140,79 @@ def _same_space(*fields):
 
 
 _RESOLVERS = []
+_DERIVED = {}             # name -> {location: compute}
+_LOCATIONS = ("points", "cells", "faces")
+
+
+def register_derived_field(name, compute, at):
+    """Lets ``name`` stand, at ``at`` (``"points"``, ``"cells"`` or ``"faces"``), for a
+    field a dataset makes from itself when it stores none of that name there:
+    ``compute(data)`` gives the ``Field`` on ``data``, made the first time it is
+    asked for and kept with the dataset while its geometry is the same. One name
+    is one quantity, registered at each location it has: Tack registers
+    ``normals`` (a surface's at its points and cells, ``surface_normals``'; every
+    mesh's at its faces, ``face_geometry``'s), ``measures`` (each cell's volume, area
+    or length; each face's area) and ``centroids`` (each cell's)."""
+    if at not in _LOCATIONS:
+        raise ValueError(f"at is one of {_LOCATIONS}, not {at!r}")
+    _DERIVED.setdefault(name, {})[at] = compute
+
+
+def derived_names():
+    """The names ``register_derived_field`` has registered."""
+    return tuple(_DERIVED)
+
+
+def _location(field):
+    from tack.dataset.spaces import H1, L2, Constant, Values
+
+    space = field.space
+    if isinstance(space, H1) or (isinstance(space, Values) and space.on == "points"):
+        return "points"
+    if isinstance(space, (Constant, L2)) or (isinstance(space, Values) and space.on == "cells"):
+        return "cells"
+    if isinstance(space, Values) and space.on == "faces":
+        return "faces"
+    return None
+
+
+def _derived(data, name, at):
+    made = data.__dict__.setdefault("_derived_fields", {})
+    hit = made.get((name, at))
+    if hit is not None and hit[0] is data.geometry:
+        return hit[1]
+    field = _DERIVED[name][at](data)
+    made[(name, at)] = (data.geometry, field)
+    return field
+
+
+def named_field(data, name, at=None):
+    """The field ``name`` names on ``data``, read at ``at`` (``"points"``, ``"cells"``,
+    ``"faces"``, or None for wherever it is). A field stored under the name there
+    comes first, then the one the dataset derives there
+    (``register_derived_field``), then one stored under the name elsewhere. With
+    no ``at``, the stored field, else the derived one at the first location it is
+    registered at that the dataset has. A ``KeyError`` for none."""
+    if not isinstance(name, str) or name not in _DERIVED:
+        return data.fields[name]
+    stored = data.fields.named(name)
+    places = _DERIVED[name]
+    if at is not None:
+        for field in stored:
+            if _location(field) == at:
+                return field
+        if at in places:
+            return _derived(data, name, at)
+        return data.fields[name]
+    if stored:
+        return data.fields[name]
+    reasons = []
+    for where in places:
+        try:
+            return _derived(data, name, where)
+        except (ValueError, TypeError, NotImplementedError) as e:
+            reasons.append(f"at {where}: {e}")
+    raise KeyError(f"{name!r} is derived from no part of this dataset ({'; '.join(reasons)})")
 
 
 def register_field_resolver(resolve):
@@ -152,8 +225,8 @@ def register_field_resolver(resolve):
 
 
 def _resolved(data, text):
-    if not isinstance(text, str) or text in data.fields:
-        return data.fields[text]
+    if not isinstance(text, str) or text in data.fields or text in _DERIVED:
+        return named_field(data, text)
     for resolve in reversed(_RESOLVERS):
         field = resolve(data, text)
         if field is not None:
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/derived.py b/packages/tack-vis/src/tack/dataset/algorithms/derived.py
new file mode 100644
index 000000000..8888b9d56
--- /dev/null
+++ b/packages/tack-vis/src/tack/dataset/algorithms/derived.py
@@ -0,0 +1,43 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""The fields a dataset derives from its own geometry where it stores none of their
+names (``register_derived_field``): ``normals``, ``measures`` and ``centroids``,
+each at the locations it has."""
+
+from tack.dataset.algorithms._common import register_derived_field
+
+
+def _surface_normals(name, on):
+    def compute(data):
+        from tack.dataset.algorithms.surface_normals import surface_normals
+        from tack.dataset.dataset import Field
+        from tack.dataset.spaces import Values
+
+        made = surface_normals(data).fields[name]
+        return Field(Values(data, on), made.values, made.unit)
+    return compute
+
+
+def _face_geometry(part):
+    def compute(data):
+        from tack.dataset.algorithms.face_geometry import face_geometry
+
+        return getattr(face_geometry(data), part)
+    return compute
+
+
+def _cell_geometry(part):
+    def compute(data):
+        from tack.dataset.algorithms.cell_geometry import cell_geometry
+
+        return getattr(cell_geometry(data), part)
+    return compute
+
+
+register_derived_field("normals", _surface_normals("normals", "points"), at="points")
+register_derived_field("normals", _surface_normals("cell_normals", "cells"), at="cells")
+register_derived_field("normals", _face_geometry("normal"), at="faces")
+register_derived_field("measures", _cell_geometry("measure"), at="cells")
+register_derived_field("measures", _face_geometry("area"), at="faces")
+register_derived_field("centroids", _cell_geometry("centroid"), at="cells")
diff --git a/packages/tack-vis/tests/test_derived_fields.py b/packages/tack-vis/tests/test_derived_fields.py
new file mode 100644
index 000000000..59ab218d8
--- /dev/null
+++ b/packages/tack-vis/tests/test_derived_fields.py
@@ -0,0 +1,70 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""Fields a dataset derives from itself where it stores none of the name: one name,
+read at the location asked for."""
+
+import numpy as np
+import pytest
+
+import tack
+import tack.dataset as td
+from tack.dataset.algorithms import _common
+
+
+def _grid():
+    return td.uniform_grid((4, 3, 3), spacing=(0.5, 0.5, 0.5))
+
+
+def test_normals_at_each_location(backend):
+    box = td.extract_surface(_grid())
+    made = td.surface_normals(box)
+    np.testing.assert_array_equal(td.named_field(box, "normals", "points").values.numpy(),
+                                  made.fields["normals"].values.numpy())
+    np.testing.assert_array_equal(td.named_field(box, "normals", "cells").values.numpy(),
+                                  made.fields["cell_normals"].values.numpy())
+    solid = _grid()
+    np.testing.assert_array_equal(td.named_field(solid, "normals", "faces").values.numpy(),
+                                  td.face_geometry(solid).normal.values.numpy())
+    assert "normals" not in box.fields
+    assert td.named_field(box, "normals", "cells") is td.named_field(box, "normals", "cells")
+    np.testing.assert_allclose(td.named_field(box, "measures", "cells").values.numpy().sum(), 8.0,
+                               rtol=1e-6)
+
+
+def test_a_field_stored_at_the_location_comes_first(backend):
+    box = td.extract_surface(_grid())
+    mine = td.Field(td.Values(box, "cells"), td.surface_normals(box).fields["cell_normals"].values)
+    box.fields["normals"] = mine
+    assert td.named_field(box, "normals", "cells") is mine
+    # Stored at the cells only: at the points, the derived ones.
+    assert td.named_field(box, "normals", "points") is not mine
+    assert td.named_field(box, "normals") is mine
+
+
+def test_no_location_takes_the_first_the_dataset_has(backend):
+    solid = _grid()
+    # A solid has no surface normals at its points or cells: its faces'.
+    assert td.named_field(solid, "normals").space.on == "faces"
+
+
+def test_registered_names_and_misses(backend):
+    calls = []
+
+    def compute(data):
+        calls.append(1)
+        return td.Field(td.Values(data, "cells"), tack.array(np.arange(data.num_cells, dtype=np.float32)))
+
+    td.register_derived_field("cell_index", compute, at="cells")
+    try:
+        box = td.extract_surface(_grid())
+        assert "cell_index" in td.derived_names()
+        assert td.statistics(box, "cell_index").max == box.num_cells - 1
+        td.statistics(box, "cell_index")
+        assert len(calls) == 1
+        with pytest.raises(KeyError):
+            td.named_field(box, "nosuch")
+        with pytest.raises(ValueError):
+            td.register_derived_field("x", compute, at="edges")
+    finally:
+        del _common._DERIVED["cell_index"]

From a8001c6388280a76d75d87c0f2eba482e08a94cc Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Thu, 5 Nov 2026 00:02:50 +0000
Subject: [PATCH 11/12] Derived fields are part of the geometry

Normals, measures and centroids are kept on the geometry, so every dataset
on the same geometry shares them. extract_surface of the same geometry keeps
its surface geometry, so a flux over many steps derives the normals once.
---
 changelog.d/derived-fields.md                 |  5 ++--
 .../src/tack/dataset/algorithms/_common.py    | 19 ++++++------
 .../dataset/algorithms/extract_surface.py     | 30 ++++++++++++++++---
 .../tack-vis/tests/test_derived_fields.py     | 11 +++++++
 4 files changed, 50 insertions(+), 15 deletions(-)

diff --git a/changelog.d/derived-fields.md b/changelog.d/derived-fields.md
index 23e3a85cb..f96b7c158 100644
--- a/changelog.d/derived-fields.md
+++ b/changelog.d/derived-fields.md
@@ -2,8 +2,9 @@
 
 - Derived fields: a name a dataset stores no field under, at the location it is read,
   may stand for a field the dataset makes from itself (`register_derived_field(name,
-  compute, at)`), made the first time it is asked for and kept while the geometry
-  is the same. One name is one quantity: `normals` is a surface's normals at its
+  compute, at)`), made the first time it is asked for and kept on the geometry.
+  `extract_surface` of the same geometry keeps its surface geometry too, so a
+  flux over many steps derives the normals once. One name is one quantity: `normals` is a surface's normals at its
   points and cells and every mesh's face normals at its faces; `measures` each
   cell's volume, area or length and each face's area; `centroids` each cell's.
   Every field argument reads them (`named_field`), and so do field expressions,
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_common.py b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
index a921ccf62..e7d41f6ef 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -147,8 +147,10 @@ def _same_space(*fields):
 def register_derived_field(name, compute, at):
     """Lets ``name`` stand, at ``at`` (``"points"``, ``"cells"`` or ``"faces"``), for a
     field a dataset makes from itself when it stores none of that name there:
-    ``compute(data)`` gives the ``Field`` on ``data``, made the first time it is
-    asked for and kept with the dataset while its geometry is the same. One name
+    ``compute(data)`` gives the ``Field`` on ``data``, a function of its geometry:
+    made the first time it is asked for and kept on the geometry, so every
+    dataset on the same geometry shares it (``extract_surface`` of the same
+    geometry gives the same surface geometry, and so the same normals). One name
     is one quantity, registered at each location it has: Tack registers
     ``normals`` (a surface's at its points and cells, ``surface_normals``'; every
     mesh's at its faces, ``face_geometry``'s), ``measures`` (each cell's volume, area
@@ -177,13 +179,12 @@ def _location(field):
 
 
 def _derived(data, name, at):
-    made = data.__dict__.setdefault("_derived_fields", {})
-    hit = made.get((name, at))
-    if hit is not None and hit[0] is data.geometry:
-        return hit[1]
-    field = _DERIVED[name][at](data)
-    made[(name, at)] = (data.geometry, field)
-    return field
+    # Part of the geometry, kept on it: any dataset on the same geometry (and so
+    # the same topology) shares them.
+    made = data.geometry.__dict__.setdefault("_derived", {})
+    if (name, at) not in made:
+        made[(name, at)] = _DERIVED[name][at](data)
+    return made[(name, at)]
 
 
 def named_field(data, name, at=None):
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py b/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py
index 1fcc50b82..ab518cbfd 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py
@@ -142,7 +142,7 @@ def _fresh_boundary(data, fields):
         return None
     if topology._boundary is not None:
         surface, cells, points = topology._boundary
-        return carry(data, surface, points=Take(points), cells=Pieces(cells), fields=fields)
+        return _on_boundary(data, surface, fields, points=Take(points), cells=Pieces(cells))
     groups = [g for g in topology.groups() if g.shape.NUM_FACES and g.shape.DIMENSION == 3]
     starts = np.concatenate([[0], np.cumsum([g.count * g.shape.NUM_FACES
                                              for g in groups])]).astype(int)
@@ -204,7 +204,21 @@ def _fresh_boundary(data, fields):
     count, points = _used_points(connectivity, data.num_points, idt)
     surface = UnstructuredTopology(types, offsets, connectivity, num_points=count)
     topology._boundary = (surface, cells, points)
-    return carry(data, surface, points=Take(points), cells=Pieces(cells), fields=fields)
+    return _on_boundary(data, surface, fields, points=Take(points), cells=Pieces(cells))
+
+
+def _on_boundary(data, surface, fields, **maps):
+    """The boundary ``surface`` with ``data``'s fields carried by ``maps``. Its
+    geometry is part of ``data``'s: kept on it, so a repeated call on the same
+    geometry gives the same surface geometry, and what is derived from that
+    (``normals``, ``measures``: ``register_derived_field``) is derived once."""
+    kept = data.geometry.__dict__.get("_boundary")
+    if kept is not None and kept[0] is surface:
+        out = carry(data, surface, geometry=kept[1], fields=fields, **kept[2])
+    else:
+        out = carry(data, surface, fields=fields, **maps)
+    data.geometry.__dict__["_boundary"] = (surface, out.geometry, maps)
+    return out
 
 
 @on_inputs_device
@@ -260,6 +274,12 @@ def extract_surface(data, faces=None, fields=None):
         fresh = _fresh_boundary(data, fields)
         if fresh is not None:
             return fresh
+        kept = data.geometry.__dict__.get("_boundary")
+        if (kept is not None and "faces_to_cells" in kept[2]
+                and not (data.topology.dimension == 3 and _has_order_two(data, fields))):
+            return _on_boundary(data, kept[0], fields, **kept[2])
+    boundary = faces is None
+    if boundary:
         from tack.dataset.algorithms.boundary_faces import boundary_faces
         faces = boundary_faces(data)
     ids = _face_ids(data, faces)
@@ -295,8 +315,10 @@ def extract_surface(data, faces=None, fields=None):
         return _on_used_points(out)
     count, points = _used_points(connectivity, data.num_points, idt)
     surface = UnstructuredTopology(types, offsets, connectivity, num_points=count)
-    return carry(data, surface, points=Take(points), cells=Pieces(cells),
-                 faces_to_cells=Take(ids), fields=fields)
+    maps = {"points": Take(points), "cells": Pieces(cells), "faces_to_cells": Take(ids)}
+    if boundary:
+        return _on_boundary(data, surface, fields, **maps)
+    return carry(data, surface, fields=fields, **maps)
 
 
 @tack.kernel
diff --git a/packages/tack-vis/tests/test_derived_fields.py b/packages/tack-vis/tests/test_derived_fields.py
index 59ab218d8..8bab79dc8 100644
--- a/packages/tack-vis/tests/test_derived_fields.py
+++ b/packages/tack-vis/tests/test_derived_fields.py
@@ -48,6 +48,17 @@ def test_no_location_takes_the_first_the_dataset_has(backend):
     assert td.named_field(solid, "normals").space.on == "faces"
 
 
+def test_a_boundary_of_the_same_geometry_derives_its_normals_once(backend):
+    solid = td.tetrahedralize(_grid())
+    for volume in (_grid(), solid):
+        first = td.extract_surface(volume)
+        again = td.extract_surface(volume)
+        assert again.geometry is first.geometry
+        assert td.named_field(again, "normals", "cells") is td.named_field(first, "normals", "cells")
+    moved = td.extract_surface(_grid())
+    assert td.named_field(moved, "normals", "cells") is not td.named_field(first, "normals", "cells")
+
+
 def test_registered_names_and_misses(backend):
     calls = []
 

From 5740b84f62f2b024f7da48038dec3bea3bda4058 Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Thu, 5 Nov 2026 00:49:48 +0000
Subject: [PATCH 12/12] Mesh metrics live with the mesh

What is computed from the mesh alone (derived fields, the boundary's
geometry) is kept with it: shared by every dataset on the same geometry, and
made again when the mesh changes, by another geometry or positions Tack has
written since.
---
 changelog.d/derived-fields.md                 |  7 ++--
 .../src/tack/dataset/algorithms/_common.py    | 32 +++++++++++++------
 .../dataset/algorithms/extract_surface.py     | 25 +++++++++------
 .../tack-vis/tests/test_derived_fields.py     | 17 ++++++++++
 4 files changed, 58 insertions(+), 23 deletions(-)

diff --git a/changelog.d/derived-fields.md b/changelog.d/derived-fields.md
index f96b7c158..00cae6ae7 100644
--- a/changelog.d/derived-fields.md
+++ b/changelog.d/derived-fields.md
@@ -2,9 +2,10 @@
 
 - Derived fields: a name a dataset stores no field under, at the location it is read,
   may stand for a field the dataset makes from itself (`register_derived_field(name,
-  compute, at)`), made the first time it is asked for and kept on the geometry.
-  `extract_surface` of the same geometry keeps its surface geometry too, so a
-  flux over many steps derives the normals once. One name is one quantity: `normals` is a surface's normals at its
+  compute, at)`), made the first time it is asked for and kept with the mesh: made
+  again when the mesh changes (another geometry, or positions Tack has written).
+  `extract_surface` of the same mesh keeps its surface geometry too, so a flux
+  over many steps derives the normals once. One name is one quantity: `normals` is a surface's normals at its
   points and cells and every mesh's face normals at its faces; `measures` each
   cell's volume, area or length and each face's area; `centroids` each cell's.
   Every field argument reads them (`named_field`), and so do field expressions,
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_common.py b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
index e7d41f6ef..8c5c76c80 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -147,10 +147,11 @@ def _same_space(*fields):
 def register_derived_field(name, compute, at):
     """Lets ``name`` stand, at ``at`` (``"points"``, ``"cells"`` or ``"faces"``), for a
     field a dataset makes from itself when it stores none of that name there:
-    ``compute(data)`` gives the ``Field`` on ``data``, a function of its geometry:
-    made the first time it is asked for and kept on the geometry, so every
-    dataset on the same geometry shares it (``extract_surface`` of the same
-    geometry gives the same surface geometry, and so the same normals). One name
+    ``compute(data)`` gives the ``Field`` on ``data``, a function of its mesh:
+    made the first time it is asked for and kept with the mesh (``mesh_kept``),
+    so every dataset on the same geometry shares it and it is made again when
+    the mesh changes. ``extract_surface`` of the same mesh gives the same surface
+    geometry, and so the same normals. One name
     is one quantity, registered at each location it has: Tack registers
     ``normals`` (a surface's at its points and cells, ``surface_normals``'; every
     mesh's at its faces, ``face_geometry``'s), ``measures`` (each cell's volume, area
@@ -178,13 +179,24 @@ def _location(field):
     return None
 
 
+def mesh_kept(data):
+    """A dict for what is computed from ``data``'s mesh alone (its derived fields,
+    its boundary's geometry), kept with the mesh: shared by every dataset on the
+    same geometry (and so the same topology), and emptied when Tack has written
+    the positions since (``arrays._revision``), so it changes when the mesh does."""
+    geometry = data.geometry
+    revision = arrays._revision(geometry.values)
+    kept = geometry.__dict__.get("_mesh_kept")
+    if kept is None or kept[0] != revision:
+        kept = geometry.__dict__["_mesh_kept"] = (revision, {})
+    return kept[1]
+
+
 def _derived(data, name, at):
-    # Part of the geometry, kept on it: any dataset on the same geometry (and so
-    # the same topology) shares them.
-    made = data.geometry.__dict__.setdefault("_derived", {})
-    if (name, at) not in made:
-        made[(name, at)] = _DERIVED[name][at](data)
-    return made[(name, at)]
+    kept = mesh_kept(data)
+    if (name, at) not in kept:
+        kept[(name, at)] = _DERIVED[name][at](data)
+    return kept[(name, at)]
 
 
 def named_field(data, name, at=None):
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py b/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py
index ab518cbfd..e0c630ae8 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py
@@ -209,15 +209,18 @@ def _fresh_boundary(data, fields):
 
 def _on_boundary(data, surface, fields, **maps):
     """The boundary ``surface`` with ``data``'s fields carried by ``maps``. Its
-    geometry is part of ``data``'s: kept on it, so a repeated call on the same
-    geometry gives the same surface geometry, and what is derived from that
-    (``normals``, ``measures``: ``register_derived_field``) is derived once."""
-    kept = data.geometry.__dict__.get("_boundary")
-    if kept is not None and kept[0] is surface:
-        out = carry(data, surface, geometry=kept[1], fields=fields, **kept[2])
+    geometry is part of ``data``'s mesh and kept with it (``mesh_kept``), so a
+    repeated call on the same mesh gives the same surface geometry, and what is
+    derived from that (``normals``, ``measures``) is derived once."""
+    from tack.dataset.algorithms._common import mesh_kept
+
+    kept = mesh_kept(data)
+    held = kept.get("boundary")
+    if held is not None and held[0] is surface:
+        out = carry(data, surface, geometry=held[1], fields=fields, **held[2])
     else:
         out = carry(data, surface, fields=fields, **maps)
-    data.geometry.__dict__["_boundary"] = (surface, out.geometry, maps)
+    kept["boundary"] = (surface, out.geometry, maps)
     return out
 
 
@@ -274,10 +277,12 @@ def extract_surface(data, faces=None, fields=None):
         fresh = _fresh_boundary(data, fields)
         if fresh is not None:
             return fresh
-        kept = data.geometry.__dict__.get("_boundary")
-        if (kept is not None and "faces_to_cells" in kept[2]
+        from tack.dataset.algorithms._common import mesh_kept
+
+        held = mesh_kept(data).get("boundary")
+        if (held is not None and "faces_to_cells" in held[2]
                 and not (data.topology.dimension == 3 and _has_order_two(data, fields))):
-            return _on_boundary(data, kept[0], fields, **kept[2])
+            return _on_boundary(data, held[0], fields, **held[2])
     boundary = faces is None
     if boundary:
         from tack.dataset.algorithms.boundary_faces import boundary_faces
diff --git a/packages/tack-vis/tests/test_derived_fields.py b/packages/tack-vis/tests/test_derived_fields.py
index 8bab79dc8..497cf1e18 100644
--- a/packages/tack-vis/tests/test_derived_fields.py
+++ b/packages/tack-vis/tests/test_derived_fields.py
@@ -59,6 +59,23 @@ def test_a_boundary_of_the_same_geometry_derives_its_normals_once(backend):
     assert td.named_field(moved, "normals", "cells") is not td.named_field(first, "normals", "cells")
 
 
+def test_mesh_metrics_change_when_the_mesh_does(backend):
+    tets = td.tetrahedralize(_grid())
+    solid = td.DataSet(tets.topology, tets.geometry.values.numpy())
+    before = td.named_field(solid, "measures", "cells")
+    surface = td.extract_surface(solid)
+    normals = td.named_field(surface, "normals", "cells")
+    stretched = solid.geometry.values.numpy() * np.float32(2.0)
+    tack.copyto(solid.geometry.values, stretched)
+    after = td.named_field(solid, "measures", "cells")
+    assert after is not before
+    np.testing.assert_allclose(after.values.numpy().sum(), 8 * before.values.numpy().sum(),
+                               rtol=1e-5)
+    again = td.extract_surface(solid)
+    assert again.geometry is not surface.geometry
+    assert td.named_field(again, "normals", "cells") is not normals
+
+
 def test_registered_names_and_misses(backend):
     calls = []