Fields evaluated where they are read, declared on Space

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.

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
 

packages/tack-vis/src/tack/dataset/algorithms/_common.py

@@ -175,6 +175,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

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

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

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())

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\"")

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)

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

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:
@@ -529,5 +533,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

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

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")

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

packages/tack-vis/src/tack/interop/vtk/__init__.py

@@ -2265,8 +2265,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