VTK field algorithms write every node of VTK's higher-order cells

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

@@ -1026,7 +1026,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 +1066,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."""
@@ -2600,9 +2650,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
 

packages/tack-vis/tests/test_vtk_algorithms.py

@@ -974,6 +974,95 @@ def test_every_algorithm_has_a_wrapper():
     assert not missing
 
 
+# ── Higher-order cells: fields at VTK's nodes ───────────────────────
+
+def _quadratic_quads():
+    """Two quad8 cells (VTK_QUADRATIC_QUAD) side by side, the shared edge's middle
+    node pushed off the line: VTK's points are corners and edge middles, 13."""
+    corners = [(0, 0), (1, 0), (2, 0), (0, 1), (1, 1), (2, 1)]
+    middles = [(0.5, 0), (1, 0.5 + 0.0), (0.5, 1), (0, 0.5), (1.5, 0), (2, 0.5), (1.5, 1)]
+    points = np.array([(x, y, 0.0) for x, y in corners + middles])
+    points[7, 0] = 1.1                                  # the shared edge, curved
+    grid = vtk.vtkUnstructuredGrid()
+    grid.points = vtk.vtkPoints(data=points)
+    grid.Allocate(2)
+    grid.InsertNextCell(vtk.VTK_QUADRATIC_QUAD, 8, [0, 1, 4, 3, 6, 7, 8, 9])
+    grid.InsertNextCell(vtk.VTK_QUADRATIC_QUAD, 8, [1, 2, 5, 4, 10, 11, 12, 7])
+    return grid
+
+
+def _quadratic_tetra():
+    corners = [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)]
+    edges = [(0, 1), (1, 2), (0, 2), (0, 3), (1, 3), (2, 3)]        # VTK's tet10 order
+    points = np.array(corners + [tuple((np.add(corners[a], corners[b])) / 2)
+                                 for a, b in edges], float)
+    points[4] += (0.0, -0.1, 0.0)                                    # one edge curved
+    grid = vtk.vtkUnstructuredGrid()
+    grid.points = vtk.vtkPoints(data=points)
+    grid.Allocate(1)
+    grid.InsertNextCell(vtk.VTK_QUADRATIC_TETRA, 10, list(range(10)))
+    return grid
+
+
+@pytest.mark.parametrize("make", [_quadratic_quads, _quadratic_tetra])
+def test_a_field_filter_writes_every_node_of_higher_order_cells(make):
+    """The new arrays have a value at each of VTK's points, the edge middles too, in
+    VTK's order: the elevation of each node's own position, and its gradient
+    constant. Tack numbers the nodes of its order-2 space otherwise, and keeps
+    values on points at the corners only."""
+    grid = make()
+    elevation = ta.Elevation(low_point=(0, 0, 0), high_point=(2, 0, 0), name="e")(grid)
+    assert elevation.GetNumberOfPoints() == grid.GetNumberOfPoints()
+    x = np.asarray(grid.points.data)
+    np.testing.assert_allclose(np.asarray(elevation.point_data["e"]),
+                               np.clip(x[:, 0] / 2, 0, 1), atol=1e-6)
+    gradient = ta.Gradients(field=("e", "points"), at="points", name="g")(elevation)
+    values = np.asarray(gradient.point_data["g"])
+    assert values.shape == (grid.GetNumberOfPoints(), 3)
+    np.testing.assert_allclose(values, np.tile([0.5, 0, 0], (len(values), 1)), atol=1e-5)
+
+
+def test_higher_order_nodes_on_another_device_and_integers_refused():
+    """A copy of the dataset on another device finds the nodes through the dataset kept
+    for the VTK one; an integer field on the points has no value at an edge's
+    middle, and is refused."""
+    from tack.runtime import dispatch
+    from tack.runtime.cpu import CPUBackend
+
+    class SecondCPU(CPUBackend):
+        display_name = "second device"
+
+        @property
+        def memory_domain(self):
+            return ("cpu", 1)
+
+    grid = _quadratic_quads()
+    other = dispatch.as_device(SecondCPU())
+    moved = ta.Elevation(low_point=(0, 0, 0), high_point=(2, 0, 0), name="e", device=other)
+    out = moved(grid)
+    assert moved.error is None, moved.error
+    np.testing.assert_allclose(np.asarray(out.point_data["e"]),
+                               np.asarray(grid.points.data)[:, 0] / 2, atol=1e-6)
+    ids = ta.PointIds()
+    ids.SetInputDataObject(0, grid)
+    assert ids.Update() == 0 and "integers" in str(ids.error)
+
+
+def test_temporal_statistics_on_higher_order_cells():
+    statistics = ta.RunningStatistics(fields=["e"])
+    for t in (0.0, 1.0):
+        step = ta.Elevation(low_point=(0, 0, 0), high_point=(2 + t, 0, 0), name="e")(
+            _quadratic_quads())
+        step.GetInformation().Set(vtk.vtkDataObject.DATA_TIME_STEP(), t)
+        statistics.SetInputDataObject(0, step)
+        statistics.Update()
+        assert statistics.error is None, statistics.error
+    x = np.asarray(_quadratic_quads().points.data)[:, 0]
+    np.testing.assert_allclose(
+        np.asarray(statistics.GetOutputDataObject(0).point_data["e_mean"]),
+        (x / 2 + x / 3) / 2, atol=1e-6)
+
+
 # ── Seeds on a port ─────────────────────────────────────────────────
 
 @pytest.mark.parametrize("make", [