VTK algorithms keep the input's active attributes, and set their own

changelog.d/vtk-algorithms-active-attributes.md

@@ -0,0 +1,11 @@
+### Fixed
+
+- The VTK algorithms' outputs had no active attributes: Tack keeps no roles, so the
+  arrays it carried came back plain, and mappers colored by nothing, glyphs and
+  lighting found no vectors or normals (VTK's tests rendered white). An output now
+  keeps its input's active scalars, vectors, normals, texture coordinates and
+  tensors where it carries arrays of those names (`Probe` its source's), as VTK's
+  filters pass them; `Elevation` and `ImplicitValues` make their array the active
+  scalars, `SurfaceNormals`, `PointNormals` and `Tubes` theirs the active normals
+  (a tube took the lines' carried `Normals` before, and was lit wrongly). Found
+  running VTK's Python tests on Tack.

changelog.d/vtk-field-wrappers-higher-order.md

@@ -0,0 +1,13 @@
+### Fixed
+
+- The VTK algorithms that add fields (`Elevation`, `Gradients`, `Project`, ... and
+  `RunningStatistics`, `TimeDerivative`) on VTK's higher-order cells wrote Tack's
+  values straight into VTK's point data: an order-2 field's in Tack's node order,
+  not VTK's (`Elevation` came out scrambled, and its `Gradients` wrong), and values
+  on Tack's points, the corners, for every VTK point, edge middles included (a
+  1193-point grid got an array of 1193 tuples of 1227 components). Each VTK point
+  now takes its node's value: the field put into the space the import made
+  (`project`, exact for an `H1` field of that order or below and for values on the
+  corners), its nodes matched to VTK's points as the import matched them. An
+  integer field on the points is refused there, having no value at an edge's
+  middle. Found running VTK's `quadraticQuadDeriv.py` on Tack.

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

@@ -223,6 +223,11 @@ def _node_places(cells, u, node_of, starts, out):
             out[starts[e] + n] = u.dof_index(c, node_of[n])
 
 
+#: The attribute of a topology made from VTK's higher-order cells that says where each
+#: VTK point is in the import's ``H1`` space: ``(order, basis, node of each point)``.
+VTK_NODES = "_vtk_nodes"
+
+
 def to_dataset(grid, types, dtype, id_dtype, shift, unit_of, float_fields):
     """A ``vtkUnstructuredGrid`` of VTK's higher-order cells -- Lagrange, Bezier, the
     classic quadratic, biquadratic and triquadratic ones, cubic lines -- of any
@@ -295,9 +300,13 @@ def to_dataset(grid, types, dtype, id_dtype, shift, unit_of, float_fields):
     point_arrays = [name for name in grid.point_data.keys() if name not in skipped]
     sources = [points] + [np.asarray(grid.point_data[name]).reshape(len(points), -1)
                           for name in point_arrays]
+    # Which node of the space each VTK point is (-1 for a point no cell uses): the
+    # field algorithms write a field back onto these points through it (_VTK_NODES).
+    node_of_point = np.full(len(points), -1, np.int64)
     if order == 1:
         space = H1(data)
         node_values = [source[corner_ids] for source in sources]
+        node_of_point[corner_ids] = np.arange(len(corner_ids))
     else:
         space = H1(data, order=order, basis=basis)
         # Where each cell's nodes live in the space, cell by cell.
@@ -330,22 +339,25 @@ def to_dataset(grid, types, dtype, id_dtype, shift, unit_of, float_fields):
                     out[at.reshape(-1)] = source[vtk_points.reshape(-1)]
                 seen = np.full(space.size, -1, np.int64)
                 seen[at.reshape(-1)] = vtk_points.reshape(-1)
+                node_of_point[vtk_points.reshape(-1)] = at.reshape(-1)
                 shared_bad |= bool((seen[at.reshape(-1)] != vtk_points.reshape(-1)).any())
             else:
                 if bezier:
                     raise NotImplementedError("Bezier cells whose nodes are not Tack's")
+                own = index >= 0
                 for out, source in zip(node_values, sources, strict=True):
                     values = np.einsum("nj,cjk->cnk", weights, source[rows])
                     # The cell's own nodes copied, as above; the others interpolated.
-                    own = index >= 0
                     values[:, own] = source[rows[:, index[own]]]
                     fresh = ~written[at]
                     out[at[fresh]] = values[fresh]
+                node_of_point[rows[:, index[own]].reshape(-1)] = at[:, own].reshape(-1)
             written[at.reshape(-1)] = True
         if shared_bad:
             raise NotImplementedError("the higher-order cells do not share their points as an "
                                       "H1 space shares its nodes (each cell's own copies)")
         data.geometry = Field(space, _vectors(node_values[0] - shift, dtype))
+    setattr(topology, VTK_NODES, (order, basis if order > 1 else None, node_of_point))
     for name, values in zip(point_arrays, node_values[1:], strict=True):
         array = grid.point_data[name]
         if not np.issubdtype(np.asarray(array).dtype, np.floating) and order > 1:

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

@@ -830,8 +830,18 @@ def RequestData(self, request, inInfo, outInfo):
             output.CompositeShallowCopy(self._tree(input))
         else:
             self._set_unit_output(output, self._unit(input))
+        _keep_active(self.GetInputData(inInfo, self.ACTIVE_FROM, 0), output,
+                     self._made_active())
         return 1
 
+    #: The input port whose active attributes the output keeps (Probe's source's).
+    ACTIVE_FROM = 0
+
+    def _made_active(self):
+        """``{(where, role): name}``: the arrays this filter makes that VTK's makes
+        active (``Scalars``, ``Normals``, ...), set over what the input's give."""
+        return {}
+
     def _by_box(self, vtk_amr):
         """The output for an AMR hierarchy this filter takes box by box: each level's
         boxes, this rank's, as the pieces of a partitioned dataset (``_level``), the
@@ -925,6 +935,56 @@ def compute(self, data):
         return self._function()(data, **self.arguments())
 
 
+#: VTK's attribute roles, as vtkDataSetAttributes names them.
+_ROLES = ("Scalars", "Vectors", "Normals", "TCoords", "Tensors")
+
+
+def _datasets_in(data):
+    """The datasets in ``data``: itself, or a composite's leaves."""
+    if data is None:
+        return []
+    if isinstance(data, vtkDataObject) and data.IsA("vtkCompositeDataSet"):
+        out, iterator = [], data.NewIterator()
+        iterator.InitTraversal()
+        while not iterator.IsDoneWithTraversal():
+            out.append(iterator.GetCurrentDataObject())
+            iterator.GoToNextItem()
+        return [d for d in out if d.IsA("vtkDataSet")]
+    return [data] if data.IsA("vtkDataSet") else []
+
+
+def _keep_active(source, out, made):
+    """Give every dataset in ``out`` the active attributes VTK's filter would: those of
+    ``source`` (its first dataset's) where an array of the name is there and the
+    role has none yet -- the scalars a contour carried, the vectors a clip kept --
+    then ``made``, the filter's own (``_made_active``). Tack keeps no roles, so the
+    arrays it carries come back plain, and mappers, glyphs and lighting read the
+    roles."""
+    names = {}
+    for leaf in _datasets_in(source)[:1]:
+        for where, attributes in (("points", leaf.GetPointData()), ("cells", leaf.GetCellData())):
+            for role in _ROLES:
+                array = getattr(attributes, f"Get{role}")()
+                if array is not None and array.GetName():
+                    names[where, role] = array.GetName()
+    if not names and not made:
+        return
+    for leaf in _datasets_in(out):
+        # Setting a role modifies the dataset: an entry valid before stays valid,
+        # the arrays being the same.
+        entry = _entry(leaf)
+        attributes = {"points": leaf.GetPointData(), "cells": leaf.GetCellData()}
+        for (where, role), name in names.items():
+            target = attributes[where]
+            if getattr(target, f"Get{role}")() is None and target.GetAbstractArray(name):
+                getattr(target, f"SetActive{role}")(name)
+        for (where, role), name in made.items():
+            if attributes[where].GetAbstractArray(name):
+                getattr(attributes[where], f"SetActive{role}")(name)
+        if entry is not None:
+            setattr(_holder(leaf), _CACHE, (_key(leaf), entry[1]))
+
+
 def _without_tracebacks(error, seen=None):
     """``error`` with its traceback, and those of the exceptions it was raised from or
     during, dropped."""
@@ -1026,7 +1086,9 @@ def _with_fields(self, block, piece):
                                       or (isinstance(space, Values) and space.on == "cells")):
                 raise TypeError(f"{name!r} is on {space!r}, which VTK keeps nowhere")
             attributes = out.point_data if on_points else out.cell_data
-            _set_values(attributes, name, _to_vtk_values(field.values, share=True, name=name))
+            values, share = _vtk_values_of(block, piece, field, name) if on_points else \
+                (field.values, True)
+            _set_values(attributes, name, _to_vtk_values(values, share=share, name=name))
             _label(attributes[name], field.unit)
         if any(name in piece.fields for name in fields):
             return out, None
@@ -1064,6 +1126,54 @@ def _output_name(self):
         raise ValueError(f"{type(self).__name__} needs a name for its new field")
 
 
+def _vtk_values_of(block, piece, field, name):
+    """``(values, share)``: a point field of ``piece``, the Tack dataset of the VTK
+    dataset ``block``, as values at ``block``'s points. On linear cells those are
+    the field's own values, shared. On VTK's higher-order cells, whose points are
+    the cells' nodes, the field goes into the space the import made (``project``:
+    exactly, for an ``H1`` field of that order or below and for values on Tack's
+    points, the corners) and each VTK point takes its node's value, a copy; the
+    import recorded which node each point is (``_lagrange.VTK_NODES``, on the
+    topology, or on the topology of the dataset kept for ``block`` where ``piece``
+    is a copy on another device)."""
+    import numpy as np
+
+    from tack import arrays
+    from tack.dataset import H1, project
+    from tack.interop.vtk._lagrange import VTK_NODES
+
+    nodes = getattr(piece.topology, VTK_NODES, None)
+    if nodes is None:
+        entry = _entry(block)
+        if entry is not None:
+            nodes = getattr(entry[1][None].topology, VTK_NODES, None)
+    count = block.GetNumberOfPoints()
+    if nodes is None:
+        if count != piece.num_points:
+            raise ValueError(f"{name!r} has values at {piece.num_points} points, and the VTK "
+                             f"dataset {count}")
+        return field.values, True
+    order, basis, node_of = nodes
+    if len(node_of) != count:
+        raise ValueError(f"{name!r}: the VTK dataset has {count} points, its cells' nodes "
+                         f"{len(node_of)}")
+    space = H1(piece, order=order, basis=basis) if order > 1 else H1(piece)
+    if not (isinstance(field.space, H1) and field.space.order == space.order
+            and field.space.basis == space.basis):
+        kind = arrays.dtype_of(field.values)
+        if getattr(kind, "scalar_dtype", kind).numpy_dtype.kind != "f":
+            raise TypeError(f"{name!r} holds integers at the cells' corners, and VTK's "
+                            "higher-order cells have nodes between them, which an integer "
+                            "field has no value at")
+        field = project(piece, field, space)
+    values = np.asarray(arrays.to_host(field.values))
+    out = values[np.maximum(node_of, 0)]
+    out[node_of < 0] = 0                    # a point no cell uses
+    import tack
+
+    return tack.array(np.ascontiguousarray(out)), False
+
+
 class _Whole(TackAlgorithm):
     """A filter whose Tack function takes a ``PartitionedDataSet`` itself: a partitioned
     input goes to it whole, not piece by piece."""
@@ -1148,6 +1258,9 @@ class SurfaceNormals(TackAlgorithm):
 
     consistent = Parameter(True, "wind the polygons consistently first")
 
+    def _made_active(self):
+        return {("points", "Normals"): "normals"}
+
 
 class ConnectedComponents(_Whole):
     """``td.connected_components``: ``region`` on the cells and ``point_region`` on the
@@ -1234,6 +1347,7 @@ class Probe(TackAlgorithm):
     OUTPUT = "structure"
     INPUTS = 2
     AMR_PORTS = (1,)
+    ACTIVE_FROM = 1                     # the probed fields are the source's
     FUNCTION = "probe"
     EXCLUDE = ("where",)
 
@@ -1306,6 +1420,9 @@ class Elevation(_FieldAlgorithm):
     range = Parameter((0.0, 1.0), "the elevation's range")
     name = Parameter("elevation", "the new array's name")
 
+    def _made_active(self):
+        return {("points", "Scalars"): self._output_name()}
+
 
 class CellQuality(_FieldAlgorithm):
     """``td.cell_quality``: each measure (Verdict's, generalized) as a cell array of its
@@ -1428,6 +1545,9 @@ class Tubes(TackAlgorithm):
     FUNCTION = "tubes"
     FIELDS = ("radius_by",)
 
+    def _made_active(self):
+        return {("points", "Normals"): "normals"}
+
 
 class Glyphs(TackAlgorithm):
     """``td.glyphs``: an arrow at each point along the ``vectors`` field, scaled by
@@ -1610,6 +1730,9 @@ class ImplicitValues(_FieldAlgorithm):
 
     name = Parameter("implicit", "the new point array's name")
 
+    def _made_active(self):
+        return {("points", "Scalars"): self._output_name()}
+
 
 class CellIds(_FieldAlgorithm):
     """``td.cell_ids``: each cell's id, as a cell array ``name``. VTK's
@@ -1843,6 +1966,9 @@ class PointNormals(_FieldAlgorithm):
     FUNCTION = "point_normals"
     RESULT_FIELDS = ("normals",)
 
+    def _made_active(self):
+        return {("points", "Normals"): "normals"}
+
 
 class Outliers(_Whole):
     """``td.outliers``: the points with ``outlier`` (1 for an outlier) and
@@ -2600,9 +2726,11 @@ def onto(block, piece):
                     attributes.AddArray(ghosts)
         for name in names:
             field = piece.fields[name]
-            attributes = out.point_data if _on_points(field.space) else out.cell_data
-            _set_values(attributes, name, _to_vtk_values(field.values, share=True,
-                                                         name=name))
+            on_points = _on_points(field.space)
+            attributes = out.point_data if on_points else out.cell_data
+            values, share = _vtk_values_of(block, piece, field, name) if on_points else \
+                (field.values, True)
+            _set_values(attributes, name, _to_vtk_values(values, share=share, name=name))
             _label(attributes[name], field.unit)
         return out