From 80ea2cab8e483a8b06ae56778d32abb180209fad Mon Sep 17 00:00:00 2001 From: Berk Geveci Date: Wed, 4 Nov 2026 20:15:47 -0500 Subject: [PATCH] VTK field algorithms write every node of VTK's higher-order cells The algorithms that add fields wrote Tack's values straight into VTK's point data. On VTK's higher-order cells that put an order-2 field's values in Tack's node order, not VTK's (Elevation scrambled, its Gradients wrong), and values on Tack's points, the corners, where VTK's points include the edge middles (a 1193-point grid got an array of 1193 tuples of 1227 components). The import now records which node of its space each VTK point is, and the writers put a point field into that space (project) and give each VTK point its node's value. Found running VTK's quadraticQuadDeriv.py on Tack. --- .../vtk-field-wrappers-higher-order.md | 13 +++ .../src/tack/interop/vtk/_lagrange.py | 14 ++- .../src/tack/interop/vtk/algorithms.py | 60 ++++++++++++- .../tack-vis/tests/test_vtk_algorithms.py | 89 +++++++++++++++++++ 4 files changed, 171 insertions(+), 5 deletions(-) create mode 100644 changelog.d/vtk-field-wrappers-higher-order.md diff --git a/changelog.d/vtk-field-wrappers-higher-order.md b/changelog.d/vtk-field-wrappers-higher-order.md new file mode 100644 index 000000000..031260f4d --- /dev/null +++ b/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. diff --git a/packages/tack-vis/src/tack/interop/vtk/_lagrange.py b/packages/tack-vis/src/tack/interop/vtk/_lagrange.py index 41521732d..14e100c5d 100644 --- a/packages/tack-vis/src/tack/interop/vtk/_lagrange.py +++ b/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: diff --git a/packages/tack-vis/src/tack/interop/vtk/algorithms.py b/packages/tack-vis/src/tack/interop/vtk/algorithms.py index 6abe1e557..e9827058f 100644 --- a/packages/tack-vis/src/tack/interop/vtk/algorithms.py +++ b/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 diff --git a/packages/tack-vis/tests/test_vtk_algorithms.py b/packages/tack-vis/tests/test_vtk_algorithms.py index 576b60187..10d255286 100644 --- a/packages/tack-vis/tests/test_vtk_algorithms.py +++ b/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", [