Vis/field resolver

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.

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.

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`).

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.

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.

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

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

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
 
 

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

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

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