Visualization Guide¶

tack.dataset (the tack-vis package) holds meshes and point clouds and the fields on them, and the algorithms that work on both: contour, slice, clip, threshold, external faces, gradients, projections between spaces, streamlines, point-cloud cleaning and gridding, SPH interpolation, halo finding and the rest. Like VTK and Viskores, it works on structured grids, unstructured meshes of the linear cells, polyhedra, polygons, polylines and point clouds. Unlike them, the algorithms are ordinary Tack kernels, so they run on every backend, and the kernels you write yourself see a mesh the way the built-in ones do.

uv sync --extra cpu      # in a clone; both packages are installed

This overview is a tour; the chapters below develop the data model and complete workflows. The Dataset API reference has the model and full support tables, Algorithms every algorithm's parameters, and Data Model every class's.

Quick start¶

For this first example, use uniform_grid, a shortcut factory that assembles a structured topology and uniform geometry for you. Add a scalar field, then extract an isosurface. Building a dataset below shows how to construct these pieces explicitly for several mesh types.

import tack
import tack.dataset as td

tack.set_default_device("cpu")

data = td.uniform_grid((17, 17, 17), origin=(-1, -1, -1), spacing=(0.125, 0.125, 0.125))
positions = data.positions()                   # corner positions as a NumPy array
values = tack.array((positions * positions).sum(axis=1), dtype=tack.f32)
data.fields["radius_squared"] = td.Field(td.H1(data), values)

surface = td.contour(data, "radius_squared", 0.5)
print(surface.num_cells, "triangles,", surface.num_points, "points")
print(td.statistics(data, "radius_squared").mean)

Save this as a .py file and run it after installing the packages above, or use uv run python from a checkout installed as in Getting Started. The scalar coefficients are sampled from x*x + y*y + z*z; the contour is an approximation to a sphere. This example needs no VTK installation.

Learn data and visualization¶

Chapter	Start here when you want to…
Function Space Fundamentals	Understand H1, H2, L2, polynomial order and degrees of freedom
Understand a Dataset	Connect raw storage to topology, geometry and field meaning
Build Meshes	Construct grids, mixed cells, polyhedra, polygons or paths
Point Clouds and Particles	Clean and grid LiDAR, resample SPH particles, find halos
AMR Hierarchies	Hold levels of refined boxes, images or not; reduce and filter them
Octrees	Hold tree-based AMR (RAMSES, HyperTreeGrids); weighted statistics
Fields and Spaces	Preserve a discontinuity, project a field or use implicit values
Physical Units	Give fields and positions units, keep them through filters, convert
Filter and Analyze	Choose between threshold, clip, contour and slice; carry fields and compute derivatives
Write Mesh Algorithms	Add cell and face kernels with correct ids, subsets and orientation
Write a Filter	Write a filter of your own that keeps origins, units, ghost marks and hidden cells as the built-ins do
Sample and Connect Simulations	Probe a solution, integrate paths or exchange datasets
Streaming Steps	Analyze a simulation step by step, in situ: running statistics, time derivatives, pathlines
Work with Data in Pieces	Filter and reduce a dataset cut into pieces, with ghost cells
Algorithms Across Pieces	Pieces from VTK, on several devices or MPI ranks; probes, particles, regions and point clouds across them
Communicate Between Blocks	Write an algorithm over pieces: exchanges, ghost plans, routing and sorting with `tack.communication`
Worked Examples	Run complete programs with analytic checks

The Dataset API reference describes every class and algorithm; Algorithms and Data Model give each one's signature, parameters and methods, from its docstring. The design records, Dataset API and Polyhedra, say how the model came to be. Kernels themselves -- the language, arrays, backends -- are in the core User's Guide.

For a first reading, follow the fundamentals, dataset and mesh chapters, then the field and filter chapters. The field-space comparison and filter comparison show why two plausible images can represent different numerical operations.

Common conventions¶

Every algorithm takes the dataset first and returns something new: a dataset, a field, or Python values. Nothing changes its input.

Fields passed to an algorithm belong to that dataset's topology. dot, cross and composite require matching value types; an explicit tack.arrays.Cast(field.values, tack.f32) converts values when needed. Integer inputs to magnitude, logarithms and statistics compute in floating point.

The model¶

A DataSet is four things:

Part	What it is	Example
`topology`	the cells and how they connect	`UnstructuredTopology`, `StructuredTopology`, `PolyhedralTopology`
`geometry`	where the points are: a field of 3-vectors	positions per point, or the implicit coordinates of a grid
`fields`	named fields, each a space and its values	`Field(H1(data), values)`
`sets`	named selections of points, cells or faces	the boundary faces

A space says which degrees of freedom define a field and how to evaluate it. H1 order 1 shares corner values; H1 order 2 adds shared higher-order nodes. Constant is a separate constant per cell; L2 can keep each cell's own coefficients. Values(data, on) describes attributes without a basis. Every field chooses its own space on the common topology, including geometry. The dataset chapter demonstrates 12 corner positions, 45 quadratic temperature values and 16 discontinuous pressure values on the same two cells.

A field's values are a tack.Array with a scalar or vector element dtype, or an implicit array that computes values on demand. Uniform coordinates need only dimensions, origin, spacing and direction metadata.

Geometry may also live in L2, giving each cell its own positions. Contour and slice preserve the resulting gaps; clipping requires shared H1 geometry. Changing an L2 geometry's values refreshes subsequent face-geometry results.

Building a dataset¶

A dataset can be assembled from simulation arrays, a file reader, generated coordinates, or another source. The construction follows the same pattern:

Create a topology describing cells and their corner incidence.
Create a geometry field: a space on that topology plus coordinate values.
Pass both to td.DataSet, then attach solution fields and named sets.

The following recipes build their topology and geometry explicitly. They assume import numpy as np, import tack, import tack.dataset as td, and an initialized backend. The complete mesh-construction program and unstructured tour include that setup.

Uniform imageRectilinearCurvilinearUnstructured

Construct structured incidence and compute corner positions from origin, spacing and logical indices. The geometry is an H1 field whose values are implicit UniformCoordinates.

    point_dims = (3, 2, 2)
    uniform_topology = td.StructuredTopology(point_dims)
    uniform_coordinates = tack.arrays.UniformCoordinates(point_dims, spacing=(0.5, 1.0, 1.0))
    uniform_geometry = td.Field(td.H1(uniform_topology), uniform_coordinates)
    uniform = td.DataSet(uniform_topology, uniform_geometry)

Construct structured incidence and compute corner positions from three coordinate axes. The geometry uses CartesianProduct values in H1.

    x_axis, y_axis, z_axis = [0.0, 0.25, 1.0], [0.0, 1.0], [0.0, 1.0]
    rectilinear_topology = td.StructuredTopology((len(x_axis), len(y_axis), len(z_axis)))
    rectilinear_coordinates = tack.arrays.CartesianProduct(x_axis, y_axis, z_axis)
    rectilinear_geometry = td.Field(td.H1(rectilinear_topology), rectilinear_coordinates)
    rectilinear = td.DataSet(rectilinear_topology, rectilinear_geometry)

Construct structured incidence and supply explicit positions at its corners. These positions bend the grid while retaining logical connectivity.

    curvilinear_topology = td.StructuredTopology((3, 2, 2))
    # Explicit corner positions in x-fast order, with a bend in z.
    positions = np.array([[0.5 * i, j, k + 0.1 * i * j]
                          for k in range(2) for j in range(2) for i in range(3)], np.float32)
    curvilinear_coordinates = tack.array(positions, dtype=tack.vec3())
    curvilinear_geometry = td.Field(td.H1(curvilinear_topology), curvilinear_coordinates)
    curvilinear = td.DataSet(curvilinear_topology, curvilinear_geometry)

Supply the cell types, offsets and corner ids explicitly, along with an H1 coordinate field. This example has a hexahedron and a pyramid sharing a face. The cell type numbers follow VTK's convention; VTK is not required.

    # A hexahedron with a pyramid on top: VTK cell types, offsets, connectivity.
    positions = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
                          [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1],
                          [0.5, 0.5, 1.6]], np.float32)
    topology = td.UnstructuredTopology(types=[12, 14], offsets=[0, 8, 13],
                                       connectivity=[0, 1, 2, 3, 4, 5, 6, 7, 4, 5, 6, 7, 8])
    coordinate_values = tack.array(positions, dtype=tack.vec3())
    geometry = td.Field(td.H1(topology), coordinate_values)
    mesh = td.DataSet(topology, geometry)

All four recipes use the same DataSet(topology, geometry) constructor. The mesh chapter explains their coordinate representations and also builds polyhedral, polygonal and polyline datasets from their own connectivity and geometry.

The rest of this tour uses mesh, the hex-and-pyramid dataset from the Unstructured tab. Attach height in H1 and a material constant per cell:

    height = tack.empty((mesh.num_points,))
    tack.copyto(height, positions[:, 2])
    mesh.fields["height"] = td.Field(td.H1(mesh), height)            # point data
    material = tack.empty((mesh.num_cells,))
    tack.copyto(material, np.array([1.0, 2.0], np.float32))
    mesh.fields["material"] = td.Field(td.Constant(mesh), material)  # cell data

Factory methods are shortcuts¶

Once the pieces are familiar, uniform_grid and rectilinear_grid save the repeated topology and geometry setup. These calls construct the same uniform and rectilinear datasets as the explicit recipes above:

    # Shortcuts for the same topology-and-geometry combinations built above.
    uniform_shortcut = td.uniform_grid((3, 2, 2), spacing=(0.5, 1.0, 1.0))
    rectilinear_shortcut = td.rectilinear_grid([0.0, 0.25, 1.0], [0.0, 1.0], [0.0, 1.0])

DataSet(topology, positions) is another convenience: when given coordinate values rather than a Field, it wraps them in H1 order-1 geometry. wavelet and tangle are sample-data factories that also initialize known scalar fields, useful for trying filters. They are examples of data sources; a dataset's topology and geometry do not depend on where its values originated.

vtk_to_dataset and tack.interop.mfem import datasets from VTK and MFEM. Neither library is required for the constructions above. In a checkout, uv sync --extra cpu --extra interop installs the optional readers where supported. See dataset interoperability.

Algorithms¶

Every algorithm is in tack.dataset.algorithms, re-exported from tack.dataset, and they share one set of conventions:

the dataset comes first;
a field argument is a Field or the name of one of the dataset's fields;
an algorithm that makes a dataset takes fields= (all by default, a name, a list, or []) and carries those fields onto its output;
one that computes a field returns a Field, or a named tuple of them;
summaries (statistics, entropy) return Python values.

    data = td.wavelet(extent=10)                     # 21^3 points, field "RTData"
    surface = td.contour(data, "RTData", 157.0)      # triangles, RTData carried
    inside = td.threshold(data, "RTData", 200.0, 300.0)
    skin = td.extract_surface(inside)                 # the kept region's boundary
    cut = td.slice_plane(data, origin=(0, 0, 0), normal=(0, 0, 1))
    stats = td.statistics(data, "RTData")            # a Statistics of Python values
    print(f"contour: {surface.num_cells} triangles; threshold: {inside.num_cells} "
          f"cells, {skin.num_cells} boundary faces; mean RTData {stats.mean:.1f}")

Fields come along by how the output's entities come from the input's. Point data is interpolated onto a contour's points, and cell data goes to the pieces each cell was cut into. A kept cell keeps its values.

Moving fields between spaces¶

project(data, field, space) gives a field in another space of the same mesh: td.H1(data) for point values, td.Constant(data) for one value per cell, td.L2(data) for each cell's own copies.

    on_points = td.project(mesh, "material", td.H1(mesh))       # Constant -> H1: averaged
    on_cells = td.project(mesh, "height", td.Constant(mesh))           # H1 -> Constant: corner means
    dg = td.project(mesh, "height", td.L2(mesh))     # each cell's own corner values
    slope = td.gradients(mesh, "height")             # a 3-vector per cell

Projecting into a space that can hold the field exactly gives it back exactly: a linear point field into L2, or into order-2 H1. Otherwise project averages the way VTK does: a point takes the mean of the cells around it, and a cell the mean of its corners.

Sets¶

A Set selects points, cells or faces. Builders such as boundary_faces return one; name it in data.sets to keep it. Filters carry sets the way they carry fields:

    mesh.sets["boundary"] = td.boundary_faces(mesh)  # a face Set
    mesh.sets["roof"] = td.Set(mesh, "cells", [1])   # a cell Set
    kept = td.threshold(mesh, "material", 1.5, 3.0)  # sets come along
    print(f"boundary: {len(mesh.sets['boundary'])} faces; the roof is cell "
          f"{kept.sets['roof'].numpy()} of the threshold")

The name of a cell or face set is also a place to run a kernel, as the next section shows.

Writing your own algorithm¶

for_each(kernel, data, domain, *args) runs a kernel over "cells", "faces", "edges" or a named set. A mesh may mix shapes, so it runs the kernel once per shape. Each launch gets that shape's view of the cells, and every Field argument becomes its view for the same cells. The kernel is compiled for each shape, so a hexahedron's loop knows it has eight corners.

for_each is for kernels that need the mesh: a cell's corners, a face's two sides, positions. A kernel that needs only a field's values, each on its own, is an ordinary kernel over field.values, looping over range(len(values)). That is also the loop over points, which are not a for_each domain: point data is one value per point. Implicit values, such as a uniform grid's coordinates, are indexed the same way.

    @tack.kernel
    def centers(cells, u, out, where):
        # Once per shape: `cells` is that shape's view, `u` the field's.
        for c in cells:
            pc = cells.parametric_center()
            out[cells.entity_id(c)] = u.value(c, pc)
            where[cells.entity_id(c)] = cells.position(c, pc)

    values = tack.empty((mesh.num_cells,))
    where = tack.empty((mesh.num_cells,), tack.vec3(tack.f32))
    td.for_each(centers, mesh, "cells", mesh.fields["height"], values, where)

A view's methods follow one naming rule per kind of answer, so a name says what comes back:

Ends in	Gives	Examples
`_id`	a global id	`entity_id(c)`, `point_id(c, j)`, `face_id(c, k)`, `side_point_id(c, k, j)`
`_position`	a position	`position(c, pc)`, `corner_position(c, j)`, `side_position(c, k, j)`
`_corner`	a shape's local corner number	`edge_corner(e, k)`, `face_corner(f, k)`
`num_` (prefix)	a count	`num_points`, `num_faces(c)`, `side_num_points(c, k)`

A field view gives value(c, pc) at parametric coordinates, parametric_gradient(c, pc) and corner_value(c, j).

Face-based kernels walk each cell's faces outward, and the same code serves every kind of cell. They need the faces, which are derived on first use, so ask for them before the launch:

    @tack.kernel
    def surface_areas(cells, out):
        # Every face of every cell, walked outward from the cell: the same code
        # for a hexahedron, a pyramid or a polyhedron.
        for c in cells:
            total = 0.0
            for k in range(cells.num_faces(c)):
                a = cells.side_position(c, k, 0)
                n = a - a
                for j in range(1, cells.side_num_points(c, k) - 1):
                    n += (cells.side_position(c, k, j) - a).cross(
                        cells.side_position(c, k, j + 1) - a)
                total += 0.5 * n.norm()
            out[cells.entity_id(c)] = total

    mesh.topology.faces()                            # the faces kernels walk
    areas = tack.empty((mesh.num_cells,))
    td.for_each(surface_areas, mesh, "cells", areas)

    polyhedra = td.as_polyhedra(mesh)                # every cell a polyhedron
    same = tack.empty((polyhedra.num_cells,))
    td.for_each(surface_areas, polyhedra, "cells", same)

Polyhedra, polygons and lines¶

A polyhedral cell has no reference element: no shape functions and no parametric coordinates. Instead, each face is fanned into triangles from the mean of its points, and the triangles into tetrahedra from the mean of the cell's points; point data is linear on each tetrahedron. Gradients, cell centers, probing, advection and streamlines run on polyhedra that way (gradients by Green–Gauss over the fanned faces), exactly for linear fields. Polygons are fanned into triangles from their points' mean the same way, so probing and advection run on them, and probing locates on a polyline's segments. The algorithms that need a reference element (L2 and order-2 spaces, side traces) refuse a polyhedral mesh and say so. The ones that work from points and faces run on both kinds: contour (by López's face walk), slice, threshold, external faces, cell and face geometry, divergence, and moving data between points and cells. A 2D mesh contours into lines, and streamlines are polylines, one per seed:

    flow = td.uniform_grid((9, 9, 3), origin=(-1, -1, 0), spacing=(0.25, 0.25, 0.5))
    x = flow.positions()
    swirl = tack.empty((flow.num_points,), tack.vec3(tack.f32))
    tack.copyto(swirl, np.c_[-x[:, 1], x[:, 0], np.zeros(len(x))])
    flow.fields["velocity"] = td.Field(td.H1(flow), swirl)
    paths = td.streamlines(flow, "velocity", seeds=[[0.5, 0, 0.5], [0, 0.75, 0.5]],
                           step_size=0.05, steps=40)
    lengths = td.cell_geometry(paths).measure        # a polyline's length

The reference has a table of which algorithm runs on which topology.

64-bit ids¶

A topology keeps its ids in i32 unless it is too large or asked for i64 (id_dtype=tack.i64). Kernels specialize on the id type, so a small mesh pays nothing for large ones. VTK's 64-bit ids come in as the narrowest type that holds them.

Running the examples¶

The whole example is docs/examples/datasets.py; --check runs its assertions.

The new chapter examples can be checked together with the rest of the guide:

uv run python docs/examples/validate.py --arch cpu

They cover mesh counts and measures, continuous/discontinuous fields, analytic gradients, flux balance, exact clipping volumes and flow integration. Their optional comparison figures use Matplotlib, independently of VTK.

VTK interop¶

tack.interop.vtk converts whole datasets: - vtk_to_dataset reads a vtkUnstructuredGrid, vtkRectilinearGrid or vtkImageData, including polyhedra and polylines; - dataset_to_vtk writes the matching VTK type back.

from tack.interop.vtk import dataset_to_vtk, vtk_to_dataset

data = vtk_to_dataset(grid)               # point data -> H1, cell data -> Constant
grid = dataset_to_vtk(td.contour(data, "pressure", 0.5))

Arrays exchange through DLPack, sharing storage where the device allows:

from tack.interop.vtk import array_to_vtk, vtk_to_array

values = vtk_to_array(vtk_array)            # a tack.Array sharing VTK's storage
vtk_array = array_to_vtk(values)

VTK's tuples and components become an array shape: a (1000, 3) scalar array represents 1000 tuples of three components, and an array of vec3 elements exports as 1000 tuples of three as well. vtk_to_array(array, flatten=True) and array_to_vtk(values, n_components=3) use the flat interleaved layout instead. This requires a VTK build exposing vtkmodules.util.dlpack_support. Not every VTK installation provides device arrays. See Sharing arrays with other libraries for lifetime and read-only rules.

Level Zero¶

Level Zero device pointers need a shared allocation context. With a compatible VTK/Viskores Kokkos-SYCL build, initialize before allocating arrays:

from tack.interop.vtk import init_level_zero

init_level_zero()     # instead of tack.set_default_device("level_zero")

The VTK DLPack support must provide level_zero_handles(). Fields created in a separate Tack context cannot be exchanged through this path. For another library, tack.device("level_zero", external_context={"driver": ..., "device": ..., "context": ...}) makes a device in its context; Tack never destroys an adopted context.