From 4f7a52d6d3f015e7bfc293b65e619afb72d7585d Mon Sep 17 00:00:00 2001
From: Jeff Lee <jeff.lee@kitware.com>
Date: Thu, 5 Nov 2026 15:25:11 +0000
Subject: [PATCH] Fields evaluated where they are read, declared on Space

Space declares evaluated, continuous, default_subdivision, integration_order, rebind and stored, each with its default, and the algorithms read them directly. integrate takes its quadrature order from the space. A toy space in the tests, a function of position computed where it is read, checks contour, clip, threshold, carry, integrate, statistics, entropy, owned_values and the writers' stored values. threshold now reads such a field at its lattice points rather than as cell data.
---
 changelog.d/evaluated-fields.md               |  12 ++
 docs/vis-guide/reference/dataset-api.md       |  20 ++-
 .../src/tack/dataset/algorithms/_common.py    |   7 +
 .../src/tack/dataset/algorithms/contour.py    |   6 +-
 .../src/tack/dataset/algorithms/entropy.py    |   6 +-
 .../src/tack/dataset/algorithms/integrate.py  |  16 +-
 .../src/tack/dataset/algorithms/statistics.py |   6 +-
 .../src/tack/dataset/algorithms/subdivide.py  |  18 +-
 .../src/tack/dataset/algorithms/threshold.py  |   2 +-
 packages/tack-vis/src/tack/dataset/carry.py   |  13 +-
 packages/tack-vis/src/tack/dataset/dataset.py |  14 ++
 packages/tack-vis/src/tack/dataset/ghosts.py  |   4 +-
 packages/tack-vis/src/tack/dataset/spaces.py  |  41 +++++
 .../tack-vis/src/tack/interop/vtk/__init__.py |   5 +
 .../tack-vis/tests/test_evaluated_fields.py   | 168 ++++++++++++++++++
 15 files changed, 314 insertions(+), 24 deletions(-)
 create mode 100644 changelog.d/evaluated-fields.md
 create mode 100755 packages/tack-vis/tests/test_evaluated_fields.py

diff --git a/changelog.d/evaluated-fields.md b/changelog.d/evaluated-fields.md
new file mode 100644
index 000000000..494b3c7ff
--- /dev/null
+++ b/changelog.d/evaluated-fields.md
@@ -0,0 +1,12 @@
+### Added
+
+- Fields evaluated where they are read: a space may declare itself `evaluated`
+  (its view computes the field by `value` and `corner_value`; it stores no values),
+  `continuous` (its crossings are shared between cells, as `H1`'s are), its
+  `default_subdivision` and its `integration_order`, and give `rebind` and
+  `stored`. `contour`, `subdivide`, `clip` and `threshold` read such a field
+  through its view (`threshold` by its values at the lattice points, not as cell
+  data), filters make it again on their output, `integrate` integrates
+  it by quadrature to its order, and `statistics`, `entropy`, `owned_values` and
+  `dataset_to_vtk` read its stored values. `Space` declares each, with its
+  default.
diff --git a/docs/vis-guide/reference/dataset-api.md b/docs/vis-guide/reference/dataset-api.md
index df81a3496..d8d4f63c9 100644
--- a/docs/vis-guide/reference/dataset-api.md
+++ b/docs/vis-guide/reference/dataset-api.md
@@ -212,7 +212,25 @@ array accesses.
   area, an id) are left alone.
 
 Attributes: `topology`, `size`, `on` (the entity kind the values are
-indexed by), `interpolated`, and `order` where it applies.
+indexed by), `interpolated`, `continuous`, and `order` where it applies.
+
+### Fields evaluated where they are read
+
+A space may compute its field at each point an algorithm reads it instead of
+storing values, as a field expression does. `Space` declares what such a space
+gives, each with its default:
+
+| Attribute | Default | Meaning |
+|---|---|---|
+| `evaluated` | `False` | the view's `value(c, pc)` and `corner_value(c, j)` compute the field (with `value_range(c)`, a bound on the cell, which may be unbounded); its `values` hold nothing to read |
+| `continuous` | `False` (`True` for `H1`) | `contour` shares crossings between cells, so its surface is watertight |
+| `default_subdivision` | `1` | the lattice steps along each cell axis that `contour`, `clip` and `threshold` cut a cell into by default |
+| `integration_order` | `None` | the polynomial order `integrate`'s quadrature is exact to; `None` and `integrate` refuses the field |
+| `rebind(field, data)` | `None` | the field on a filter's output, made from the fields carried there; `None` and the filter leaves it out |
+| `stored(field, data)` | `field` | the field as stored values on `data`'s device, which `statistics`, `entropy`, `owned_values` and `dataset_to_vtk` read |
+
+An unbounded `value_range` cuts every cell, so `clip` and `threshold` give the
+same region as for a stored field of the same values, in more pieces.
 
 ### Moving a field between spaces
 
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/_common.py b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
index d50425db3..3ee8a3884 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/_common.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/_common.py
@@ -158,6 +158,13 @@ def _field(data, field):
     return field
 
 
+def _stored_field(data, field):
+    """``_field``, read as stored values (its space's ``stored``): for an algorithm
+    that reduces a field's values rather than reading it through its view."""
+    field = _field(data, field)
+    return field.space.stored(field, data)
+
+
 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/contour.py b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
index aa2ee8cec..d77167b3e 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/contour.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/contour.py
@@ -461,8 +461,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 one of a
+    space that declares itself ``continuous``), 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
@@ -520,7 +520,7 @@ 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
+    per_cell = 0 if space.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/entropy.py b/packages/tack-vis/src/tack/dataset/algorithms/entropy.py
index 85b82d999..0d1de2a81 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
@@ -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 0c7d014da..0c7e91f23 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/integrate.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/integrate.py
@@ -206,7 +206,12 @@ 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: an evaluated
+    space's ``integration_order``."""
+    if space.evaluated:
+        if space.integration_order is None:
+            raise TypeError(f"integrate needs {space!r} to give its integration_order")
+        return int(space.integration_order)
     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)
@@ -327,6 +332,8 @@ def _piece(data, names):
     for (name, _), field in zip(names, given):
         field = _field(data, field)
         field_units[name] = field.unit
+        if not getattr(data.topology, "reference_cells", True):
+            field = field.space.stored(field, data)
         if not getattr(data.topology, "reference_cells", True):
             space = _with_basis(data, field).space
             if _on_points(space):
@@ -366,7 +373,10 @@ def integrate(data, fields=None):
     its own map (``quadrature``: collapsed onto simplices and wedges), with enough
     points for the field's order and the geometry's, so a curved geometry and a
     field of any order (``H1``, ``L2``, point and cell data) integrate exactly, to
-    rounding. Polyhedra, polygons and polylines take their measure from
+    rounding. An evaluated field (an expression's) is integrated through its view
+    in the same way, with points for its space's ``integration_order``: exact where
+    the field is a polynomial of that order, a quadrature approximation where it
+    is not. Polyhedra, polygons and polylines take their measure from
     ``cell_geometry`` and integrate cell data and point data, which is linear on
     their pieces (sub-tetrahedra, triangles, segments), exactly. Only the cells of the mesh's dimension count, as in VTK: a volume's
     integrals leave out its surfaces, lines and vertices. Only owned cells count
@@ -408,7 +418,7 @@ def _names(data, fields):
     gives the field is resolved once a piece, where it is integrated."""
     if fields is None:
         return [(name, name) for name, field in data.fields.items()
-                if isinstance(field.space, (H1, L2, Constant))
+                if isinstance(field.space, (H1, L2, Constant)) or field.space.evaluated
                 or (isinstance(field.space, Values) and field.space.on in ("points", "cells"))]
     if hasattr(fields, "items"):
         names = list(fields.items())
diff --git a/packages/tack-vis/src/tack/dataset/algorithms/statistics.py b/packages/tack-vis/src/tack/dataset/algorithms/statistics.py
index c5f584ffa..20d09f9d2 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)
@@ -150,7 +150,7 @@ def _weights(data, field, weights):
         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..18790133b 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/subdivide.py
@@ -865,7 +865,7 @@ def _order(field):
     space = field.space
     if isinstance(space, (H1, L2)) and not space.varies:
         return space.order
-    return 1
+    return space.default_subdivision
 
 
 def default_level(data, *fields):
@@ -1029,7 +1029,7 @@ 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:
+        if not (isinstance(field.space, (H1, L2)) or field.space.evaluated) or field.space.varies:
             continue
         out = tack.empty((points,), _element(field.values))
         if count:
@@ -1370,13 +1370,15 @@ def level_for(data, field, flags, tolerance):
 
 
 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,
+    """Whether cutting ``field`` on ``data`` takes subdivision: a curved geometry, a
+    field above order 1, or an evaluated field whose ``default_subdivision`` is above
+    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 space.default_subdivision > 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/algorithms/threshold.py b/packages/tack-vis/src/tack/dataset/algorithms/threshold.py
index d562f9363..efa14b600 100644
--- a/packages/tack-vis/src/tack/dataset/algorithms/threshold.py
+++ b/packages/tack-vis/src/tack/dataset/algorithms/threshold.py
@@ -76,7 +76,7 @@ def threshold(data, field, lower, upper, all_points=True, fields=None, subdivisi
         raise TypeError(f"threshold needs a field on the points or cells, not {space!r}")
     from tack.dataset.algorithms.subdivide import split_by_range, subdivides
 
-    by_cells = isinstance(space, Constant) or (space.on == "cells"
+    by_cells = isinstance(space, Constant) or (space.on == "cells" and not space.evaluated
                                                and not isinstance(space, (H1, L2)))
     if not by_cells and getattr(data.topology, "reference_cells", True) and (
             subdivides(data, field) or tolerance is not None
diff --git a/packages/tack-vis/src/tack/dataset/carry.py b/packages/tack-vis/src/tack/dataset/carry.py
index f1eea220b..ddb687b9e 100644
--- a/packages/tack-vis/src/tack/dataset/carry.py
+++ b/packages/tack-vis/src/tack/dataset/carry.py
@@ -502,8 +502,12 @@ 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):
         source = data.fields[name]
+        if source.space.evaluated:
+            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 source.gradient is not None:
@@ -527,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)
+    # An evaluated field has no values to move: its space makes it again on the
+    # output, from the fields carried there.
+    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 a4e18e0d5..0d76b7514 100644
--- a/packages/tack-vis/src/tack/dataset/dataset.py
+++ b/packages/tack-vis/src/tack/dataset/dataset.py
@@ -737,6 +737,20 @@ def __repr__(self):
         return f"Fields({{{', '.join(f'{k!r}: {f.space!r}' for k, f in self.items())}}})"
 
 
+def _with_stored_values(data):
+    """``data`` with each evaluated field replaced by its stored values (its space's
+    ``stored``), for writers, which take arrays; ``data`` itself when it has none."""
+    entries = data.fields.entries()
+    if not any(f.space.evaluated for _, f in entries):
+        return data
+    import copy
+
+    out = copy.copy(data)
+    out.fields = [(name, f.space.stored(f, data) if f.space.evaluated 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/dataset/ghosts.py b/packages/tack-vis/src/tack/dataset/ghosts.py
index 913d4c198..3add96987 100644
--- a/packages/tack-vis/src/tack/dataset/ghosts.py
+++ b/packages/tack-vis/src/tack/dataset/ghosts.py
@@ -223,7 +223,8 @@ 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``, or its name; an evaluated one's
+    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
@@ -234,6 +235,7 @@ def owned_values(data, field):
 
     if is_field_key(field):
         field = data.fields[field]
+    field = field.space.stored(field, data)
     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/spaces.py b/packages/tack-vis/src/tack/dataset/spaces.py
index 39c8e747a..e20e05d6b 100644
--- a/packages/tack-vis/src/tack/dataset/spaces.py
+++ b/packages/tack-vis/src/tack/dataset/spaces.py
@@ -70,6 +70,45 @@ class Space:
     #: Whether cells differ in how their values are read (``cell_keys``).
     varies = False
 
+    #: Whether a field of this space is continuous across cells, so ``contour``
+    #: shares its crossings between them and its surface is watertight (``H1``).
+    continuous = False
+
+    #: Whether a field of this space is evaluated where it is read, not stored: an
+    #: interpolated space whose view computes the field at the points an algorithm
+    #: asks for, by ``value(c, pc)`` and ``corner_value(c, j)`` (and ``value_range(c)``,
+    #: a bound on the cell's values, which may be unbounded). Its ``values`` hold
+    #: nothing to read. Algorithms that read values rather than the view (the
+    #: reductions, ``owned_values``, writers) read ``stored``; filters make it again on
+    #: their output with ``rebind``; ``subdivide`` samples it at each lattice point.
+    evaluated = False
+
+    #: For an evaluated space: the lattice steps along each cell axis that
+    #: ``contour``, ``clip`` and ``threshold`` cut a cell into by default to sample it
+    #: (``default_level``), as an ``H1`` field's order does; 1 samples the cell's
+    #: corners only.
+    default_subdivision = 1
+
+    #: For an evaluated space: the polynomial order along each axis ``integrate``
+    #: integrates it exactly to. Its quadrature takes enough points for that order
+    #: and the geometry's, so the integral is exact where the field is a polynomial
+    #: of at most that order in the reference coordinates, and a quadrature
+    #: approximation where it is not. ``None``: ``integrate`` refuses the field.
+    integration_order = None
+
+    def rebind(self, field, data):
+        """For an evaluated space: ``field`` on ``data``, a filter's output, made
+        again from what the filter carried there, or ``None`` when it cannot be (a
+        field it reads was not carried); the filter then leaves it out. Other
+        spaces' fields are carried by their values and never asked."""
+
+    def stored(self, field, data):
+        """``field`` (of this space, on ``data``) as a field of stored values on
+        ``data``'s device, for an algorithm that reads values rather than the view:
+        ``field`` itself for a space that stores its values; an evaluated space gives
+        its values at the points or cells, in a space that stores them."""
+        return field
+
     def __new__(cls, where, *args, **kwargs):
         topology = getattr(where, "topology", where)
         params = cls._params(*args, **cls._defaults(topology, args, kwargs))
@@ -445,6 +484,8 @@ class H1(Space):
     is its point's), and the space holds as many values as the largest index and
     the points need. A filter's output numbers its values by entity again."""
 
+    continuous = True
+
     on = "points"
     interpolated = True
     basis = None
diff --git a/packages/tack-vis/src/tack/interop/vtk/__init__.py b/packages/tack-vis/src/tack/interop/vtk/__init__.py
index 06aa96a53..b1a167926 100644
--- a/packages/tack-vis/src/tack/interop/vtk/__init__.py
+++ b/packages/tack-vis/src/tack/interop/vtk/__init__.py
@@ -2264,8 +2264,13 @@ 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 _with_stored_values
     from tack.interop.vtk import _lagrange
 
+    if isinstance(data, DataSet):
+        data = _with_stored_values(data)
+
     if cell_grid:
         from tack.interop.vtk import _cellgrid
 
diff --git a/packages/tack-vis/tests/test_evaluated_fields.py b/packages/tack-vis/tests/test_evaluated_fields.py
new file mode 100755
index 000000000..01bd01e0a
--- /dev/null
+++ b/packages/tack-vis/tests/test_evaluated_fields.py
@@ -0,0 +1,168 @@
+# SPDX-FileCopyrightText: Copyright (c) Kitware Inc.
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""Fields evaluated where they are read, through a toy space: f = x y + z computed at
+each point an algorithm asks for, with nothing stored. Each algorithm that reads such
+a field is checked against the same function stored in H1 of order 2, which holds
+it exactly on these hexahedra."""
+
+import itertools
+
+import numpy as np
+import pytest
+
+import tack
+import tack.dataset as td
+from tack.dataset.spaces import Space
+
+_keys = itertools.count(1)
+
+
+@tack.func
+def _f(x):
+    return x[0] * x[1] + x[2]
+
+
+class _View:
+    @tack.func
+    def value(self, c, pc):
+        return _f(self.dom.position(c, pc))
+
+    @tack.func
+    def corner_value(self, c, j):
+        return self.value(c, self.dom.parametric_point(j))
+
+    @tack.func
+    def value_range(self, c):
+        return -3.0e38, 3.0e38
+
+    @tack.func
+    def dof(self, c, k):
+        return 0.0 * self.dom.parametric_center()[0]
+
+
+class Analytic(Space):
+    """``_f``, evaluated where it is read. ``carried`` says whether a filter's output
+    can make it again (``rebind``)."""
+
+    on = "cells"
+    interpolated = True
+    evaluated = True
+    continuous = True
+    default_subdivision = 2
+    integration_order = 2
+
+    @classmethod
+    def _params(cls, key=0):
+        return {"key": key}
+
+    @property
+    def size(self):
+        return self.topology.num_cells
+
+    def mixin_for(self, key, shape=None):
+        return _View
+
+    def attributes(self, group=None, dtype=None):
+        return {} if group is None else {"dom": self.data.domain_view("cells", group)}
+
+    def rebind(self, field, data):
+        return analytic(data, self.carried) if self.carried else None
+
+    def stored(self, field, data):
+        return td.interpolate(data, _f, td.H1(data))
+
+
+def analytic(data, carried=True):
+    space = Analytic(data, key=next(_keys))
+    space.data, space.carried = data, carried
+    return td.Field(space, tack.zeros((data.num_cells,), tack.f32))
+
+
+def _box():
+    data = td.uniform_grid((5, 4, 4), spacing=(0.25, 1 / 3, 1 / 3))
+    data.fields["f"] = analytic(data)
+    data.fields["exact"] = td.interpolate(data, _f, td.H1(data, 2))
+    return data
+
+
+def test_space_declares_the_contract_with_its_defaults():
+    from tack.dataset.spaces import H1, Constant
+    data = td.uniform_grid((2, 2, 2))
+    stored = td.Field(Constant(data), tack.zeros((data.num_cells,), tack.f32))
+    for space in (H1(data), Constant(data)):
+        assert not space.evaluated and space.default_subdivision == 1
+        assert space.integration_order is None
+        assert space.rebind(stored, data) is None
+    assert H1(data).continuous and not Constant(data).continuous
+    assert Constant(data).stored(stored, data) is stored
+
+
+def test_contour_reads_it_through_its_view(backend):
+    data = _box()
+    mine = td.contour(data, "f", 0.6)
+    exact = td.contour(data, "exact", 0.6)
+    assert mine.num_cells == exact.num_cells
+    # Continuous: crossings shared between cells, as H1's are.
+    assert mine.num_points == exact.num_points
+    assert td.integrate(mine).measure == pytest.approx(td.integrate(exact).measure, rel=1e-5)
+
+
+def _cells(out):
+    return sum(piece.num_cells for piece in out) if isinstance(
+        out, td.PartitionedDataSet) else out.num_cells
+
+
+def test_clip_and_threshold_subdivide_by_its_default(backend):
+    data = _box()
+    # Its value_range is unbounded, so every cell is cut where H1's bounds keep some
+    # whole: the same region, in more pieces.
+    mine, exact = (td.clip(data, f, 0.6, exact=False) for f in ("f", "exact"))
+    assert _cells(mine) >= _cells(exact)
+    assert td.integrate(mine).measure == pytest.approx(td.integrate(exact).measure, rel=1e-5)
+    mine, exact = (td.threshold(data, f, 0.0, 0.6) for f in ("f", "exact"))
+    assert isinstance(mine, td.PartitionedDataSet)
+    assert td.integrate(mine).measure == pytest.approx(td.integrate(exact).measure, rel=1e-5)
+
+
+def test_filters_make_it_again_on_their_output(backend):
+    data = _box()
+    data.fields["dropped"] = analytic(data, carried=False)
+    x = td.cell_centers(data).values.numpy()[:, 0].astype(np.float32)
+    data.fields["x"] = td.Field(td.Constant(data), tack.array(x))
+    out = td.threshold(data, "x", 0.0, 0.6)
+    assert out.fields["f"].space.evaluated
+    assert out.fields["f"].space.data is out
+    assert "dropped" not in out.fields
+
+
+def test_integrate_to_its_order(backend):
+    data = _box()
+    # Over [0, 1]^3: the integral of x y + z is 1/4 + 1/2.
+    assert td.integrate(data, "f").values["f"] == pytest.approx(0.75, rel=1e-6)
+    data.fields["g"] = analytic(data)
+    data.fields["g"].space.integration_order = None
+    with pytest.raises(TypeError, match="integration_order"):
+        td.integrate(data, "g")
+
+
+def test_reductions_read_its_stored_values(backend):
+    data = _box()
+    points = td.interpolate(data, _f, td.H1(data))
+    by_values = td.statistics(data, points)
+    mine = td.statistics(data, "f")
+    assert mine.n == data.num_points
+    assert mine.mean == pytest.approx(by_values.mean)
+    assert td.entropy(data, "f") == pytest.approx(td.entropy(data, points))
+    np.testing.assert_allclose(td.owned_values(data, "f").numpy(), points.values.numpy())
+
+
+def test_writers_take_its_stored_values():
+    from tack.dataset.dataset import _with_stored_values
+
+    data = _box()
+    written = _with_stored_values(data)
+    assert not written.fields["f"].space.evaluated
+    assert written.fields["exact"] is data.fields["exact"]
+    plain = td.uniform_grid((2, 2, 2))
+    assert _with_stored_values(plain) is plain