diff --git a/changelog.d/derived-fields.md b/changelog.d/derived-fields.md
new file mode 100644
index 000000000..00cae6ae7
--- /dev/null
+++ b/changelog.d/derived-fields.md
@@ -0,0 +1,12 @@
+### 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 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,
+  so `integrate(surface, ["dot(U, normals)"])` is a flux with no filter first.
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/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/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/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/docs/vis-guide/reference/data-model.md b/docs/vis-guide/reference/data-model.md
index 4a29b6c13..ac22b254a 100644
--- a/docs/vis-guide/reference/data-model.md
+++ b/docs/vis-guide/reference/data-model.md
@@ -205,6 +205,14 @@ do, so pieces sharing an edge compute it alike.
 
 ::: tack.dataset.edge_crossing
 
+::: 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 9e9e1b36e..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
+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,8 @@
     "probe",
     "project",
     "rectilinear_grid",
+    "register_derived_field",
+    "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 5de4ccd2c..8c5c76c80 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -139,18 +139,137 @@ def _same_space(*fields):
     return space, n
 
 
+_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``, 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
+    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 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):
+    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):
+    """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):
+    """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 or text in _DERIVED:
+        return named_field(data, 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
 
 
+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
@@ -160,9 +279,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/contour.py b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
index aa2ee8cec..c1bf9009b 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/contour.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
@@ -461,14 +461,15 @@ 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
     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)
@@ -520,7 +521,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/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/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/extract_surface.py b/packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py
index 1fcc50b82..e0c630ae8 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,24 @@ 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 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)
+    kept["boundary"] = (surface, out.geometry, maps)
+    return out
 
 
 @on_inputs_device
@@ -260,6 +277,14 @@ def extract_surface(data, faces=None, fields=None):
         fresh = _fresh_boundary(data, fields)
         if fresh is not None:
             return fresh
+        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, held[0], fields, **held[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 +320,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/src/tack/dataset/algorithms/gradients.py b/packages/tack-vis/src/tack/dataset/algorithms/gradients.py
index 39fe44868..6d713e6e4 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/gradients.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/gradients.py
@@ -20,20 +20,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
@@ -49,36 +62,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
@@ -278,6 +312,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]):
@@ -355,6 +407,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``
@@ -365,6 +420,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 "
@@ -394,7 +452,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/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/algorithms/subdivide.py b/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
index e5eb57a50..a67cac2b0 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
@@ -865,7 +865,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):
@@ -1029,7 +1030,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:
@@ -1373,10 +1377,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 29545fb6b..1a3ffac2c 100644
--- a/packages/tack-vis/src/tack/dataset/carry.py
+++ b/packages/tack-vis/src/tack/dataset/carry.py
@@ -502,9 +502,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:
@@ -521,5 +531,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 b01f7b307..279f6d650 100644
--- a/packages/tack-vis/src/tack/dataset/dataset.py
+++ b/packages/tack-vis/src/tack/dataset/dataset.py
@@ -54,9 +54,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
@@ -67,6 +74,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):
@@ -96,7 +113,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):
@@ -133,6 +151,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))
@@ -142,6 +162,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
@@ -696,6 +736,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
@@ -781,7 +836,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/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/src/tack/interop/vtk/__init__.py b/packages/tack-vis/src/tack/interop/vtk/__init__.py
index e84fcd433..83e30c2ba 100644
--- a/packages/tack-vis/src/tack/interop/vtk/__init__.py
+++ b/packages/tack-vis/src/tack/interop/vtk/__init__.py
@@ -2247,8 +2247,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/_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
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..497cf1e18
--- /dev/null
+++ b/packages/tack-vis/tests/test_derived_fields.py
@@ -0,0 +1,98 @@
+# 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_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_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 = []
+
+    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"]
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..698f18090
--- /dev/null
+++ b/packages/tack-vis/tests/test_field_resolver.py
@@ -0,0 +1,65 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""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)
+
+
+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"])
+
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..4134e5d18
--- /dev/null
+++ b/packages/tack-vis/tests/test_given_gradients.py
@@ -0,0 +1,132 @@
+# 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."""
+
+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)
+
+
+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])
+