Data Model

The classes and constructors of tack.dataset, and the interop with VTK, MFEM and Warp, with their parameters, methods and properties. The Dataset API reference explains how they fit together; the algorithms that work on them are in Algorithms.

Contents

• Datasets and fields: DataSet, Field, the constructors and sources
• Data in pieces and hierarchies: PartitionedDataSet, AMRDataSet, OctreeDataSet
• Topologies: structured, unstructured, polyhedral, polygonal, polyline; cell shapes
• Spaces: where a field's values live
• Sets, ghost marks and global ids
• Implicit functions
• Locating points: CellLocator, PointLocator
• Steps in time: temporal filters and particle paths
• Kernels over meshes: for_each
• Filters of your own: carry and its maps, ghost marks, hidden cells, owned values
• Units: parameters as quantities, units of what a filter computes
• Interoperability: VTK, MFEM, Warp

Datasets and fields

tack.dataset.DataSet

DataSet(
    topology,
    geometry,
    fields=None,
    sets=None,
    dtype=tack.f32,
    ghosts=None,
    global_ids=None,
    origin=None,
    unit=None,
    time=None,
    time_unit=None,
)

A topology, a geometry, named fields and named sets. The geometry is a Field in H1 or L2 of 3-vectors -- or positions per point (a host array, or any array of 3-vectors, such as a device array or a CartesianProduct), which become an H1 field -- kept apart from fields as data.geometry: algorithms move it with the cells, never as one of the fields they carry or select. ghosts, a tack.dataset.Ghosts or None, marks the cells and points another piece of a partitioned dataset owns; global_ids, a tack.dataset.GlobalIds or None, numbers them across the pieces.

origin, three floats kept in double precision on the host, is where the geometry's coordinates start: a position is origin + value, the value in the geometry's type (tack.f32, or tack.f64 on a device that has it). Kernels see only the values, so coordinates far from zero -- a LiDAR survey's eastings and northings -- keep their precision in single precision. Filters keep it; what takes or gives absolute positions (implicit functions, probe points, seeds, world_positions, the VTK interop) converts on the host in double precision.

unit is the spatial unit, the geometry field's (DataSet.unit), in which origin is too; a geometry Field brings its own.

time is when the data are, a double on the host or None, in time_unit (a unit of time, or None for unknown): one step of a simulation as it streams through filters, which keep both, as they keep the origin. Tack holds no series of steps; filters that need earlier steps keep what they need (tack.dataset.temporal).

Source code in packages/tack-vis/src/tack/dataset/dataset.py

@on_inputs_device
def __init__(self, topology, geometry, fields=None, sets=None, dtype=tack.f32,
             ghosts=None, global_ids=None, origin=None, unit=None, time=None,
             time_unit=None):
    self.topology = topology
    self._geometry = _as_geometry(geometry, topology, dtype, unit)
    self.origin = origin
    self.time = time
    self.time_unit = time_unit
    self.fields = dict(fields or {})
    for name, field in self.fields.items():
        if field.space.topology is not topology:
            raise ValueError(f"field {name!r} is on another topology")
    self.sets = sets
    # Set by a filter that dropped its input's ghost marks (ghosts.drop_marks).
    self._ghosts_lost = False
    self.ghosts = ghosts
    self.global_ids = global_ids
    _ = self._device                    # refuses parts on two devices

origin property writable

origin

Where the geometry's coordinates start, (x, y, z) in double precision: a position is origin plus the geometry's value. (0.0, 0.0, 0.0) unless set.

time property writable

time

When the data are, a float, or None: one step's time, in time_unit.

time_unit property writable

time_unit

The unit of time, a unit of time, or None if unknown.

sets property writable

sets

Named Sets of this dataset's points, cells, faces or edges, checked as they are put in: a set of another topology, or named for a domain, is refused.

ghosts property writable

ghosts

The ghost marks of its cells and points (tack.dataset.Ghosts), or None: which cells and points another piece owns, in VTK's vtkGhostType bits.

Filters that keep the cells and points (warp, transform, shrink, as_polyhedra, probe's target, ...) keep the marks; others do not carry them yet, and their output remembers that it lost them, so statistics and entropy refuse it rather than count shared cells and points once per piece. Setting marks again clears that.

global_ids property writable

global_ids

The global ids of its points and cells (tack.dataset.GlobalIds), or None: each one's number across every piece of a partitioned dataset. Filters carry them where an output entity is an input entity (kept points, kept or extracted cells); what a filter makes has none.

device property

device

The tack.Device it is on: where its algorithms run and allocate.

geometry property writable

geometry

The geometry, a field of positions in H1 or L2, relative to origin.

unit property

unit

The spatial unit, the geometry's (origin is in it too), a string, or None if unknown.

dtype property

dtype

The geometry's floating-point type.

id_dtype property

id_dtype

The topology's integer type for ids: tack.i32 or tack.i64.

to_device

to_device(device)

This dataset on device: itself if it is already there, else a copy -- the topology made from its arrays copied there, each field in the same space on it, the sets by the same ids. What the topology derives (faces, groups, layouts) is derived there afresh, not copied.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def to_device(self, device):
    """This dataset on ``device``: itself if it is already there, else a copy --
    the topology made from its arrays copied there, each field in the same
    space on it, the sets by the same ids. What the topology derives (faces,
    groups, layouts) is derived there afresh, not copied."""
    from tack.dataset.sets import Set

    target = as_device(backend_of(device))
    here = self._device
    if here is not None and target._backend.memory_domain == here.memory_domain:
        return self
    topology = self.topology.to_device(target)

    def moved(field):
        return Field(field.space._rebuilt_on(topology), field.values.to_device(target),
                     field.unit)

    sets = {name: Set._from_sorted_ids(topology, s.on, s.ids.to_device(target))
            for name, s in self.sets.items()}
    ghosts = None if self.ghosts is None else self.ghosts.to_device(target)
    global_ids = None if self.global_ids is None else self.global_ids.to_device(target)
    out = DataSet(topology, moved(self.geometry),
                  {name: moved(field) for name, field in self.fields.items()},
                  sets, dtype=self.dtype, ghosts=ghosts, global_ids=global_ids,
                  **self._frame())
    out._ghosts_lost = self._ghosts_lost
    return out

positions

positions()

Every point's position as a host array, (num_points, 3): an H1 geometry's values, relative to origin (world_positions adds it). An L2 geometry has none per point.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def positions(self):
    """Every point's position as a host array, ``(num_points, 3)``: an ``H1``
    geometry's values, relative to ``origin`` (``world_positions`` adds it).
    An ``L2`` geometry has none per point."""
    if isinstance(self.geometry.space, L2):
        raise ValueError("an L2 geometry has positions per cell corner, not per point")
    # An order-2 geometry's first values are the points'.
    return arrays.to_host(self.geometry.values)[:self.num_points]

rebased

rebased(origin)

This dataset with its geometry relative to origin instead: the same world positions, each value moved by the difference of the two origins (taken in double precision, added in the geometry's type). An image stays an image, its corner moved in double precision, so an image at survey coordinates rebased to a nearby origin holds its points exactly; a rectilinear grid moves its axes. Other implicit geometry becomes positions per point.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def rebased(self, origin):
    """This dataset with its geometry relative to ``origin`` instead: the same
    world positions, each value moved by the difference of the two origins
    (taken in double precision, added in the geometry's type). An image
    stays an image, its corner moved in double precision, so an image at
    survey coordinates rebased to a nearby origin holds its points exactly;
    a rectilinear grid moves its axes. Other implicit geometry becomes
    positions per point."""
    from tack.arrays import CartesianProduct, UniformCoordinates
    from tack.dataset.algorithms._moving import _moved, _shifted

    origin = as_origin(origin)
    if origin == self.origin:
        return self
    offset = np.asarray(self.origin) - np.asarray(origin)
    values = self.geometry.values
    if isinstance(values, UniformCoordinates):
        moved = UniformCoordinates(values.dims, np.asarray(values._origin) + offset,
                                   values._spacing,
                                   values._matrix if values.ORIENTED else None,
                                   dtype=values.dtype, device=self.device)
    elif isinstance(values, CartesianProduct):
        axes = [arrays.to_host(axis).astype(np.float64).reshape(-1) + offset[a]
                for a, axis in enumerate((values.xs, values.ys, values.zs))]
        moved = CartesianProduct(*axes, dtype=values.dtype, device=self.device)
    else:
        moved = _shifted(self, offset)
    out = _moved(self, moved)
    out.origin = origin
    return out

world_positions

world_positions()

Every point's position as a float64 host array, (num_points, 3): origin plus positions(), added in double precision.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def world_positions(self):
    """Every point's position as a float64 host array, ``(num_points, 3)``:
    ``origin`` plus ``positions()``, added in double precision."""
    return self.positions().astype(np.float64) + np.asarray(self.origin)

domain_groups

domain_groups(domain)

The DomainGroups of domain: "cells", "faces", "edges", or the name of a cell, face or edge set in sets -- a cell set's as the cell groups' subgroups of its members.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

@on_inputs_device
def domain_groups(self, domain):
    """The ``DomainGroup``s of ``domain``: ``"cells"``, ``"faces"``, ``"edges"``, or
    the name of a cell, face or edge set in ``sets`` -- a cell set's as the
    cell groups' subgroups of its members."""
    kind = self._kind(domain)
    if domain == "cells":
        return self.topology.groups()
    if domain == "faces":
        return self.topology.faces().groups()
    if domain == "edges":
        return self.topology.edges().groups()
    if kind == "faces":
        return self.topology.faces().groups(self.sets[domain].ids)
    if kind == "edges":
        return self.topology.edges().groups(self.sets[domain].ids)
    return self.launch_groups(domain)

domain_view

domain_view(domain, group, extra=(), **extra_attributes)

group's view: kind, shape, geometry, and for cells their incidence when faces and edges have been derived; extra mixins and attributes last (a space's node mixins, node_views).

Source code in packages/tack-vis/src/tack/dataset/dataset.py

@on_inputs_device
def domain_view(self, domain, group, extra=(), **extra_attributes):
    """``group``'s view: kind, shape, geometry, and for cells their incidence when
    faces and edges have been derived; ``extra`` mixins and attributes last
    (a space's node mixins, ``node_views``)."""
    mixins = []
    kind = self._kind(domain)
    geometry, attributes = _geometry_parts(self, group, kind)
    mixins.extend(geometry)
    if kind == "cells" and issubclass(group.kind, views._CellsById):
        mixins.append(views._PolyhedralWalk)
        if group.shape is shapes.Polyhedron:
            mixins.append(views._PolyhedralLocation)
        elif group.shape is shapes.Polygon:
            mixins.append(views._PolygonLocation)
        elif group.shape is shapes.PolyLine:
            mixins.append(views._PolylineLocation)
    if (kind == "cells" and issubclass(group.kind, views._PolylineCells)
            and self.topology._faces is not None):
        mixins.append(views._PolylineIncidence)
        attributes.update(cell_faces=self.topology.cell_faces,
                          cell_face_sides=self.topology.cell_face_sides)
    if kind == "cells":
        # A subgroup's index(c) is its position in the topology's group, so
        # it shares that group's incidence.
        faces = self.topology._faces
        # The cells the faces were derived from: those of the top dimension.
        if faces is not None and group.root in getattr(faces, "group_starts", {}):
            mixins.extend((views._FaceIncidence, views._FaceWalk))
            attributes.update(side_face=faces.side_face, side_slot=faces.side_slot,
                              side_orientation=faces.side_orientation,
                              face_start=faces.group_starts[group.root])
        elif group.shape.NUM_FACES and not issubclass(group.kind, (views._CellsById,
                                                                   views._PolylineCells)):
            # A shape's faces walk from its tables alone, with no derived
            # faces: a fresh mesh's volumes need not sort them.
            mixins.append(views._FaceWalk)
        edges = self.topology._edges
        if edges is not None and group.shape.NUM_EDGES:
            mixins.append(views._EdgeIncidence)
            attributes.update(side_edge=edges.side_edge, side_sign=edges.side_sign,
                              edge_start=edges.group_starts[group.root])
    # A key that shapes the kernel -- a size bucket's cap -- brings its mixin.
    for key, value in group.keys.items():
        if hasattr(key, "domain_mixin"):
            mixins.append(key.domain_mixin(value))
    mixins.extend(extra)
    attributes.update(extra_attributes)
    return group.view(*mixins, **attributes)

node_views

node_views(space)

Views over space's nodes, each visited once: entity views -- points, edges, quad faces, cells -- whose node_count(i) nodes have node_dof(i, j), where the node's value is in space's values, and node_position(i, j), where the geometry puts it. for_each(kernel, data, space, ...) runs a kernel over them:

for i in nodes:
    for j in range(nodes.node_count(i)):
        values[nodes.node_dof(i, j)] = f(nodes.node_position(i, j))

H1's nodes are its points and, at order 2, its edges' middles, its quad faces' centers and its hexahedra's and quads' interiors; a shared node's position is the mean of its entity's points on a straight geometry, the geometry's own value there on an order-2 one. L2's are its cells' corners (order 1) or parametric centers (order 0), Constant's and values on cells' the parametric centers, values on points' the points. Shared nodes need positions per point: an L2 geometry has none, and is refused for them. On polyhedra, polygons and polylines only point nodes exist.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

@on_inputs_device
def node_views(self, space):
    """Views over ``space``'s nodes, each visited once: entity views -- points,
    edges, quad faces, cells -- whose ``node_count(i)`` nodes have
    ``node_dof(i, j)``, where the node's value is in ``space``'s values, and
    ``node_position(i, j)``, where the geometry puts it. ``for_each(kernel,
    data, space, ...)`` runs a kernel over them:

        for i in nodes:
            for j in range(nodes.node_count(i)):
                values[nodes.node_dof(i, j)] = f(nodes.node_position(i, j))

    ``H1``'s nodes are its points and, at order 2, its edges' middles, its
    quad faces' centers and its hexahedra's and quads' interiors; a shared
    node's position is the mean of its entity's points on a straight
    geometry, the geometry's own value there on an order-2 one. ``L2``'s
    are its cells' corners (order 1) or parametric centers (order 0),
    ``Constant``'s and values on cells' the parametric centers, values on
    points' the points. Shared nodes need positions per point: an ``L2``
    geometry has none, and is refused for them. On polyhedra, polygons and
    polylines only point nodes exist."""
    from tack.dataset.spaces import Constant, Values, reference_cells

    if space.topology is not self.topology:
        raise ValueError("a space of another topology")
    geometry = self.geometry
    point_nodes = ((isinstance(space, H1) and space.order == 1)
                   or (isinstance(space, Values) and space.on == "points"))
    shared = point_nodes or isinstance(space, H1)
    if shared and not isinstance(geometry.space, H1):
        raise NotImplementedError(
            f"the nodes of {space!r} are shared by cells, and an L2 geometry gives a "
            "shared node no single position")
    out = []
    if shared and self.num_points:
        group = views.DomainGroup(views._Points, shapes.Vertex, (self.num_points,),
                            self.num_points, 0)
        mixins, attributes = _geometry_parts(self, group, "points")
        out.append(group.view(*mixins, views._PointDof, views._MeanPosition, **attributes))
    if point_nodes:
        return out
    if isinstance(space, Values) and space.on not in ("points", "cells"):
        raise TypeError(f"values on {space.on} have no nodes to visit as such; "
                        f"for_each over {space.on!r} visits them")
    reference_cells(self.topology, f"the nodes of {space!r}")
    if isinstance(space, H1):
        face_dofs, cell_dofs, _ = space._layout()
        curved = geometry.space.order == 2
        where = views._GeometryAtDof if curved else views._MeanPosition
        for group in self.topology.edges().groups():
            if group.count:
                out.append(self.domain_view("edges", group, (views._EdgeDof, where),
                                            node_edge_base=self.num_points))
        if any(g.count and g.shape.QUADRATIC_FACES for g in self.topology.groups()):
            for group in self.topology.faces().groups():
                if group.count and group.shape is shapes.Quad:
                    out.append(self.domain_view("faces", group, (views._FaceDof, where),
                                                node_face_dofs=face_dofs))
        for group in self.topology.groups():
            if group.count and group.shape.QUADRATIC_INTERIOR:
                out.append(self.domain_view("cells", group, (views._InteriorDof, where),
                                            node_cell_dofs=cell_dofs))
        return out
    groups = (self.launch_groups("cells", keys=(space,)) if space.varies
              else self.topology.groups())
    for group in groups:
        if not group.count:
            continue
        if isinstance(space, L2):
            order = group.keys[space] if space.varies else space.order
            dof = views._CornerDof if order == 1 else views._OffsetDof
            place = views._CornerPosition if order == 1 else views._CenterPosition
            out.append(self.domain_view("cells", group, (dof, place),
                                        node_offsets=space.offsets))
        elif isinstance(space, (Constant, Values)):
            out.append(self.domain_view("cells", group,
                                        (views._CellDof, views._CenterPosition)))
        else:
            raise TypeError(f"no nodes for {space!r}")
    return out

launch_groups

launch_groups(domain, fields=(), keys=())

The groups a kernel over domain launches once each, given the Fields it reads: the domain's groups, each split by the keys of the spaces among fields that vary from cell to cell, and by any other keys -- objects that give each cell a key, as polyhedra.SizeBuckets does.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

@on_inputs_device
def launch_groups(self, domain, fields=(), keys=()):
    """The groups a kernel over ``domain`` launches once each, given the ``Field``s
    it reads: the domain's groups, each split by the keys of the spaces among
    ``fields`` that vary from cell to cell, and by any other ``keys`` -- objects
    that give each cell a key, as ``polyhedra.SizeBuckets`` does."""
    kind = self._kind(domain)
    members = self.sets[domain] if kind == "cells" and domain != "cells" else None
    groups = self.topology.groups() if members is not None else self.domain_groups(domain)
    spaces = [members] if members is not None else []
    for field in fields:
        if field.space.varies and field.space not in spaces:
            spaces.append(field.space)
    spaces.extend(k for k in keys if k not in spaces)
    if kind != "cells" or not spaces:
        return groups
    # A cell set splits the groups by membership, and keeps its members.
    return [sub for group in groups if group.count
            for sub in _subgroups(group, tuple(spaces))
            if members is None or sub.keys[members] == 1]

tack.dataset.Field

Field(space, values, unit=None)

A space and its values: any array, explicit or implicit, of scalars or vectors, with one value per degree of freedom of the space. Everything else a kernel needs to find them belongs to the space.

unit is the values' physical unit, a string such as "m/s" kept in unyt's canonical spelling (tack.dataset.units), or None for unknown: meaning, which kernels never see. Filters that keep a quantity keep its unit; anything else gives None rather than a guess.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def __init__(self, space, values, unit=None):
    if not isinstance(space, Space):
        raise TypeError(f"a field's space is a Space on a topology, not {space!r}")
    views.storage(values)                       # refuses what is not an array
    if arrays.size_of(values) != space.size:
        raise ValueError(f"{space!r} holds {space.size} values, not "
                         f"{arrays.size_of(values)}")
    self.space = space
    self.values = values
    self._unit = units.canonical(unit)
    _ = self._device                    # refuses values on another device

unit property

unit

The values' unit, a string, or None if unknown.

device property

device

The tack.Device its values are on.

view

view(group)

This field's view for group: the group's kind and shape, the space, and the storage its values are read through. A space that varies reads its key for the group from group.keys: launch it through for_each.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def view(self, group):
    """This field's view for ``group``: the group's kind and shape, the space, and
    the storage its values are read through. A space that varies reads its key
    for the group from ``group.keys``: launch it through ``for_each``."""
    mixins, attributes = _space_parts(self, group)
    return group.view(*mixins, **attributes)

Making datasets

tack.dataset.uniform_grid

uniform_grid(
    dims,
    origin=(0.0, 0.0, 0.0),
    spacing=(1.0, 1.0, 1.0),
    direction=None,
    dtype=tack.f32,
    device=None,
    unit=None,
)

A dataset of dims points (one to three sizes) at origin + direction @ (ijk * spacing): VTK's image data, Viskores' uniform data set. Nothing per point is stored (UniformCoordinates). A size of 1 adds no cells, in any position, and a grid with no axis of two points has none (grid_topology). origin and spacing take a value per size given (or three), and direction is 3 x 3 (2 x 2 for a 2D image). It is on device (by default the default device).

The image's coordinates are in the geometry's type, so one far from zero -- a raster at survey coordinates -- loses digits in single precision: rebased(origin) keeps it an image with its corner relative to a nearby dataset origin, moved in double precision, as a point cloud's origin keeps its points. unit is the spatial unit (DataSet.unit), which origin and spacing are in.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def uniform_grid(dims, origin=(0.0, 0.0, 0.0), spacing=(1.0, 1.0, 1.0), direction=None,
                 dtype=tack.f32, device=None, unit=None):
    """A dataset of ``dims`` points (one to three sizes) at ``origin + direction @
    (ijk * spacing)``: VTK's image data, Viskores' uniform data set. Nothing per
    point is stored (``UniformCoordinates``). A size of 1 adds no cells, in any
    position, and a grid with no axis of two points has none (``grid_topology``).
    ``origin`` and ``spacing`` take a value per size given (or three), and
    ``direction`` is 3 x 3 (2 x 2 for a 2D image). It is on ``device`` (by
    default the default device).

    The image's coordinates are in the geometry's type, so one far from zero
    -- a raster at survey coordinates -- loses digits in single precision:
    ``rebased(origin)`` keeps it an image with its corner relative to a nearby
    dataset origin, moved in double precision, as a point cloud's origin
    keeps its points. ``unit`` is the spatial unit (``DataSet.unit``), which
    ``origin`` and ``spacing`` are in."""
    topology = grid_topology(dims, device=device)   # the sizes as given, before padding
    return DataSet(topology, UniformCoordinates(dims, origin, spacing, direction, dtype=dtype,
                                                device=device), unit=unit)

tack.dataset.rectilinear_grid

rectilinear_grid(
    x,
    y=(0.0,),
    z=(0.0,),
    dtype=tack.f32,
    device=None,
    unit=None,
)

A dataset of the grid whose lines run through x, y and z; see grid_topology for axes of one point or none. It is on device, by default that of axes given as device arrays, else the default device. unit is the spatial unit (DataSet.unit).

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def rectilinear_grid(x, y=(0.0,), z=(0.0,), dtype=tack.f32, device=None, unit=None):
    """A dataset of the grid whose lines run through ``x``, ``y`` and ``z``; see
    ``grid_topology`` for axes of one point or none. It is on ``device``, by
    default that of axes given as device arrays, else the default device.
    ``unit`` is the spatial unit (``DataSet.unit``)."""
    coordinates = CartesianProduct(x, y, z, dtype=dtype, device=device)
    return DataSet(grid_topology(coordinates.dims, device=coordinates.device), coordinates,
                   unit=unit)

tack.dataset.point_cloud

point_cloud(
    positions,
    fields=None,
    origin=None,
    dtype=tack.f32,
    device=None,
    unit=None,
    units=None,
)

A point cloud: a dataset of positions, one vertex cell per point (types, offsets and connectivity implicit, stored nowhere), with fields (a dict of names to one value per point, scalars or rows of components) as point data.

positions is an (n, 3) host array of world positions, or a Tack array of 3-vectors already relative to origin. Host positions are taken relative to origin in double precision before they become dtype; without one, float64 positions get the lower corner of their box, rounded down to whole units, so a LiDAR survey's eastings and northings keep their precision in single precision, and other positions get zero. dtype is tack.f32, or tack.f64 on a device that has it. Floating-point fields become dtype too; integer fields (a LAS classification, intensities) keep their type. unit is the spatial unit (DataSet.unit), which origin is in, and units maps field names to their units.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

def point_cloud(positions, fields=None, origin=None, dtype=tack.f32, device=None, unit=None,
                units=None):
    """A point cloud: a dataset of ``positions``, one vertex cell per point (types,
    offsets and connectivity implicit, stored nowhere), with ``fields`` (a dict of
    names to one value per point, scalars or rows of components) as point data.

    ``positions`` is an ``(n, 3)`` host array of world positions, or a Tack array
    of 3-vectors already relative to ``origin``. Host positions are taken
    relative to ``origin`` in double precision before they become ``dtype``;
    without one, float64 positions get the lower corner of their box, rounded
    down to whole units, so a LiDAR survey's eastings and northings keep their
    precision in single precision, and other positions get zero. ``dtype`` is
    ``tack.f32``, or ``tack.f64`` on a device that has it. Floating-point fields
    become ``dtype`` too; integer fields (a LAS classification, intensities) keep
    their type. ``unit`` is the spatial unit (``DataSet.unit``), which ``origin``
    is in, and ``units`` maps field names to their units."""
    from tack.dataset import ids
    from tack.dataset.algorithms._extract import _vertex_cells
    from tack.dataset.topology import UnstructuredTopology

    target = as_device(backend_of(device))
    if dtype not in target._backend.supported_dtypes:
        raise TypeError(f"the {target} device has no {dtype} arrays; a point cloud there "
                        f"takes positions of {tack.f32} relative to an origin")
    if isinstance(positions, tack.Array):
        n = int(positions.shape[0])
        local = positions
    else:
        world = np.asarray(positions, np.float64 if np.asarray(positions).dtype == np.float64
                           else None).reshape(-1, 3)
        n = len(world)
        if origin is None and world.dtype == np.float64 and n:
            origin = np.floor(world.min(axis=0))
        local = world.astype(np.float64) - np.asarray(as_origin(origin))
    with target:
        topology = UnstructuredTopology(*_vertex_cells(n, ids.choose(None, n)), num_points=n)
        data = DataSet(topology, local, dtype=dtype, origin=origin, unit=unit)
        units = dict(units or {})
        unknown = set(units) - set(fields or {})
        if unknown:
            raise KeyError(f"units for fields not given: {', '.join(map(repr, sorted(unknown)))}")
        for name, values in (fields or {}).items():
            if not isinstance(values, tack.Array):
                host = np.asarray(values)
                if host.shape[:1] != (n,):
                    raise ValueError(f"field {name!r} has {host.shape[:1]} values for "
                                     f"{n} points")
                if host.dtype.kind == "f":
                    host = host.astype(dtype.numpy_dtype)
                width = host.shape[1] if host.ndim == 2 else None
                if width:
                    from tack.lang.types import from_numpy_dtype

                    values = tack.array(np.ascontiguousarray(host),
                                        tack.vec(width, from_numpy_dtype(host.dtype)))
                else:
                    values = tack.array(np.ascontiguousarray(host))
            data.fields[name] = Field(H1(data), values, units.get(name))
    return data

Sources

tack.dataset.sources

Datasets made from formulas, for tests and examples: Viskores' sources (viskores::source), on rectilinear grids with uniform spacing.

wavelet is VTK's and Viskores' wavelet (vtkRTAnalyticSource): a Gaussian plus three periodic terms, the point field RTData. tangle is the tangle cube, a quartic with a surface of several handles at 0.5, the point field tangle.

tack.dataset.wavelet

wavelet(
    extent=10,
    spacing=1.0,
    center=(0.0, 0.0, 0.0),
    maximum=255.0,
    frequency=(60.0, 30.0, 40.0),
    magnitude=(10.0, 18.0, 5.0),
    standard_deviation=0.5,
    dtype=tack.f32,
    device=None,
)

VTK's and Viskores' wavelet on the grid of point indices extent (n: -n to n on every axis, or x0, x1, y0, y1, z0, z1), points at index times spacing: maximum * exp(-|a|^2 / (2 sd^2)) plus magnitude times sin, sin and cos of frequency * a, where a is center minus the point, each axis divided by its extent's length. It is on device (by default the default device).

Source code in packages/tack-vis/src/tack/dataset/sources.py

def wavelet(extent=10, spacing=1.0, center=(0.0, 0.0, 0.0), maximum=255.0,
            frequency=(60.0, 30.0, 40.0), magnitude=(10.0, 18.0, 5.0),
            standard_deviation=0.5, dtype=tack.f32, device=None):
    """VTK's and Viskores' wavelet on the grid of point indices ``extent`` (``n``:
    -n to n on every axis, or ``x0, x1, y0, y1, z0, z1``), points at index times
    ``spacing``: ``maximum * exp(-|a|^2 / (2 sd^2))`` plus ``magnitude`` times
    sin, sin and cos of ``frequency * a``, where ``a`` is ``center`` minus the
    point, each axis divided by its extent's length. It is on ``device`` (by
    default the default device)."""
    lo, hi = _extent(extent)
    s = np.broadcast_to(np.asarray(spacing, float), (3,))
    data = uniform_grid([int(hi[d] - lo[d] + 1) for d in range(3)],
                        [lo[d] * s[d] for d in range(3)], s, dtype=dtype, device=device)
    n = data.num_points
    out = tack.empty((n,), dtype, device=data.device)
    scale = [1.0 / (hi[d] - lo[d]) if hi[d] > lo[d] else 1.0 for d in range(3)]
    _wavelet(out, int(hi[0] - lo[0] + 1), int(hi[1] - lo[1] + 1), *map(float, lo),
             *map(float, s), *map(float, center), *scale, *map(float, frequency),
             *map(float, magnitude), float(maximum),
             1.0 / (2.0 * standard_deviation * standard_deviation))
    data.fields["RTData"] = Field(H1(data), out)
    return data

tack.dataset.tangle

tangle(dims=(16, 16, 16), dtype=tack.f32, device=None)

Viskores' tangle on dims points per axis over [0, 1]^3: the quartic (x^4 - 5x^2 + y^4 - 5y^2 + z^4 - 5z^2 + 11.8) / 5 + 1/2 of the point mapped to [-3, 3]^3, the point field tangle. It is on device (by default the default device).

Source code in packages/tack-vis/src/tack/dataset/sources.py

def tangle(dims=(16, 16, 16), dtype=tack.f32, device=None):
    """Viskores' tangle on ``dims`` points per axis over [0, 1]^3: the quartic
    ``(x^4 - 5x^2 + y^4 - 5y^2 + z^4 - 5z^2 + 11.8) / 5 + 1/2`` of the point
    mapped to [-3, 3]^3, the point field ``tangle``. It is on ``device`` (by
    default the default device)."""
    dims = np.broadcast_to(np.asarray(dims, int), (3,))
    if (dims < 2).any():
        raise ValueError("tangle needs at least two points per axis")
    cells = dims - 1
    data = uniform_grid(dims, (0.0, 0.0, 0.0), 1.0 / cells, dtype=dtype, device=device)
    out = tack.empty((data.num_points,), dtype, device=data.device)
    _tangle(out, int(dims[0]), int(dims[1]), *map(float, cells))
    data.fields["tangle"] = Field(H1(data), out)
    return data

tack.dataset.amr_gaussian_pulse

amr_gaussian_pulse(
    dimension=3,
    root_spacing=0.5,
    refinement_ratio=2,
    pulse_origin=(0.0, 0.0, 0.0),
    pulse_width=(0.5, 0.5, 0.5),
    amplitude=0.0001,
    dtype=tack.f32,
    device=None,
)

A two-level AMR hierarchy with a Gaussian pulse, as VTK's vtkAMRGaussianPulseSource makes it: a root image of 6 x 6 (x 6) points from (-2, -2, 0), and two boxes refining parts of it by refinement_ratio.

Each box has VTK's two cell fields: Gaussian-Pulse, amplitude * exp(-sum(((c - pulse_origin) / pulse_width)^2)) over the first dimension axes at each cell's centre c, and Centroid, the centres. An AMRDataSet, its refined cells and points marked.

Source code in packages/tack-vis/src/tack/dataset/sources.py

def amr_gaussian_pulse(dimension=3, root_spacing=0.5, refinement_ratio=2,
                       pulse_origin=(0.0, 0.0, 0.0), pulse_width=(0.5, 0.5, 0.5),
                       amplitude=0.0001, dtype=tack.f32, device=None):
    """A two-level AMR hierarchy with a Gaussian pulse, as VTK's
    ``vtkAMRGaussianPulseSource`` makes it: a root image of 6 x 6 (x 6) points from
    ``(-2, -2, 0)``, and two boxes refining parts of it by ``refinement_ratio``.

    Each box has VTK's two cell fields: ``Gaussian-Pulse``, ``amplitude *
    exp(-sum(((c - pulse_origin) / pulse_width)^2))`` over the first
    ``dimension`` axes at each cell's centre ``c``, and ``Centroid``, the
    centres. An ``AMRDataSet``, its refined cells and points marked."""
    from tack.dataset.amr import AMRDataSet
    from tack.dataset.spaces import Constant

    if dimension not in (2, 3):
        raise ValueError(f"dimension is 2 or 3, not {dimension!r}")
    r = int(refinement_ratio)
    h = float(root_spacing)
    origin = np.array([-2.0, -2.0, 0.0])
    # VTK's patches, as point ranges of the root: x, y (and z) inclusive.
    patches = ([(0, 2), (0, 3), (0, 5)], [(3, 5), (2, 5), (0, 5)])

    def box(points, spacing, corner, offset):
        dims = tuple(points)[:dimension]
        grid = uniform_grid(dims, origin=tuple(corner)[:dimension],
                            spacing=(spacing,) * dimension, dtype=dtype, device=device)
        topology = type(grid.topology)(dims, id_dtype=grid.topology.id_dtype,
                                       device=device, offset=tuple(offset)[:dimension])
        centres = _cell_centres(dims, corner[:dimension], spacing)
        d = (centres - np.asarray(pulse_origin, np.float64)[:dimension]) / \
            np.asarray(pulse_width, np.float64)[:dimension]
        pulse = float(amplitude) * np.exp(-(d * d).sum(axis=1))
        full = np.zeros((len(centres), 3))
        full[:, :dimension] = centres
        from tack.dataset.dataset import DataSet

        out = DataSet(topology, grid.geometry.values, dtype=dtype)
        with out.device:
            out.fields["Gaussian-Pulse"] = Field(Constant(out), tack.array(
                pulse.astype(dtype.numpy_dtype)))
            out.fields["Centroid"] = Field(Constant(out), tack.array(
                full.astype(dtype.numpy_dtype), tack.vec3(dtype)))
        return out

    root = box((6, 6, 6), h, origin, (0, 0, 0))
    fine = []
    for patch in patches:
        lo = np.array([a for a, _ in patch])
        points = [r * (b - a) + 1 for a, b in patch]
        fine.append(box(points, h / r, origin + lo * h, lo * r))
    return AMRDataSet([[root], fine], ratios=[r])

Data in pieces and hierarchies

tack.dataset.PartitionedDataSet

PartitionedDataSet(partitions=(), world=None)

A dataset in pieces: an ordered list of DataSets (partitions), each with its own topology, fields, sets and device, as VTK's vtkPartitionedDataSet. A piece may have no cells. Pieces may be on different devices; an algorithm given the partitioned dataset runs on each piece where the piece is.

It is a container, not a dataset: iterate over it (for piece in pds), index it (pds[i]) or take its length.

Across MPI ranks, it holds this rank's pieces, and world (a tack.communication.World; this process alone by default) is where the others are, for the algorithms that combine pieces.

Source code in packages/tack-vis/src/tack/dataset/partitioned.py

def __init__(self, partitions=(), world=None):
    partitions = tuple(partitions)
    for i, piece in enumerate(partitions):
        if not isinstance(piece, DataSet):
            raise TypeError(f"partition {i} is a {type(piece).__name__}, not a "
                            "tack.dataset.DataSet")
    self._partitions = partitions
    if world is None:
        from tack.communication import World
        world = World()
    self.world = world

partitions property

partitions

The pieces, a tuple of DataSets.

tack.dataset.AMRDataSet

AMRDataSet(levels, ratios=(), domain=None, world=None)

Levels of structured boxes, each level finer than the one before: a block-structured AMR hierarchy, as VTK's vtkOverlappingAMR and AMReX's.

levels holds this rank's boxes of each level, a PartitionedDataSet or a list of DataSets each. Every box has a StructuredTopology (all of one dimension) whose offset is its first cell's index in its level, and any structured geometry: an image, a rectilinear grid or a curvilinear grid. ratios[l] refines level l into l + 1: a whole number, or one per axis. domain is level 0's cell-index box (lo, hi), hi exclusive; by default its boxes' bounding box. world is where the other ranks' boxes are (tack.communication.World; the levels' own, else this process).

Refused when made: boxes of a level that overlap; a box of level l + 1 that, coarsened, is not within level l's boxes (proper nesting); image boxes of a level not on one lattice (one spacing, one origin, axis-aligned: VTK's model); boxes with ghost cells (DUPLICATECELL).

The boxes it holds are copies of those given with their marks computed (REFINEDCELL, REFINEDPOINT, siblings' DUPLICATEPOINT; other bits kept), sharing their topology, geometry and fields. Iterating over it gives (level, box) for this rank's boxes, coarsest level first.

Source code in packages/tack-vis/src/tack/dataset/amr.py

def __init__(self, levels, ratios=(), domain=None, world=None):
    levels = list(levels)
    if not levels:
        raise ValueError("an AMRDataSet has at least one level")
    if world is None:
        world = next((lvl.world for lvl in levels if isinstance(lvl, PartitionedDataSet)),
                     None)
    if world is None:
        from tack.communication import World
        world = World()
    self.world = world
    given = [list(lvl) for lvl in levels]
    self.dimension = self._dimension(given)
    self.ratios = self._ratios(ratios, len(given))
    boxes, self._first = [], []
    for level, pieces in enumerate(given):
        mine = np.array([_box(piece, self.dimension) for piece in pieces],
                        np.int64).reshape(-1, 2, 3)
        every = world._allgather(mine)
        self._first.append(int(sum(len(part) for part in every[:world.rank])))
        boxes.append(np.concatenate(every).reshape(-1, 2, 3))
    self._boxes = boxes
    #: Each box's arrays' cell-index box: the box itself, or grown by ghost
    #: layers (``ghost_cells``).
    self._extents = [b.copy() for b in boxes]
    self.layers = 0
    self._fills = None
    self.domain = self._domain(domain)
    self._check(given)
    self.levels = tuple(PartitionedDataSet([self._marked(level, k, piece)
                                            for k, piece in enumerate(pieces)], world)
                        for level, pieces in enumerate(given))

level_domain

level_domain(level)

domain refined to level: its cell-index box there.

Source code in packages/tack-vis/src/tack/dataset/amr.py

def level_domain(self, level):
    """``domain`` refined to ``level``: its cell-index box there."""
    scale = np.ones(3, np.int64)
    for r in self.ratios[:level]:
        scale *= np.asarray(r, np.int64)
    lo, hi = (np.asarray(c, np.int64) for c in self.domain)
    return lo * scale, hi * scale

boxes

boxes(level)

Every rank's boxes of level, an (n, 2, 3) integer array of cell-index boxes [lo, hi), in order of global index: rank 0's first, each rank's in the order given. Axes a box does not have are [0, 1).

Source code in packages/tack-vis/src/tack/dataset/amr.py

def boxes(self, level):
    """Every rank's boxes of ``level``, an ``(n, 2, 3)`` integer array of cell-index
    boxes ``[lo, hi)``, in order of global index: rank 0's first, each rank's
    in the order given. Axes a box does not have are ``[0, 1)``."""
    return self._boxes[level].copy()

tack.dataset.OctreeDataSet

OctreeDataSet(
    roots,
    tree,
    level,
    key,
    fields=None,
    mask=None,
    world=None,
    units=None,
)

A grid of octrees, Morton-sorted: VTK's vtkHyperTreeGrid with branching factor 2, as a linear octree.

roots is the grid of root cells, a dataset of a StructuredTopology over an image's or a rectilinear grid's coordinates, 2D or 3D (uniform_grid, rectilinear_grid); its cell t (x fastest) is tree t. tree, level and key give each node: its tree, its level (0 at the root) and its Morton key at that level (integer coordinates in the tree interleaved, x lowest). Every node's ancestors must be among the nodes. fields maps names to one value per node (cell data), in the order the nodes are given; mask (one per node, nonzero hides) becomes HIDDENCELL. The nodes are sorted in pre-order and the fields with them. world is where other ranks' trees are, for reductions: each rank holds whole trees. units maps field names to their units (tack.dataset.units); the spatial unit is the roots' (unit).

OctreeDataSet.from_cells makes one from cells' levels and positions, as a RAMSES reader gives them; tack.interop.vtk converts vtkHyperTreeGrids.

Source code in packages/tack-vis/src/tack/dataset/octree.py

def __init__(self, roots, tree, level, key, fields=None, mask=None, world=None,
             units=None):
    from tack.dataset.topology import StructuredTopology

    if not isinstance(roots.topology, StructuredTopology) or \
            not isinstance(roots.geometry.values, (UniformCoordinates, CartesianProduct)):
        raise TypeError("an octree's roots are an image or a rectilinear grid")
    self.roots = roots
    self.dimension = len(roots.topology.point_dims)
    if self.dimension not in (2, 3):
        raise ValueError("an octree is 2D or 3D")
    if world is None:
        from tack.communication import World
        world = World()
    self.world = world
    tree = np.asarray(arrays.to_host(tree) if isinstance(tree, tack.Array) else tree,
                      np.int64).reshape(-1)
    level = np.asarray(arrays.to_host(level) if isinstance(level, tack.Array) else level,
                       np.int64).reshape(-1)
    key = np.asarray(arrays.to_host(key) if isinstance(key, tack.Array) else key,
                     np.uint64).reshape(-1)
    if not len(tree) == len(level) == len(key):
        raise ValueError("tree, level and key give each node: one of each per node")
    d = self.dimension
    if len(level) and (level.min() < 0 or level.max() > MAX_LEVEL[d]):
        raise ValueError(f"levels run from 0 to {MAX_LEVEL[d]} in {d}D")
    if len(tree) and (tree.min() < 0 or tree.max() >= roots.num_cells):
        raise ValueError(f"trees are the roots' cells, 0 to {roots.num_cells - 1}")
    if len(key) and (key >> (d * level).astype(np.uint64)).any():
        raise ValueError("a key holds dimension bits a level: one is too large for its "
                         "level")
    order = _preorder(tree, level, key, d)
    self._tree, self._level, self._key = tree[order], level[order], key[order]
    _check_ancestors(self._tree, self._level, self._key, d)
    leaf = _leaves(self._tree, self._level, self._key, d)
    marks = np.where(leaf, 0, Ghosts.REFINEDCELL).astype(np.uint8)
    if mask is not None:
        mask = np.asarray(arrays.to_host(mask) if isinstance(mask, tack.Array) else mask
                          ).reshape(-1)
        marks[np.asarray(mask[order] != 0)] |= Ghosts.HIDDENCELL
    device = roots.device
    with device:
        self.tree = tack.array(self._tree, tack.i64)
        self.level = tack.array(self._level.astype(np.int32), tack.i32)
        self.key = tack.array(self._key, tack.u64)
        self.marks = tack.array(marks, tack.u8)
    self.fields = {}
    for name, values in (fields or {}).items():
        self.fields[name] = _reordered(values, order, device, len(order), name)
    unknown = set(units or {}) - set(self.fields)
    if unknown:
        raise KeyError(f"units for fields not given: {', '.join(map(repr, sorted(unknown)))}")
    self.units = {name: _units.canonical(unit) for name, unit in (units or {}).items()}
    self._geometry = None
    self._carried_keys = None

num_trees property

num_trees

The trees this rank holds.

unit property

unit

The spatial unit, the roots' (DataSet.unit).

time property

time

The step's time, the roots' (DataSet.time).

time_unit property

time_unit

The unit of time, the roots'.

is_leaf property

is_leaf

Which nodes are leaves: a bool host array, in the nodes' order.

finest property

finest

The deepest level of this rank's nodes (0 with none).

geometry

geometry()

Each node's lower corner and sides, relative to origin: two arrays of 3-vectors (f64 where the device has it, else the roots' type), from the keys, made once.

Source code in packages/tack-vis/src/tack/dataset/octree.py

def geometry(self):
    """Each node's lower corner and sides, relative to ``origin``: two arrays of
    3-vectors (f64 where the device has it, else the roots' type), from the
    keys, made once."""
    if self._geometry is None:
        dtype = _compute_type(self.device, self.roots.dtype)
        lo, size = _root_cells(self.roots)
        with self.device:
            root_lo = tack.array(lo.astype(dtype.numpy_dtype), tack.vec3(dtype))
            root_size = tack.array(size.astype(dtype.numpy_dtype), tack.vec3(dtype))
            corner = tack.empty((self.num_nodes,), tack.vec3(dtype))
            sides = tack.empty((self.num_nodes,), tack.vec3(dtype))
            if self.num_nodes:
                _node_geometry(self.tree, self.level, self.key, root_lo, root_size,
                               self.dimension, 1.0, corner, sides)
        self._geometry = (corner, sides)
    return self._geometry

centers

centers()

Each node's centre, relative to origin: a host (n, 3) array.

Source code in packages/tack-vis/src/tack/dataset/octree.py

def centers(self):
    """Each node's centre, relative to ``origin``: a host ``(n, 3)`` array."""
    corner, sides = self.geometry()
    return corner.numpy().reshape(-1, 3) + 0.5 * sides.numpy().reshape(-1, 3)

measures

measures()

Each node's volume (area in 2D): a host array.

Source code in packages/tack-vis/src/tack/dataset/octree.py

def measures(self):
    """Each node's volume (area in 2D): a host array."""
    _, sides = self.geometry()
    s = sides.numpy().reshape(-1, 3)
    return s[:, 0] * s[:, 1] * (s[:, 2] if self.dimension == 3 else 1.0)

leaves

leaves()

This rank's leaves as a DataSet of voxels (pixels in 2D), relative to origin: cell c is the c-th leaf in the nodes' order (np.flatnonzero(octree.is_leaf)[c]), with every field's value there as cell data, and hidden leaves HIDDENCELL. Corners that coincide are one point, found by sorting their indices on the finest level's lattice; hanging corners stay hanging, as the tree has them. Every filter that reads cell data works on it (threshold, cell_centers, probe).

Each rank gives its own trees' leaves; a corner on a face between two ranks' trees is a point on both, unmarked.

Source code in packages/tack-vis/src/tack/dataset/octree.py

def leaves(self):
    """This rank's leaves as a ``DataSet`` of voxels (pixels in 2D), relative to
    ``origin``: cell ``c`` is the ``c``-th leaf in the nodes' order
    (``np.flatnonzero(octree.is_leaf)[c]``), with every field's value there as
    cell data, and hidden leaves ``HIDDENCELL``. Corners that coincide are one
    point, found by sorting their indices on the finest level's lattice;
    hanging corners stay hanging, as the tree has them. Every filter that
    reads cell data works on it (``threshold``, ``cell_centers``, ``probe``).

    Each rank gives its own trees' leaves; a corner on a face between two
    ranks' trees is a point on both, unmarked."""
    from tack.algorithms import argsort, gather, inclusive_scan
    from tack.dataset import ids, shapes
    from tack.dataset.dataset import DataSet, Field
    from tack.dataset.spaces import Constant
    from tack.dataset.topology import UnstructuredTopology, uniform_cells

    d = self.dimension
    is_leaf = self.is_leaf
    count = int(is_leaf.sum())
    corners = count << d
    finest = self.finest
    nx, ny, nz = ([n - 1 for n in self.roots.topology.point_dims] + [1])[:3]
    idt = ids.choose(None, corners)
    dtype = _compute_type(self.device, self.roots.dtype)
    with self.device:
        leaf = tack.array(np.flatnonzero(is_leaf).astype(np.int64), tack.i64)
        columns = [tack.empty((corners,), tack.i64) for _ in range(3)]
        if count:
            _leaf_corners(leaf, self.tree, self.level, self.key, nx, ny, d, finest,
                          *columns)
        # Sorted by z, then y, then x (stable sorts from the fastest axis), so
        # points come x fastest; the runs of equal corners are the points.
        order = argsort(columns[0], index_dtype=tack.i64) if corners else \
            tack.empty((0,), tack.i64)
        for column in columns[1:d]:
            if corners:
                order = gather(order, argsort(gather(column, order), index_dtype=tack.i64))
        si, sj, sk = (gather(c, order) for c in columns)
        flags = tack.empty((corners,), tack.i32)
        numbers = tack.empty((corners,), tack.i64)
        connectivity = tack.empty((corners,), idt)
        points = 0
        if corners:
            _run_starts(si, sj, sk, flags)
            points = int(inclusive_scan(flags, numbers, corners))
            _numbered(order, numbers, connectivity)
        axes, plane = self._axes(dtype)
        positions = tack.empty((points,), tack.vec3(self.roots.dtype))
        if points:
            _lattice_points(si, sj, sk, flags, numbers, *axes, nx, ny, nz, finest, d,
                            plane, 1.0, positions)
        types, offsets = uniform_cells(shapes.VOXEL if d == 3 else shapes.PIXEL, count,
                                       idt)
        topology = UnstructuredTopology(types, offsets, connectivity, num_points=points)
        fields = {name: Field(Constant(topology), gather(values, leaf),
                              self.units.get(name))
                  for name, values in self.fields.items()}
        marks = gather(self.marks, leaf)
        ghosts = None
        if count and (marks.numpy() & Ghosts.HIDDENCELL).any():
            hidden = marks.numpy() & Ghosts.HIDDENCELL
            ghosts = Ghosts(cells=tack.array(hidden.astype(np.uint8), tack.u8))
    return DataSet(topology, positions, fields=fields, dtype=self.roots.dtype,
                   ghosts=ghosts, **self._frame())

from_cells classmethod

from_cells(
    roots,
    level,
    position,
    fields=None,
    world=None,
    units=None,
)

An octree of cells given by their level and position (an (n, 3) world position inside each cell, its centre say), as a RAMSES reader or yt gives them -- whole octs: every refined cell's children. Each cell's tree is the root cell holding the position; its key, the integer coordinates there at its level. Ancestors no cell gives are added, their fields NaN (integer fields 0): reductions of the leaves never read them.

Source code in packages/tack-vis/src/tack/dataset/octree.py

@classmethod
def from_cells(cls, roots, level, position, fields=None, world=None, units=None):
    """An octree of cells given by their ``level`` and ``position`` (an ``(n, 3)``
    world position inside each cell, its centre say), as a RAMSES reader or
    yt gives them -- whole octs: every refined cell's children. Each cell's
    tree is the root cell holding the position; its key, the integer
    coordinates there at its level. Ancestors no cell gives are added, their
    fields NaN (integer fields 0): reductions of the leaves never read
    them."""
    level = np.asarray(level, np.int64).reshape(-1)
    position = np.asarray(position, np.float64).reshape(-1, 3) - np.asarray(roots.origin)
    d = len(roots.topology.point_dims)
    axes = _root_axes(roots)
    tree = np.zeros(len(level), np.int64)
    stride = 1
    index = np.zeros((len(level), 3), np.int64)
    for a in range(d):
        edges = axes[a]
        cell = np.clip(np.searchsorted(edges, position[:, a], side="right") - 1, 0,
                       len(edges) - 2)
        tree += cell * stride
        stride *= len(edges) - 1
        fraction = (position[:, a] - edges[cell]) / (edges[cell + 1] - edges[cell])
        index[:, a] = np.clip(np.floor(fraction * (1 << level)), 0,
                              (1 << level) - 1).astype(np.int64)
    key = _interleave(index, d)
    tree, level, key, fields = _with_ancestors(tree, level, key, fields or {}, d)
    return cls(roots, tree, level, key, fields=fields, world=world, units=units)

Topologies

tack.dataset.StructuredTopology

StructuredTopology(
    point_dims, id_dtype=None, device=None, offset=None
)

Bases: _Topology

The cells of a structured grid of point_dims points, x fastest: lines, quads or hexahedra in vtkStructuredGrid's order, iterated by (i, j, k).

offset is the grid index of its first point in a larger grid, one per axis (VTK's extent's start; zeros by default): a piece split cut, or an image whose extent starts elsewhere. What depends on a cell's place in the whole grid -- tetrahedralize's alternating splits -- reads it, so neighbouring pieces agree.

Source code in packages/tack-vis/src/tack/dataset/topology.py

def __init__(self, point_dims, id_dtype=None, device=None, offset=None):
    point_dims = tuple(int(d) for d in point_dims)
    if not 1 <= len(point_dims) <= 3 or min(point_dims) < 2:
        raise ValueError(f"point_dims must be 1 to 3 sizes of at least 2, not {point_dims}")
    offset = (0,) * len(point_dims) if offset is None else tuple(int(o) for o in offset)
    if len(offset) != len(point_dims):
        raise ValueError(f"offset needs one index per axis ({len(point_dims)}), not "
                         f"{offset}")
    self.point_dims = point_dims
    self.offset = offset
    self.num_cells = int(np.prod([d - 1 for d in point_dims]))
    self.num_points = int(np.prod(point_dims))
    self.id_dtype = ids.choose(id_dtype, max(self.num_points, 12 * self.num_cells))
    kind, shape = self._KINDS[len(point_dims)]
    dims = point_dims + (1,) * (3 - len(point_dims))
    # One group, kept: derived incidence is laid out per group.
    self._groups = [DomainGroup(kind, shape, (*dims, self.num_cells), self.num_cells, 0)]
    # It holds no arrays, so its device is the one it is made for: ``device``,
    # or the default device when it is made.
    self._found_device = backend_of(device)

tack.dataset.UnstructuredTopology

UnstructuredTopology(
    types,
    offsets,
    connectivity,
    num_points=None,
    id_dtype=None,
)

Bases: _Topology

Cells of any of the linear shapes: a VTK type per cell (u8), and each cell's point ids at connectivity[offsets[c]:offsets[c + 1]], over num_points points -- by default one past the highest id used, but a dataset may have points no cell uses, as VTK's may.

offsets and connectivity are fields or arrays of ids, kept in id_dtype: i32 unless asked for tack.i64, given i64 fields, or too large for i32 -- a mesh whose ids, connectivity or derived arrays (twelve edges a cell) pass 2**31 - 1.

Any of the three may be an implicit array (tack.arrays), kept as it is: the types of cells all of one shape a ConstantArray, the offsets of cells all of one size a CountingArray (uniform_cells makes both), the connectivity of cells that each have their own points a CountingArray too. VTK's constant and affine arrays map onto them.

Source code in packages/tack-vis/src/tack/dataset/topology.py

@on_inputs_device
def __init__(self, types, offsets, connectivity, num_points=None, id_dtype=None):
    self.types = _as_field(types, tack.u8)
    self.num_cells = size_of(self.types)
    offsets, connectivity = (v if isinstance(v, _Implicit) else ids.host_or_field(v)
                             for v in (offsets, connectivity))
    if size_of(offsets) != self.num_cells + 1:
        raise ValueError(f"offsets must have num_cells + 1 = {self.num_cells + 1} "
                         f"entries, not {size_of(offsets)}")
    if num_points is None:
        used = (connectivity.numpy() if isinstance(connectivity, (Array, _Implicit))
                else connectivity)
        num_points = int(used.max()) + 1 if used.size else 0
    self.num_points = int(num_points)
    extent = max(self.num_points, size_of(connectivity), 12 * self.num_cells)
    if id_dtype is None and any(isinstance(v, _Implicit) and v.dtype == tack.i64
                                for v in (offsets, connectivity)):
        id_dtype = tack.i64
    self.id_dtype = ids.choose(id_dtype, extent,
                               given=[v for v in (offsets, connectivity)
                                      if isinstance(v, Array)])
    self.offsets = _ids_of(offsets, self.id_dtype)
    self.connectivity = _ids_of(connectivity, self.id_dtype)
    self._groups = None
    if _checks_ids():
        _refuse("unstructured topology", self._id_problems())

tack.dataset.polyhedra

Polyhedral topologies: every cell a list of faces, every face a polygon stored once.

docs/vis-guide/design/polyhedra.md, phase 1. A PolyhedralTopology holds

• face_offsets/face_points: each face's points (CSR), in the order that winds the face out of its side 0 -- the cell its normal points away from;
• cell_offsets/cell_faces: each cell's faces (CSR);
• cell_face_sides: one u8 per cell -> face entry, 0 if the face is wound out of that cell (it is the face's side 0), 1 if into it. This is the shape path's side_slot under the same meaning.

Face numbering is the input's: a producer's face ids (an MPAS edge id) are the face ids, so face data lines up with the file, and nothing re-derives them. Everything else is derived on demand by sorting, and kept:

• faces(): the shape path's two-sided record per face, [cell0, local0, cell1, local1] (-1, -1 for a missing side 1), and the boundary;
• cell_points(): each cell's points, its face points sorted and unique -- exactly what VTK keeps as a polyhedron's connectivity;
• edges(): the faces' consecutive point pairs as (low, high) keys, sorted into runs -- the same numbering the shape path gives the same edges.

Ids are the topology's id_dtype (tack.dataset.ids), as on the shape path.

Orientation comes from the input, never from geometry. check_winding lists the cells that do not walk each of their edges once in each direction; orient turns a face wound into its only cell around. as_polyhedra converts any shape-based dataset, reusing its derived faces as they stand.

tack.dataset.PolyhedralTopology

PolyhedralTopology(
    face_offsets,
    face_points,
    cell_offsets,
    cell_faces,
    cell_face_sides,
    num_points=None,
    id_dtype=None,
)

Bases: _Topology

Cells given by their faces, faces stored once (docs/vis-guide/design/polyhedra.md).

face_offsets/face_points (CSR) are each face's points, in the order that winds it out of its side 0; cell_offsets/cell_faces (CSR) each cell's faces; cell_face_sides (u8, one per cell_faces entry) 0 where the face is wound out of that cell and 1 where into it. num_points defaults to one past the highest point id used. Ids are kept in id_dtype, as UnstructuredTopology keeps them.

Every cell is a Polyhedron: no reference element, so no spaces with a basis (L2, H1 above order 1) and no parametric coordinates.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def __init__(self, face_offsets, face_points, cell_offsets, cell_faces, cell_face_sides,
             num_points=None, id_dtype=None):
    face_offsets, face_points, cell_offsets, cell_faces = (
        _given(a) for a in (face_offsets, face_points, cell_offsets, cell_faces))
    if num_points is None:
        num_points = _points_used(face_points)
    # Cell points and directed edges hold each face's points once per side.
    extent = max(int(num_points), 2 * size_of(face_points),
                 2 * size_of(face_offsets), size_of(cell_faces))
    given = (face_offsets, face_points, cell_offsets, cell_faces)
    self.id_dtype = idt = ids.choose(id_dtype, extent, given)
    self.face_offsets = _ids_of(face_offsets, idt)
    self.face_points = _ids_of(face_points, idt)
    self.cell_offsets = _ids_of(cell_offsets, idt)
    self.cell_faces = _ids_of(cell_faces, idt)
    self.cell_face_sides = _as_field(cell_face_sides, tack.u8)
    self.num_faces = size_of(self.face_offsets) - 1
    self.num_cells = size_of(self.cell_offsets) - 1
    if self.num_faces < 0 or self.num_cells < 0:
        raise ValueError("offsets have one entry more than there are faces or cells")
    if size_of(self.cell_face_sides) != size_of(self.cell_faces):
        raise ValueError("cell_face_sides needs one entry per cell_faces entry")
    self.num_points = int(num_points)
    self._cell_points = None
    self._groups = None
    if _checks_ids():
        _refuse("polyhedral topology", self._id_problems())

faces

faces()

The faces' sides and boundary (PolygonFaces), derived on first use.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def faces(self):
    """The faces' sides and boundary (``PolygonFaces``), derived on first use."""
    if self._faces is None:
        self._faces = PolygonFaces(self)
    return self._faces

edges

edges()

The faces' edges (PolygonEdges), derived on first use.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def edges(self):
    """The faces' edges (``PolygonEdges``), derived on first use."""
    if self._edges is None:
        self._edges = PolygonEdges(self)
    return self._edges

point_links

point_links()

Each point's entries in cell_points(), derived on first use and kept: (offsets, entries), CSR by point, entries in increasing order (so in cell order).

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

def point_links(self):
    """Each point's entries in ``cell_points()``, derived on first use and kept:
    ``(offsets, entries)``, CSR by point, entries in increasing order (so in
    cell order)."""
    if self._point_links is None:
        _, point_ids = self.cell_points()
        order, offsets = bucket_order(point_ids, self.num_points, self.id_dtype)
        self._point_links = (offsets, order)
    return self._point_links

cell_points

cell_points()

Each cell's points, sorted and unique: (offsets, point_ids), CSR.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def cell_points(self):
    """Each cell's points, sorted and unique: ``(offsets, point_ids)``, CSR."""
    if self._cell_points is None:
        self._cell_points = _derive_cell_points(self)
    return self._cell_points

groups

groups()

The cells as one Polyhedron group.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def groups(self):
    """The cells as one ``Polyhedron`` group."""
    if self._groups is None:
        point_offsets, point_ids = self.cell_points()
        self._groups = [DomainGroup(
            _PolyhedralCells, self.cell_shape,
            (self.cell_offsets, self.cell_faces, self.cell_face_sides, self.face_offsets,
             self.face_points, point_offsets, point_ids, self.num_cells),
            self.num_cells, 0)]
    return self._groups

from_cell_faces classmethod

from_cell_faces(
    cells,
    num_points,
    positions=None,
    orient=False,
    id_dtype=None,
)

A topology from each cell's own list of faces (point-id rings), as VTK and many readers give them: copies of a face are matched by point set, the first copy found is the face -- its cell side 0 -- and a second copy, run the other way, is side 1. Two copies run the same way are refused.

With orient, the copies are made consistent first, on the host: within each cell, faces walk each shared edge in opposite directions; across cells, the two copies of a shared face run opposite ways (whole cells are reversed to make it so); then each connected component is turned outward as a whole, by the sign of its total volume -- the one geometric test, made once per component on a sum over many cells, which thin cells cannot fool as they fool a per-cell test. positions (one row per point) are needed for that last step.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@classmethod
def from_cell_faces(cls, cells, num_points, positions=None, orient=False, id_dtype=None):
    """A topology from each cell's own list of faces (point-id rings), as VTK and many
    readers give them: copies of a face are matched by point set, the first copy
    found is the face -- its cell side 0 -- and a second copy, run the other way,
    is side 1. Two copies run the same way are refused.

    With ``orient``, the copies are made consistent first, on the host:
    within each cell, faces walk each shared edge in opposite directions;
    across cells, the two copies of a shared face run opposite ways (whole
    cells are reversed to make it so); then each connected component is
    turned outward as a whole, by the sign of its total volume -- the one
    geometric test, made once per component on a sum over many cells, which
    thin cells cannot fool as they fool a per-cell test. ``positions`` (one
    row per point) are needed for that last step.
    """
    cells = [[[int(p) for p in ring] for ring in faces] for faces in cells]
    if orient:
        if positions is None:
            raise ValueError("orient needs the points' positions, for the sign of each "
                             "component's volume")
        cells = _orient_cells(cells, np.asarray(positions, float))
    known, face_points, face_offsets = {}, [], [0]
    cell_offsets, cell_faces, sides = [0], [], []
    same_way = 0
    for faces in cells:
        for ring in faces:
            key = tuple(sorted(ring))
            f = known.get(key)
            if f is None:
                f = known[key] = len(face_offsets) - 1
                face_points.extend(ring)
                face_offsets.append(len(face_points))
                side = 0
            else:
                if not _opposite(face_points[face_offsets[f]:face_offsets[f + 1]], ring):
                    same_way += 1
                side = 1
            cell_faces.append(f)
            sides.append(side)
        cell_offsets.append(len(cell_faces))
    if same_way:
        raise ValueError(f"{same_way} shared faces are wound the same way by both their "
                         "cells; from_cell_faces(..., orient=True) makes them consistent")
    return cls(np.array(face_offsets, np.int64), np.array(face_points, np.int64),
               np.array(cell_offsets, np.int64), np.array(cell_faces, np.int64),
               np.array(sides, np.uint8), num_points=num_points, id_dtype=id_dtype)

arrays

arrays()

The five defining arrays as host arrays, in constructor order.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

def arrays(self):
    """The five defining arrays as host arrays, in constructor order."""
    return (self.face_offsets.numpy(), self.face_points.numpy(),
            self.cell_offsets.numpy(), self.cell_faces.numpy(),
            self.cell_face_sides.numpy())

tack.dataset.PolygonalTopology

PolygonalTopology(
    loop_offsets,
    loop_points,
    num_points=None,
    id_dtype=None,
)

Bases: PolyhedralTopology

Polygons of any number of points, the polyhedral topology one dimension down: cells are polygons, their facets edges.

Built from point loops -- loop_points[loop_offsets[c]:loop_offsets[c + 1]] is polygon c's points in order -- as producers and isosurfaces give them. Facet k of a polygon is its edge from loop point k to the next. Edges are matched by sorting: the first polygon to use an edge is its side 0 and the edge runs its way; a second must run it the other way, or the winding is inconsistent (or the surface not orientable) and is refused, as is an edge of three polygons. loops() gives the loops back.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def __init__(self, loop_offsets, loop_points, num_points=None, id_dtype=None):
    loop_offsets, loop_points = _given(loop_offsets), _given(loop_points)
    total = size_of(loop_points)
    if num_points is None:
        num_points = _points_used(loop_points)
    # Its facets' points are two per edge, and there are as many edges as loop points.
    idt = ids.choose(id_dtype, max(int(num_points), 4 * total + 2),
                     given=(loop_offsets, loop_points))
    loop_offsets = _ids_of(loop_offsets, idt)
    loop_points = _ids_of(loop_points, idt)
    if _checks_ids():
        # Before the facets are derived from the loops, which index by point.
        _refuse("polygonal topology", _csr_problems(
            loop_offsets, loop_points, int(num_points), "loop_offsets", "loop_points",
            "num_points"))
    n = size_of(loop_offsets) - 1
    keys = tack.empty((total,), tack.vec(2, idt))
    directions = tack.empty((total,), tack.i32)
    cell_faces = tack.empty((total,), idt)
    sides = tack.empty((total,), tack.u8)
    count = 0
    if total:
        _loop_edge_keys(loop_offsets, loop_points, keys, directions, n)
        order = key_order(keys, num_points, idt)
        offsets, count = run_offsets(keys, order)
    face_points = tack.empty((2 * count,), idt)
    if count:
        problems = tack.zeros((2,), tack.i32)
        _facets_from_runs(order, offsets, keys, directions, face_points, cell_faces,
                          sides, problems, count)
        if problems[0]:
            raise ValueError(f"{problems[0]} edges are shared by more than two polygons")
        if problems[1]:
            raise ValueError(f"{problems[1]} edges are walked the same way by both their "
                             "polygons: the winding is inconsistent, or the surface is "
                             "not orientable")
    face_offsets = tack.empty((count + 1,), idt)
    _every_other(face_offsets)
    super().__init__(face_offsets, face_points, loop_offsets, cell_faces, sides,
                     num_points=num_points, id_dtype=idt)
    self.loop_offsets = loop_offsets
    self.loop_points = loop_points

loops

loops()

Each polygon's points in order: (loop_offsets, loop_points).

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

def loops(self):
    """Each polygon's points in order: ``(loop_offsets, loop_points)``."""
    return self.loop_offsets, self.loop_points

tack.dataset.PolylineTopology

PolylineTopology(
    line_offsets,
    line_points,
    num_points=None,
    id_dtype=None,
)

Bases: _Topology

Polylines of any number of points, the polyhedral topology two dimensions down: each cell a chain of points, line_points[line_offsets[c]:line_offsets[c + 1]], two or more, as streamlines and VTK's VTK_POLY_LINE give them.

Its facets (faces(), codimension 1) are the chains' ends: a facet per point where chains end, its normal along the chain out of its side 0 -- the first chain ending there, in cell order -- and into its side 1, the second. A vertex has no winding, so any two chains meeting end to end are consistent. Three or more ending at one point (a branch) the faces cannot hold, and faces() refuses them, as an edge of three polygons is refused; the topology itself, and everything that does not need its faces, works. A chain may close (its last point its first): its two ends are then one interior facet. lines() gives the chains back; a cell's points (cell_points()) are its chain's, sorted and unique.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def __init__(self, line_offsets, line_points, num_points=None, id_dtype=None):
    line_offsets, line_points = _given(line_offsets), _given(line_points)
    total = size_of(line_points)
    if num_points is None:
        num_points = _points_used(line_points)
    n = size_of(line_offsets) - 1
    if n < 0:
        raise ValueError("line_offsets has one entry more than there are polylines")
    self.id_dtype = idt = ids.choose(id_dtype, max(int(num_points), total, 2 * n + 1),
                                     given=(line_offsets, line_points))
    self.line_offsets = _ids_of(line_offsets, idt)
    self.line_points = _ids_of(line_points, idt)
    self.num_cells = n
    self.num_points = int(num_points)
    if _checks_ids():
        # Before the ends are read from the lines.
        self.validate()
    self._ends = tack.empty((2 * n,), idt)
    if n:
        short = tack.zeros((1,), tack.i32)
        _line_ends(self.line_offsets, self.line_points, self._ends, short, n)
        if short[0]:
            raise ValueError(f"{short[0]} polylines have fewer than two points")
    self._cell_points = None
    self._groups = None

lines

lines()

Each polyline's points in order: (line_offsets, line_points).

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

def lines(self):
    """Each polyline's points in order: ``(line_offsets, line_points)``."""
    return self.line_offsets, self.line_points

cell_points

cell_points()

Each cell's points, sorted and unique: (offsets, point_ids), CSR.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def cell_points(self):
    """Each cell's points, sorted and unique: ``(offsets, point_ids)``, CSR."""
    if self._cell_points is None:
        # Each chain as its cell's one "face".
        identity = tack.empty((self.num_cells + 1,), self.id_dtype)
        _counting(identity)
        self._cell_points = _derive_cell_points(self, identity, identity, self.line_offsets,
                                                self.line_points)
    return self._cell_points

groups

groups()

The cells as one PolyLine group.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def groups(self):
    """The cells as one ``PolyLine`` group."""
    if self._groups is None:
        point_offsets, point_ids = self.cell_points()
        self._groups = [DomainGroup(
            _PolylineCells, self.cell_shape,
            (self.line_offsets, self.line_points, point_offsets, point_ids, self.num_cells),
            self.num_cells, 0)]
    return self._groups

faces

faces()

The chains' ends as faces (PolygonFaces of one-point faces), derived on first use and kept; refuses a point where more than two chains end.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def faces(self):
    """The chains' ends as faces (``PolygonFaces`` of one-point faces), derived on
    first use and kept; refuses a point where more than two chains end."""
    if self._faces is None:
        n, idt = self.num_cells, self.id_dtype
        cell_faces = tack.empty((2 * n,), idt)
        sides = tack.empty((2 * n,), tack.u8)
        count = 0
        if n:
            order = key_order(self._ends, self.num_points, idt)
            offsets, count = run_offsets(self._ends, order)
        face_points = tack.empty((count,), idt)
        if count:
            crowded = tack.zeros((1,), tack.i32)
            _ends_from_runs(order, offsets, self._ends, face_points, cell_faces, sides,
                            crowded, count)
            if crowded[0]:
                raise ValueError(
                    f"more than two polylines end at {crowded[0]} points: the faces of a "
                    "branching network are not defined (split it at the branches)")
        self.face_offsets = tack.empty((count + 1,), idt)
        _counting(self.face_offsets)
        self.face_points = face_points
        self.cell_offsets = tack.empty((n + 1,), idt)
        _every_other(self.cell_offsets)
        self.cell_faces = cell_faces
        self.cell_face_sides = sides
        self.num_faces = count
        self._faces = PolygonFaces(self)
    return self._faces

edges

edges()

Not derived: a polyline's segments are parts of its cells, and nothing needs them as global entities yet.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def edges(self):
    """Not derived: a polyline's segments are parts of its cells, and nothing
    needs them as global entities yet."""
    raise NotImplementedError("a polyline topology's edges (its segments) are not "
                              "derived")

tack.dataset.SizeBuckets

Bases: _OnItsArrays

Groups a polyhedral topology's cells by size, for kernels that keep per-cell scratch: a cell's key is the first cap in caps at least half its facet points (an upper bound on its edges, so on its crossings of any isovalue). for_each/launch_groups given the buckets split the cells by them, and each launch's view carries its cap as the class constant MAX_SCRATCH, so tack.local_array(dtype, cells.MAX_SCRATCH) has a compile-time size -- as launches split by shape and order. One per topology and caps.

cell_keys

cell_keys()

Each cell's bucket, an i32 array.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def cell_keys(self):
    """Each cell's bucket, an i32 array."""
    if "_keys" not in self.__dict__:
        t = self.topology
        n = t.num_cells
        field = tack.empty((n,), tack.i32)
        if n:
            sizes = tack.empty((n,), tack.i32)
            _cell_sizes(t.cell_offsets, t.cell_faces, t.face_offsets, sizes, n)
            caps = _as_field(np.asarray(self.caps, np.int32), tack.i32)
            largest = tack.zeros((1,), tack.i32)
            _bucket_keys(sizes, caps, field, largest, n)
            if largest[0]:
                raise ValueError(f"a cell needs scratch for {largest[0]} entries, more "
                                 f"than the largest bucket, {self.caps[-1]}")
        self._keys = field
    return self._keys

domain_mixin

domain_mixin(key)

The mixin that gives a launch of bucket key its MAX_SCRATCH.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

def domain_mixin(self, key):
    """The mixin that gives a launch of bucket ``key`` its ``MAX_SCRATCH``."""
    return _bucket_mixin(self.caps[key])

tack.dataset.orient

orient(topology)

A topology whose faces wound into their only cell are turned around: each such face's points reversed and its one reference made side 0. Other faces are left as they are; check_winding judges the rest.

Source code in packages/tack-vis/src/tack/dataset/polyhedra.py

@on_inputs_device
def orient(topology):
    """A topology whose faces wound into their only cell are turned around: each such
    face's points reversed and its one reference made side 0. Other faces are left as
    they are; ``check_winding`` judges the rest."""
    face_offsets, face_points, cell_offsets, cell_faces, sides = topology.arrays()
    counts, bad = _side_counts(topology)
    if bad:
        raise ValueError(f"{bad} cell -> face entries name no face, or a bad side")
    counts = counts.numpy().reshape(-1, 2)
    inward = (counts[:, 0] == 0) & (counts[:, 1] == 1)
    face_points = face_points.copy()
    for f in np.flatnonzero(inward):
        a, b = face_offsets[f], face_offsets[f + 1]
        face_points[a:b] = face_points[a:b][::-1]
    sides = np.where(inward[cell_faces], 0, sides).astype(np.uint8)
    return PolyhedralTopology(face_offsets, face_points, cell_offsets, cell_faces, sides,
                              num_points=topology.num_points, id_dtype=topology.id_dtype)

Cell shapes

tack.dataset.shapes

Cell shapes: VTK's linear cells as template classes.

Each shape is a @tack.device_class class whose @tack.func methods are its parametric functions. A kernel takes a shape object as a template argument and calls its methods directly, so every compiled variant is specialized to one shape and contains no dispatch on the shape::

@tack.kernel
def centers(cell, conn, xyz, out):
    for c in range(out.shape[0]):
        pc = cell.parametric_center()
        x = tack.vec3(0.0, 0.0, 0.0)
        for j in range(cell.NUM_POINTS):
            x += cell.shape_function(j, pc) * xyz[conn[c, j]]
        out[c] = x

centers(shapes.Hexahedron(), conn, xyz, out)

A mesh of several shapes runs such a kernel once per shape present, over that shape's cells.

Every shape carries VTK's cell type id (Hexahedron.ID is 12, as in vtkCellType.h), and VTK's point order, parametric coordinates, edges and faces, so VTK serves as the reference without any remapping:

=========== ==== ====== === ============================================== class ID points dim parametric domain =========== ==== ====== === ============================================== Vertex 1 1 0 the point (0, 0, 0) Line 3 2 1 r in [0, 1] Triangle 5 3 2 r, s >= 0, r + s <= 1 Pixel 8 4 2 r, s in [0, 1]; points in x-fastest order Quad 9 4 2 r, s in [0, 1]; points counterclockwise Tetra 10 4 3 r, s, t >= 0, r + s + t <= 1 Voxel 11 8 3 r, s, t in [0, 1]; points in x-fastest order Hexahedron 12 8 3 r, s, t in [0, 1]; each face counterclockwise Wedge 13 6 3 a triangle in (r, s) times t in [0, 1] Pyramid 14 5 3 r, s, t in [0, 1]; the apex is (0, 0, 1) =========== ==== ====== === ==============================================

Class constants: ID, NUM_POINTS, DIMENSION, NUM_EDGES and NUM_FACES. Faces are codimension-1: a 3D shape's faces are polygons, a 2D shape's are its edges (lines, in the order going round the cell), and a line has none. VTK counts a 2D cell's faces as 0; Tack does not, so a 2D mesh has faces, as a polygonal one does.

Methods, all callable only from kernels and device functions:

• parametric_point(j), parametric_center(): 3-vectors.
• shape_function(j, pc), shape_gradient(j, pc): VTK's InterpolateFunctions and InterpolateDerivs for point j at parametric coordinates pc, a 3-vector; the gradient is (d/dr, d/ds, d/dt).
• is_inside(pc, tol): 1 if pc is in the parametric domain widened by tol, else 0.
• edge_corner(e, k), face_corner(f, k), face_shape(f), face_num_points(f): VTK's GetEdge/GetFace as indices into the cell's points, and the face's shape id.
• interpolate(values, pc), interpolate_point(pts, pc), jacobian(pts, pc), world_to_parametric(pts, x), nearest_point(pts, x, pc): the geometry of a cell whose points are given, defined once for every shape.

Coordinates a shape does not use are ignored: a quad's functions do not read t. A point index must satisfy 0 <= j < NUM_POINTS, and edge and face indices likewise; others give unspecified values.

The geometric methods take the cell's points from a local array the caller gathered, three coordinates per point, and per-point values from another::

pts = tack.local_array(tack.f32, 3 * cell.NUM_POINTS)
for j in range(cell.NUM_POINTS):
    pts[3 * j], pts[3 * j + 1], pts[3 * j + 2] = xyz[conn[c, j]]
pc, found = cell.world_to_parametric(pts, x)
if found and cell.is_inside(pc, 1e-4):
    ...

Spaces

tack.dataset.spaces

Spaces: where a field's values live, on one topology, and the layout that follows.

A space is constructed on a topology (or a dataset, for its topology) and parameters, as MFEM's FiniteElementSpace is on a mesh::

td.H1(data)                 # one value per point, interpolated linearly
td.H1(data, order=2)        # quadratic: values on points, edges, quad faces, cells
td.L2(data)                 # each cell's own corner values: linear DG
td.L2(data, order=orders)   # an order per cell (0 or 1 here): p-adaptive DG
td.Constant(data)           # one value per cell
td.Values(data, "faces")    # one value per face, no basis
td.SideTraces(data)         # values on each side of each face, at its points

Spaces are interned on their topology: the same kind and parameters on the same topology is the same object, so fields that share a space share its layout -- and equality is identity. The topology holds them weakly: a space holds its topology, and a cycle would keep a dropped dataset's arrays (on a GPU, its device memory) until Python's cycle collector ran. The layout is what a kernel needs beyond a field's values to find them:

• H1, order 1: nothing of its own. Its values are the points', found through the topology's connectivity.
• H1, order 2: its numbering of values -- the points', then one per edge (edge_base + edge), one per quad face (face_dofs) and one per hexahedron (cell_dofs) -- built from the derived edges and faces. Its views read a cell's through the cell's edge and face ids.
• L2, order 1: offsets, where each cell's values start. An unstructured topology's are its connectivity offsets, so they are shared, not copied. With an order per cell, the offsets come from a scan of each cell's count of values (1 at order 0, its corners at order 1).
• Constant and Values: nothing; value k belongs to entity k.

Every space knows size, the number of values a field on it holds, and on, the entities a loop must run over to read it through view.

A space whose cells differ -- an order per cell -- varies: it gives each cell a key (cell_keys, here the order), and tack.dataset.for_each launches each shape group once per key, with the view mixin_for(key), so every launch is specialized as a uniform space's is.

tack.dataset.Space

What every space has: topology, size, on and a field-view mixin, and the layout attributes that mixin reads.

attributes

attributes(group=None)

The layout a view of this space reads for group, beside the values.

Source code in packages/tack-vis/src/tack/dataset/spaces.py

def attributes(self, group=None):
    """The layout a view of this space reads for ``group``, beside the values."""
    return {}

extra_mixins

extra_mixins(group)

Mixins a field view of this space needs for group beyond its own and the storage.

Source code in packages/tack-vis/src/tack/dataset/spaces.py

def extra_mixins(self, group):
    """Mixins a field view of this space needs for ``group`` beyond its own and
    the storage."""
    return ()

mixin_for

mixin_for(key, shape=None)

The field-view mixin for cells of key (None for a uniform space) and shape.

Source code in packages/tack-vis/src/tack/dataset/spaces.py

def mixin_for(self, key, shape=None):
    """The field-view mixin for cells of ``key`` (``None`` for a uniform space) and
    ``shape``."""
    return self.mixin

cell_keys

cell_keys()

For a space that varies: an i32 array of each cell's key, by cell id.

Source code in packages/tack-vis/src/tack/dataset/spaces.py

def cell_keys(self):
    """For a space that ``varies``: an i32 array of each cell's key, by cell id."""
    return

tack.dataset.H1

Bases: Space

Continuous: values shared by the cells around them, interpolated by the shape's functions. Order 1 holds one value per point -- today's point data; order 2 adds one per edge, per quad face and per hexahedron or quad, for the shape's quadratic Lagrange functions (tetrahedra, hexahedra, voxels, wedges, triangles, quads and pixels).

tack.dataset.L2

Bases: Space

Discontinuous: each cell's own values, at offsets[cell]. Order 1 holds a value at each of the cell's corners -- a linear DG field -- and order 0 one value per cell. order is one order for every cell, or an array of one per cell (p-adaptive), read when the space is made.

offsets property

offsets

Where each cell's values start: value j of cell c is value offsets[c] + j, and offsets[num_cells] is size.

cell_keys

cell_keys()

Each cell's order, an i32 array (for a space with an order per cell).

Source code in packages/tack-vis/src/tack/dataset/spaces.py

@on_inputs_device
def cell_keys(self):
    """Each cell's order, an i32 array (for a space with an order per cell)."""
    if not self.varies:
        return None
    if "_orders" not in self.__dict__:
        host = self.order.numpy() if hasattr(self.order, "numpy") else self.order
        host = np.asarray(host).reshape(-1)
        if host.shape[0] != self.topology.num_cells:
            raise ValueError(f"an order per cell: {self.topology.num_cells} orders, "
                             f"not {host.shape[0]}")
        if not np.isin(host, (0, 1)).all():
            bad = sorted(set(host.tolist()) - {0, 1})
            raise NotImplementedError(f"orders {bad} are not supported yet: L2 "
                                      "has orders 0 and 1")
        orders = tack.empty(host.shape, tack.i32)
        if host.size:
            tack.copyto(orders, host.astype(np.int32))
        self._orders = orders
    return self._orders

tack.dataset.Constant

Bases: Space

One value per cell (L2, order 0). Today's cell data.

tack.dataset.Values

Bases: Space

One value per entity of kind on -- points, edges, faces or cells -- with no basis: data about the entity, not a function to interpolate.

On faces, oriented=True says the value is measured along the face's stored normal (out of its side 0): a normal flux or velocity component, as MPAS's normalVelocity or upwind_flux's output. Turning a face around -- as threshold must when it keeps only a face's side-1 cell -- negates it. A plain face value (an area, an id, a quality) is left alone.

tack.dataset.SideTraces

Bases: Space

Per side: each face holds, for each of its sides, a value at each of its points (in the face's own row) -- a field's trace from each cell onto the face, MFEM's double-valued face E-vector. Laid out (face * 2 + side) * 4 + point, triangles leaving the fourth slot unused and boundary faces side 1. tack.dataset.algorithms.traces fills one from any field with a basis.

Sets, ghost marks and global ids

tack.dataset.Set

Set(data, on, ids)

Some of a topology's on entities -- "points", "cells", "faces" or "edges" -- by id: Set(data, on, ids), ids a field, an array or a list. Duplicates collapse and the ids are kept sorted (ids, a field of the topology's id_dtype); an id outside the topology is refused.

A face set is also a side set: on the boundary a face has one side, so its face names its cell's side; Set.from_sides takes (cell, local face) pairs, as Exodus side sets give them. Put a set in data.sets to name it.

Source code in packages/tack-vis/src/tack/dataset/sets.py

@on_inputs_device
def __init__(self, data, on, ids):
    topology = getattr(data, "topology", data)
    if on not in ("points", "cells", "faces", "edges"):
        raise ValueError(f"a set selects points, cells, faces or edges, not {on!r}")
    self.topology = topology
    self.on = on
    n = _count(topology, on)
    idt = topology.id_dtype
    given = id_types.as_ids(ids, id_types.at_least(idt, n))
    flags = tack.zeros((n,), tack.i32)
    if given.shape[0]:
        bad = tack.zeros((1,), tack.i32)
        _mark(given, flags, bad, n)
        if bad[0]:
            raise ValueError(f"{bad[0]} ids are outside the topology's {n} {on}")
    self._adopt(flags)

device property

device

The tack.Device its ids are on.

from_flags classmethod

from_flags(data, on, flags)

The on entities whose flag is not zero: flags holds one value per entity -- a Tack array of any integer or floating type, a field (its values), a NumPy array or a list -- as a kernel that tests each entity writes them. extract_cells of a cell set keeps those cells, so a filter that selects cells needs no compaction of its own.

Source code in packages/tack-vis/src/tack/dataset/sets.py

@classmethod
@on_inputs_device
def from_flags(cls, data, on, flags):
    """The ``on`` entities whose flag is not zero: ``flags`` holds one value per
    entity -- a Tack array of any integer or floating type, a field (its
    values), a NumPy array or a list -- as a kernel that tests each entity
    writes them. ``extract_cells`` of a cell set keeps those cells, so a
    filter that selects cells needs no compaction of its own."""
    from tack.arrays import materialize, size_of
    from tack.dataset.dataset import Field

    topology = getattr(data, "topology", data)
    if on not in ("points", "cells", "faces", "edges"):
        raise ValueError(f"a set selects points, cells, faces or edges, not {on!r}")
    n = _count(topology, on)
    if isinstance(flags, Field):
        flags = flags.values
    if isinstance(flags, (list, tuple, np.ndarray)):
        host = np.asarray(flags)
        flags = tack.array(host.astype(np.int32) if host.dtype == bool else host)
    values = materialize(flags)
    if size_of(values) != n:
        raise ValueError(f"from_flags() needs a flag for each of the topology's {n} "
                         f"{on}, not {size_of(values)}")
    out = tack.empty((n,), tack.i32)
    if n:
        _as_flags(values, out)
    return cls._from_flags(topology, on, out)

from_sides classmethod

from_sides(data, cells, local_faces)

The faces named by (cell, local face) pairs -- local faces numbered as the cell's view numbers them (VTK's order on the shape path, the cell's own face list on the polyhedral one). Which side was named is not kept: on the boundary, a face has only the one.

Source code in packages/tack-vis/src/tack/dataset/sets.py

@classmethod
@on_inputs_device
def from_sides(cls, data, cells, local_faces):
    """The faces named by (cell, local face) pairs -- local faces numbered as the
    cell's view numbers them (VTK's order on the shape path, the cell's own
    face list on the polyhedral one). Which side was named is not kept: on
    the boundary, a face has only the one."""
    from tack.dataset.dataset import for_each

    cells = np.asarray(cells, np.int64).reshape(-1)
    local = np.asarray(local_faces, np.int64).reshape(-1)
    if cells.shape != local.shape:
        raise ValueError("one local face per cell")
    n = data.num_cells
    if cells.size and (cells.min() < 0 or cells.max() >= n):
        raise ValueError(f"cells outside the topology's {n}")
    if local.size and local.min() < 0:
        raise ValueError("local faces are numbered from 0")
    topology = data.topology
    faces = topology.faces()
    flags = tack.zeros((faces.num_faces,), tack.i32)
    bad = tack.zeros((1,), tack.i32)
    if hasattr(topology, "cell_faces"):
        # Cells list their faces, as many as they like: look each pair up.
        if cells.size:
            _listed_faces(tack.array(cells), tack.array(local),
                          topology.cell_offsets, topology.cell_faces, flags, bad)
    else:
        # Shapes have at most six faces: a bit per local face, a mask per cell.
        if local.size and local.max() >= 63:
            bad = tack.full((1,), int((local >= 63).sum()), tack.i32)
        else:
            host = np.zeros(n, np.int64)
            np.bitwise_or.at(host, cells, np.left_shift(np.int64(1), local))
            masks = tack.empty((n,), tack.i64)
            if n:
                tack.copyto(masks, host)
            if cells.size:
                for_each(_side_faces, data, "cells", masks, flags, bad)
    if bad[0]:
        raise ValueError(f"{bad[0]} cells are given a local face they do not have")
    return cls._from_flags(data.topology, "faces", flags)

flags

flags()

Membership: an i32 array, 1 for each entity in the set, else 0.

Source code in packages/tack-vis/src/tack/dataset/sets.py

@on_inputs_device
def flags(self):
    """Membership: an i32 array, 1 for each entity in the set, else 0."""
    if self._flags is None:
        n = _count(self.topology, self.on)
        flags = tack.zeros((n,), tack.i32)
        if self.ids.shape[0]:
            _mark(self.ids, flags, tack.zeros((1,), tack.i32), n)
        self._flags = flags
    return self._flags

cell_keys

cell_keys()

Each cell's key, for splitting the cells' groups: its membership.

Source code in packages/tack-vis/src/tack/dataset/sets.py

def cell_keys(self):
    """Each cell's key, for splitting the cells' groups: its membership."""
    return self.flags()

numpy

numpy()

The ids, sorted, on the host.

Source code in packages/tack-vis/src/tack/dataset/sets.py

def numpy(self):
    """The ids, sorted, on the host."""
    return self.ids.numpy()

__array__

__array__(dtype=None, copy=None)

np.asarray(set): the ids, sorted, on the host.

Source code in packages/tack-vis/src/tack/dataset/sets.py

def __array__(self, dtype=None, copy=None):
    """``np.asarray(set)``: the ids, sorted, on the host."""
    if copy is False:
        raise ValueError("a Set cannot be viewed by NumPy without a copy")
    ids = self.numpy()
    return ids if dtype is None else ids.astype(dtype, copy=False)

tack.dataset.Ghosts

Ghosts(cells=None, points=None)

Ghost marks of a dataset's cells and points: one u8 per cell (cells) and per point (points), either of which may be None, holding VTK's vtkGhostType bits.

A cell or point marked DUPLICATECELL or DUPLICATEPOINT is a copy of one another piece owns. HIDDENCELL and HIDDENPOINT mark data that is not valid there (a blanked cell). REFINEDCELL marks a cell a finer AMR level covers, whose data are valid but superseded, and REFINEDPOINT (a Tack bit; VTK has no such point bit) a point at whose location a finer level has one. Filters that make cells leave hidden and refined cells out (MASKED_CELLS); reductions leave out every skipped bit. The other bits are carried, not read.

Source code in packages/tack-vis/src/tack/dataset/ghosts.py

def __init__(self, cells=None, points=None):
    self.cells = _marks(cells, "cells")
    self.points = _marks(points, "points")
    _ = self._device                    # refuses marks on two devices

device property

device

The tack.Device the marks are on.

to_device

to_device(device)

These marks on device.

Source code in packages/tack-vis/src/tack/dataset/ghosts.py

def to_device(self, device):
    """These marks on ``device``."""
    target = as_device(backend_of(device))
    return Ghosts(None if self.cells is None else self.cells.to_device(target),
                  None if self.points is None else self.points.to_device(target))

tack.dataset.GlobalIds

GlobalIds(cells=None, points=None)

Global ids of a dataset's points (points) and cells (cells): one i64 each, either of which may be None. A point or cell several pieces hold has the same id in each; ids are unique across the pieces otherwise.

Source code in packages/tack-vis/src/tack/dataset/global_ids.py

def __init__(self, cells=None, points=None):
    self.cells = _ids(cells, "cells")
    self.points = _ids(points, "points")
    _ = self._device                    # refuses ids on two devices

device property

device

The tack.Device the ids are on.

to_device

to_device(device)

These ids on device.

Source code in packages/tack-vis/src/tack/dataset/global_ids.py

def to_device(self, device):
    """These ids on ``device``."""
    target = as_device(backend_of(device))
    return GlobalIds(None if self.cells is None else self.cells.to_device(target),
                     None if self.points is None else self.points.to_device(target))

Implicit functions

tack.dataset.implicit

Implicit functions: regions of space by the sign of a function.

value(p) is negative inside, zero on the surface and positive outside, in the forms VTK (vtkImplicitFunction) and Viskores (viskores::ImplicitFunction) use, so a slice by one lands on the same points as theirs: a sphere and a cylinder are quadrics (squared distance minus squared radius), a plane and a set of planes signed distances (a set's the largest of its planes'), a box the signed distance to it.

Each is a @tack.device_class object whose parameters are instance values: a kernel takes it as a template and calls f.value(p) on a 3-vector, and moving a sphere or resizing a box does not compile anything anew. algorithms.implicit_values evaluates one on a dataset's geometry; slice, clip, extract_geometry and extract_points take one. Union, Intersection and Difference combine them, as VTK's vtkImplicitBoolean does (the smaller value, the larger, the first's or the second's negated), on the device: a kernel inlines the whole tree.

Parameters are positions in world coordinates. On a dataset with an origin (DataSet.origin) a filter evaluates the function at the geometry's values -- positions relative to that origin -- through relative_to: each function here moves its parameters by the origin in double precision (_relative_to), so nothing is lost far from zero; any other object with a value gets the origin added to the position it is given, in the geometry's type.

Every product in a value is named, so the float contract rounds it on its own and no backend fuses it into a sum: a function gives the same bits on every backend. A slice through mesh points -- whose values are then exactly the level, and merge (contour) -- cuts alike everywhere; with a fused multiply-add, Metal's noise around zero cut a plane through grid points into 1015 triangles where the CPU made 883.

tack.dataset.Plane

Plane(origin, normal)

The half-space below a plane: n . (p - origin), n the unit normal.

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, origin, normal):
    self.ox, self.oy, self.oz = map(float, _vector(origin, "origin"))
    self.nx, self.ny, self.nz = map(float, _unit(normal, "normal"))

tack.dataset.Sphere

Sphere(center, radius)

A ball: |p - center|^2 - radius^2.

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, center, radius):
    self.cx, self.cy, self.cz = map(float, _vector(center, "center"))
    if not radius > 0:
        raise ValueError(f"radius must be positive, not {radius!r}")
    self.r2 = float(radius) ** 2

tack.dataset.Cylinder

Cylinder(center, axis, radius)

An infinite cylinder about the line through center along axis: the squared distance from that line minus radius^2.

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, center, axis, radius):
    self.cx, self.cy, self.cz = map(float, _vector(center, "center"))
    self.ax, self.ay, self.az = map(float, _unit(axis, "axis"))
    if not radius > 0:
        raise ValueError(f"radius must be positive, not {radius!r}")
    self.r2 = float(radius) ** 2

tack.dataset.Box

Box(lower, upper)

An axis-aligned box from lower to upper: the signed distance to it (negative inside, by the distance to the nearest face).

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, lower, upper):
    lo, hi = _vector(lower, "lower"), _vector(upper, "upper")
    if not (lo < hi).all():
        raise ValueError(f"lower {lo} must be below upper {hi} on every axis")
    self.cx, self.cy, self.cz = map(float, (lo + hi) / 2)
    self.hx, self.hy, self.hz = map(float, (hi - lo) / 2)

tack.dataset.Planes

Planes(origins, normals)

The convex region below every one of several planes (a frustum, a clipping box at any angle): the largest of their signed distances.

The planes follow the data's precision, as the other functions' Python numbers do: implicit_values, and every filter that takes an implicit function, evaluate a copy holding them in the geometry's floating type (f64 for f64 data) on its device (_matching). Made, it holds them in f32 on the default device, which a kernel of one's own calling value reads.

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, origins, normals):
    origins = np.asarray(origins, float).reshape(-1, 3)
    normals = np.asarray(normals, float).reshape(-1, 3)
    if len(origins) == 0 or origins.shape != normals.shape:
        raise ValueError("Planes needs one origin and one normal per plane, at least one")
    lengths = np.linalg.norm(normals, axis=1)
    if not (lengths > 0).all():
        raise ValueError("a plane's normal must not be zero")
    # Kept as given (underscored: no kernel receives them), and stored per
    # precision and device as they are asked for.
    self._rows = np.hstack([origins, normals / lengths[:, None]]).reshape(-1)
    self._kept = {}
    self.planes = self._stored(tack.f32, None)

tack.dataset.Union

Union(first, second, *more)

Bases: _Boolean

Inside either function's region: the smaller of their values.

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, first, second, *more):
    self.a = _operand(first)
    self.b = type(self)(second, *more) if more else _operand(second)

tack.dataset.Intersection

Intersection(first, second, *more)

Bases: _Boolean

Inside both functions' regions: the larger of their values.

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, first, second, *more):
    self.a = _operand(first)
    self.b = type(self)(second, *more) if more else _operand(second)

tack.dataset.Difference

Difference(first, second, *more)

Bases: _Boolean

Inside the first function's region and outside the others': the first's value or the second's negated, whichever is larger.

Source code in packages/tack-vis/src/tack/dataset/implicit.py

def __init__(self, first, second, *more):
    # first - (second | more...): the rest as one union.
    self.a = _operand(first)
    self.b = Union(second, *more) if more else _operand(second)

Locating points

tack.dataset.CellLocator

CellLocator(data, density=1.0, levels=2)

Locates points in data's cells; build once, find many times. density is the bins per cell of the uniform grid a search from the cells buckets the points in; the search from the points reads the cell bins (cell_bins): with levels=2 (the default) a two-level grid, as Viskores' CellLocatorTwoLevel, which holds meshes graded toward a body in a large far field; with levels=1 that uniform grid, as Viskores' CellLocatorUniformBins and VTK's vtkStaticCellLocator, fewer candidates a search on a small, even mesh, and a few bins holding most of the cells on a graded one.

The bins follow the cells' bounding boxes when it is made. It records the geometry it was made from and that geometry's state (serves): probe and advection take it for any dataset of the same topology and geometry, and refuse it once Tack has written the positions since (a kernel, fill or copyto into them). Without one they use the locator the topology keeps (cached_locator), made again when the geometry changes. Writes Tack does not see -- through memory shared by DLPack, export_memory or wrap_ptr -- need a new locator.

Source code in packages/tack-vis/src/tack/dataset/locator.py

@on_inputs_device
def __init__(self, data, density=1.0, levels=2):
    from tack.dataset.dataset import DataSet

    if levels not in (1, 2):
        raise ValueError(f"a cell locator has 1 or 2 levels of bins, not {levels!r}")
    self.levels = levels
    #: Cells with no reference element -- polyhedra, polygons, polylines --
    #: have no parametric coordinates: points are located in a piece of each
    #: cell's fan (a polyline's segment), which ``find`` says (``tets``).
    self.polyhedral = not getattr(data.topology, "reference_cells", True)
    #: The cells and their positions, without the fields: what the bins and
    #: searches read, held by a locator a topology keeps.
    self.data = DataSet(data.topology, data.geometry, **data._frame())
    self.density = density
    self._geometry = data.geometry.values
    self._revision = arrays._revision(self._geometry)
    n = data.num_cells
    dtype = data.dtype
    self.lo = tack.empty((n,), tack.vec(3, dtype))
    self.hi = tack.empty((n,), tack.vec(3, dtype))
    if n and self.polyhedral:
        bad = tack.zeros((n,), tack.i32)
        dimension = data.topology.dimension
        if dimension >= 2:
            # A polyline needs no fan, so has no such failure.
            for_each(_not_star_shaped if dimension == 3 else _not_star_shaped_polygons,
                     data, "cells", bad)
        count = int(bad.numpy().sum())
        if count:
            raise NotImplementedError(
                f"{count} of the {n} cells are not star-shaped from their cell point "
                "(the mean of their points), or have no volume: locating points in "
                "them needs another inside test, not built yet")
        for_each(_polyhedral_boxes, data, "cells", self.lo, self.hi)
    elif n:
        for_each(_cell_boxes, data, "cells", self.lo, self.hi)
    if n:
        # On the device: reading the boxes back to take their bounds there
        # was 1.0 of the M6's 1.4 s build on an M1 Max's GPU (NumPy's
        # strided reductions, and the copy, over PCIe elsewhere).
        if self.lo.device.type == "cpu":
            count = -(-n // min(n, _BOUNDS_BLOCKS))
            k, first, step = -(-n // count), count, 1
        else:
            k = min(n, _BOUNDS_THREADS)
            first, step, count = 1, k, -(-n // k)
        lows = tack.empty((k,), tack.vec(3, dtype))
        highs = tack.empty((k,), tack.vec(3, dtype))
        _bounds(self.lo, self.hi, first, step, count, lows, highs)
        low = lows.numpy().min(axis=0).astype(float)
        high = highs.numpy().max(axis=0).astype(float)
    else:
        low = high = np.zeros(3)
    extent = high - low
    size = max(float(extent.max()), 1e-30)
    live = extent > 1e-9 * size
    # About density * n bins, cubes along the axes the cells span.
    side = (np.prod(extent[live]) / max(density * n, 1)) ** (1.0 / max(live.sum(), 1))
    dims = np.where(live, np.clip(np.ceil(extent / max(side, 1e-30)), 1, 1024), 1)
    self.dims = dims.astype(np.int64)
    self.origin = low
    self.inverse = np.where(live, self.dims / np.where(live, extent, 1.0), 0.0)
    # A point this close to the grid can be within a cell's tolerance.
    self.slack = _TOLERANCE * float(np.linalg.norm(extent))
    # The coarse grid of the two-level one (cell_bins), as the uniform grid
    # above with about _CELLS_PER_BIN cells a bin; with one level, that grid.
    per_bin = _CELLS_PER_BIN if levels == 2 else 1.0 / density
    side = (np.prod(extent[live]) / max(n / per_bin, 1)) ** (1.0 / max(live.sum(), 1))
    coarse = np.where(live, np.clip(np.ceil(extent / max(side, 1e-30)), 1, 1024), 1)
    self.coarse_dims = coarse.astype(np.int64)
    self.coarse_inverse = np.where(live, coarse / np.where(live, extent, 1.0), 0.0)
    # Its bins' sides over the sides-th root of their volume, for the leaves.
    widths = np.where(live, extent / coarse, 0.0)
    sides = max(int(live.sum()), 1)
    root = float(np.prod(widths[live])) ** (1.0 / sides) if live.any() else 1.0
    self._units = (*map(float, widths / root), float(sides))

serves

serves(data)

Whether it locates in data's cells as they are: the topology and the geometry's values it was made from, with nothing written to them by Tack since.

Source code in packages/tack-vis/src/tack/dataset/locator.py

def serves(self, data):
    """Whether it locates in ``data``'s cells as they are: the topology and the
    geometry's values it was made from, with nothing written to them by Tack
    since."""
    return (data.topology is self.data.topology
            and data.geometry.values is self._geometry
            and arrays._revision(self._geometry) == self._revision)

cell_bins

cell_bins()

The cells each leaf may hold, for queries that run from the points (a particle asking for its cell): (offsets, cells), CSR by leaf, the cells whose boxes, widened by their tolerance, overlap the leaf. With two levels the leaves are a two-level grid: a coarse uniform grid of about _CELLS_PER_BIN cells a bin, each bin a uniform grid of leaves sized by its count (leaves); with one, the uniform grid of about density bins a cell, each bin one leaf. Made on first use and kept.

Source code in packages/tack-vis/src/tack/dataset/locator.py

def cell_bins(self):
    """The cells each leaf may hold, for queries that run from the points (a
    particle asking for its cell): ``(offsets, cells)``, CSR by leaf, the
    cells whose boxes, widened by their tolerance, overlap the leaf. With
    two levels the leaves are a two-level grid: a coarse uniform grid of
    about ``_CELLS_PER_BIN`` cells a bin, each bin a uniform grid of leaves
    sized by its count (``leaves``); with one, the uniform grid of about
    ``density`` bins a cell, each bin one leaf. Made on first use and kept."""
    if "_cell_bins" not in self.__dict__:
        data = self.data
        n = data.num_cells
        grid = self._coarse_grid()
        nbins = int(np.prod(self.coarse_dims))
        # Each coarse bin's cells, then its leaves (sizes, by a scan their
        # first leaf; with one level every bin is its own leaf), then each
        # cell's leaves, entered and grouped.
        dims = tack.empty((3 * nbins,), tack.i32)
        sizes = tack.empty((nbins,), tack.i64)
        if self.levels == 2:
            counts = tack.zeros((nbins,), tack.i32)
            if n:
                _coarse_counts(self.lo, self.hi, *grid, _TOLERANCE, counts)
            _leaf_dims(counts, *self._units, _CELLS_PER_LEAF, dims, sizes)
            del counts
        else:
            dims.fill(1)
            sizes.fill(1)
        leaf_starts = tack.empty((nbins,), tack.i64)
        nleaves = exclusive_scan(sizes, leaf_starts, nbins)
        del sizes
        per_cell = tack.empty((n,), tack.i64)
        starts = tack.empty((n,), tack.i64)
        total = 0
        if n:
            _leaves_per_cell(self.lo, self.hi, *grid, dims, _TOLERANCE, per_cell)
            total = exclusive_scan(per_cell, starts, n)
        del per_cell
        index = ids.at_least(data.id_dtype, max(total, nleaves))
        keys = tack.empty((total,), ids.at_least(tack.i32, nleaves))
        cells = tack.empty((total,), data.id_dtype)
        if total:
            _leaf_entries(self.lo, self.hi, *grid, dims, leaf_starts, _TOLERANCE, starts,
                          keys, cells)
        # A leaf's cells in any order: the search ranks its candidates itself.
        order, offsets = bucket_groups(keys, nleaves, index)
        self._leaves = (dims, leaf_starts)
        self._cell_bins = (offsets, gather(cells, order) if total else cells)
    return self._cell_bins

leaves

leaves()

The two-level grid's leaves (cell_bins): (dims, starts), each coarse bin's leaves along x, y and z (dims[3 * b + axis]) and its first leaf's index.

Source code in packages/tack-vis/src/tack/dataset/locator.py

def leaves(self):
    """The two-level grid's leaves (``cell_bins``): ``(dims, starts)``, each
    coarse bin's leaves along x, y and z (``dims[3 * b + axis]``) and its
    first leaf's index."""
    self.cell_bins()
    return self._leaves

find

find(points)

Locations(cells, pcs, tets) for points (an (n, 3) host array or a field of 3-vectors): each point's cell (-1 for none) and its coordinates there, and in polyhedra its sub-tetrahedron (tets is None otherwise). See Locations.

Source code in packages/tack-vis/src/tack/dataset/locator.py

@on_inputs_device
def find(self, points):
    """``Locations(cells, pcs, tets)`` for ``points`` (an ``(n, 3)`` host array or
    a field of 3-vectors): each point's cell (-1 for none) and its
    coordinates there, and in polyhedra its sub-tetrahedron (``tets`` is
    ``None`` otherwise). See ``Locations``."""
    return self._find(points)[0]

tack.dataset.Locations

Bases: NamedTuple

Where CellLocator.find put each point: cells, the cell holding it (-1 for none), in the dataset's id_dtype; pcs, its coordinates there; and tets, in polyhedra, its sub-tetrahedron in the cell (i32, -1 for none), None on cells with a reference element. On those pcs are parametric coordinates; in a polyhedron they are the weights of the sub-tetrahedron's face point and edge points (views._PolyhedralLocation), which a field's view evaluates with tet_value(c, t, pc).

tack.dataset.PointLocator

PointLocator(
    data, bin_size=None, points_per_bin=4, periodic=None
)

Uniform bins over data's points, for finding points near points; build once, query many times. bin_size is a bin's side (in the geometry's units); by default bins hold about points_per_bin points each, cubes along the axes the points span.

periodic is a periodic box, in world coordinates: its side (a number, the box from zero), its three sides, or (lower, upper) corners. The bins then tile the box exactly, points outside it are wrapped in, and every distance is to a point's nearest image.

It records the geometry it was made from and that geometry's state (serves), as CellLocator does: filters take it for any dataset of the same topology and geometry, and refuse it once Tack has written the positions since.

Source code in packages/tack-vis/src/tack/dataset/point_locator.py

@on_inputs_device
def __init__(self, data, bin_size=None, points_per_bin=4, periodic=None):
    from tack.dataset.algorithms._common import _length
    from tack.dataset.carry import _point_values
    from tack.dataset.dataset import DataSet
    from tack.dataset.spaces import H1

    bin_size = _length(data, bin_size, "the bin size")
    periodic = _length(data, periodic, "the periodic box")
    if not isinstance(data.geometry.space, H1):
        raise TypeError("a PointLocator needs an H1 geometry: one position per point")
    if bin_size is not None and not bin_size > 0:
        raise ValueError(f"bin_size must be positive, not {bin_size!r}")
    #: The points, without the fields.
    self.data = DataSet(data.topology, data.geometry, **data._frame())
    self._geometry = data.geometry.values
    self._revision = arrays._revision(self._geometry)
    n = data.num_points
    self.positions = arrays.materialize(_point_values(data.geometry.values, n))
    #: The periodic box, ``(lower, upper)`` in world coordinates, or None.
    self.periodic = periodic_box(periodic)
    if self.periodic is not None:
        low = np.asarray(self.periodic[0]) - np.asarray(data.origin)
        high = np.asarray(self.periodic[1]) - np.asarray(data.origin)
    elif n:
        x = arrays.to_host(self.positions).reshape(-1, 3).astype(np.float64)
        low, high = x.min(axis=0), x.max(axis=0)
    else:
        low = high = np.zeros(3)
    extent = high - low
    size = max(float(extent.max()), 1e-30)
    live = extent > 1e-9 * size
    if bin_size is None:
        volume = float(np.prod(extent[live])) if live.any() else 0.0
        side = (volume * points_per_bin / max(n, 1)) ** (1.0 / max(int(live.sum()), 1))
    else:
        side = float(bin_size)
    side = max(side, 1e-30)
    # Bins that tile a periodic box exactly (at most the side asked), or
    # cover the points.
    rounding = np.floor if self.periodic is not None else np.ceil
    dims = np.where(live, np.clip(rounding(extent / side), 1, _AXIS_BINS), 1)
    while live.any() and np.prod(dims) > max(_BINS_PER_POINT * n, 1):
        side *= 2.0                  # too fine for the points: coarser bins
        dims = np.where(live, np.clip(rounding(extent / side), 1, _AXIS_BINS), 1)
    self.dims = dims.astype(np.int64)
    self.origin = low
    self.inverse = np.where(live, self.dims / np.where(live, extent, 1.0), 0.0)
    sides = np.where(live, extent / self.dims, np.inf)
    self.side = float(sides.min()) if live.any() else 1e30
    self._box = None if self.periodic is None else extent
    nbins = int(np.prod(self.dims))
    keys = tack.empty((n,), tack.i64)
    self.order = self.starts = tack.empty((1,), tack.i64)    # for the binning grid
    if n:
        _bins(self.grid(), keys)
    self.order, self.starts = bucket_order(keys, nbins)

serves

serves(data, periodic=None)

Whether it finds data's points as they are -- the topology and the geometry's values it was made from, with nothing written to them by Tack since -- in the periodic box given (or none).

Source code in packages/tack-vis/src/tack/dataset/point_locator.py

def serves(self, data, periodic=None):
    """Whether it finds ``data``'s points as they are -- the topology and the
    geometry's values it was made from, with nothing written to them by Tack
    since -- in the periodic box given (or none)."""
    return (data.topology is self.data.topology
            and data.geometry.values is self._geometry
            and arrays._revision(self._geometry) == self._revision
            and self.periodic == periodic_box(periodic))

grid

grid()

The sorted points and the grid, as a template for kernels (_PointGrid).

Source code in packages/tack-vis/src/tack/dataset/point_locator.py

def grid(self):
    """The sorted points and the grid, as a template for kernels (``_PointGrid``)."""
    return _PointGrid(self.positions, self.order, self.starts, self.origin,
                      self.inverse, self.dims, self.side, self._box)

check_radius

check_radius(radius)

Refuse a search radius of half a periodic box or more: a point's nearest image is then not the only one within it.

Source code in packages/tack-vis/src/tack/dataset/point_locator.py

def check_radius(self, radius):
    """Refuse a search radius of half a periodic box or more: a point's nearest
    image is then not the only one within it."""
    if self._box is not None and radius >= 0.5 * float(np.min(self._box)):
        raise ValueError(f"a radius of {radius} reaches half the periodic box "
                         f"({float(np.min(self._box))} across): search less far")

neighbors

neighbors(radius, points=None)

Every point within radius (inclusive) of each query point, CSR: (offsets, ids), the neighbours of query q at ids[offsets[q]: offsets[q + 1]], by bin and then by index. The queries are the points themselves (each leaving itself out) unless points are given.

Source code in packages/tack-vis/src/tack/dataset/point_locator.py

@on_inputs_device
def neighbors(self, radius, points=None):
    """Every point within ``radius`` (inclusive) of each query point, CSR:
    ``(offsets, ids)``, the neighbours of query ``q`` at ``ids[offsets[q]:
    offsets[q + 1]]``, by bin and then by index. The queries are the points
    themselves (each leaving itself out) unless ``points`` are given."""
    from tack.dataset.algorithms._common import _length

    radius = _length(self.data, radius, "the radius")
    if not radius >= 0:
        raise ValueError(f"radius must be at least zero, not {radius!r}")
    self.check_radius(radius)
    queries, own = self._queries(points)
    m = int(queries.shape[0])
    id_dtype = self.data.id_dtype
    counts = tack.empty((m,), id_dtype)
    offsets = tack.empty((m + 1,), id_dtype)
    if m == 0:
        offsets.fill(0)
        return offsets, tack.empty((0,), id_dtype)
    grid = self.grid()
    _count_within(grid, queries, float(radius), own, counts)
    total = exclusive_scan(counts, offsets, m, close=True)
    ids = tack.empty((int(total),), id_dtype)
    if total:
        _fill_within(grid, queries, float(radius), own, offsets, ids)
    return offsets, ids

nearest

nearest(k, points=None)

The k nearest points (1 <= k <= MAX_K) to each query point: (ids, squared_distances), each (m, k), nearest first, ties by index; -1 (and a squared distance of -1) where there are fewer points than k. The queries are the points themselves (each leaving itself out) unless points are given.

Source code in packages/tack-vis/src/tack/dataset/point_locator.py

@on_inputs_device
def nearest(self, k, points=None):
    """The ``k`` nearest points (``1 <= k <= MAX_K``) to each query point:
    ``(ids, squared_distances)``, each ``(m, k)``, nearest first, ties by
    index; -1 (and a squared distance of -1) where there are fewer points
    than ``k``. The queries are the points themselves (each leaving itself
    out) unless ``points`` are given."""
    if not (isinstance(k, (int, np.integer)) and 1 <= k <= MAX_K):
        raise ValueError(f"k is a whole number from 1 to {MAX_K}, not {k!r}")
    queries, own = self._queries(points)
    m = int(queries.shape[0])
    ids = tack.empty((m * k,), self.data.id_dtype)
    d2s = tack.empty((m * k,), self.data.dtype)
    if m:
        _nearest(self.grid(), queries, int(k), own, ids, d2s)
    return ids.reshape((m, k)), d2s.reshape((m, k))

Steps in time

tack.dataset.RunningStatistics

RunningStatistics(fields)

Bases: _Streamed

Per point or cell, statistics of fields (a name or names) over the steps so far: each update(data) adds the step's values and returns the step with <name>_mean, <name>_variance (the population variance, M2 / n), <name>_min and <name>_max beside its fields, computed by Welford's method on the device (in f64 where the device has it). Vector fields are taken component by component.

The steps must be on one mesh -- the same topology, or one of the same kind, sizes and arrays -- or the update is refused: remesh onto the first step's mesh first. Nothing of earlier steps is kept but the sums, so nothing is copied. The mean, minimum and maximum are in the field's unit, the variance in its square (kelvin squared for degrees Celsius). count is the steps added.

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def __init__(self, fields):
    self.fields = [fields] if isinstance(fields, str) else list(fields)
    if not self.fields:
        raise ValueError("name at least one field")
    self._states = None

tack.dataset.TimeDerivative

TimeDerivative(fields, copy=True)

Bases: _Streamed

Each of fields' rate of change since the previous step: update(data) returns the step with <name>_dt, (value - previous) / (time - previous time) per point or cell (a backward difference, in f64 where the device has it), in the field's unit over the time unit (a field in degrees Celsius gives kelvin per second). The first step has no previous one: its update returns None.

The previous step's values are kept as copies (copy=False keeps references, where the caller's buffers are not overwritten between steps). The steps must be on one mesh and in increasing time, with times in units that convert (_elapsed).

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def __init__(self, fields, copy=True):
    super().__init__(fields)
    self.copy = copy

tack.dataset.Advection

Advection(
    field,
    seeds=None,
    integrator="rk4",
    step_size=None,
    copy=True,
)

Particles carried by field (a velocity field's name) through the steps: each update(data) moves every particle from its time to the step's, with the velocity linear in time between the previous step (kept, as a copy unless copy is false) and this one, and returns the particles at their current points -- a point cloud at the step's time with seed (each particle's number), source (the seed point it came from, or -1), release (its release time), time and status (ParticleStatus bits, as advect's).

seeds (an (n, 3) array of world positions) are released at the first step's time; inject(points, time=None) adds particles at any later time (by default the last step's). The velocity is the probe of each step's own mesh, so the mesh may change between steps, and over a PartitionedDataSet each point is probed on the rank whose cells hold it: the particles stay on the rank that released them. integrator is "rk4" or "euler"; step_size (in the steps' time unit) bounds each substep, by default one substep from step to step. A particle leaving the mesh stops where it last was, marked SPATIAL_BOUNDS; one released outside has SUCCESS cleared. Times convert between the steps' unit and the velocity's (a velocity in m/s over steps in minutes).

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def __init__(self, field, seeds=None, integrator="rk4", step_size=None, copy=True):
    if integrator not in ("rk4", "euler"):
        raise ValueError(f"integrator is 'rk4' or 'euler', not {integrator!r}")
    if step_size is not None and not step_size > 0:
        raise ValueError(f"step_size must be positive, not {step_size!r}")
    self.field, self.integrator = field, integrator
    self.step_size, self.copy = step_size, copy
    self.x = np.zeros((0, 3))
    self.time = np.zeros(0)
    self.release = np.zeros(0)
    self.seed = np.zeros(0, np.int64)
    self.source = np.zeros(0, np.int64)
    self.status = np.zeros(0, np.int32)
    self._pending = []                       # (points, time or None, sources)
    self._previous = None                    # the kept previous step
    self._last = None                        # its time, in the stream's unit
    self._unit = None                        # the stream's time unit
    self._released = 0
    if seeds is not None:
        points = np.asarray(seeds, np.float64).reshape(-1, 3)
        self.inject(points, None, np.arange(len(points)))

inject

inject(points, time=None, sources=None)

Release particles at world points at time (in the steps' time unit; by default the last step's, or the first step's before any); the next update carries each from its release time. sources numbers the seed point each came from (-1 by default). A time before the last step is refused: that interval has been integrated.

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def inject(self, points, time=None, sources=None):
    """Release particles at world ``points`` at ``time`` (in the steps' time
    unit; by default the last step's, or the first step's before any); the next
    update carries each from its release time. ``sources`` numbers the seed
    point each came from (``-1`` by default). A time before the last step is
    refused: that interval has been integrated."""
    points = np.asarray(points, np.float64).reshape(-1, 3)
    if time is not None and self._last is not None and time < self._last:
        raise ValueError(f"particles released at {time}, before the last step's time "
                         f"{self._last}: that interval is integrated")
    sources = np.full(len(points), -1, np.int64) if sources is None else \
        np.asarray(sources, np.int64).reshape(-1)
    if len(sources) != len(points):
        raise ValueError("one source per point")
    self._pending.append((points, None if time is None else float(time), sources))

update

update(data)

Carry the particles to data's time and return them (see the class).

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def update(self, data):
    """Carry the particles to ``data``'s time and return them (see the class)."""
    t1 = self._stream_time(data)
    if self._last is not None and not t1 > self._last:
        raise ValueError(f"the steps' times do not increase: {self._last} then {t1}")
    self._release(data, t1)
    if self._previous is not None:
        self._advance(self._previous, data, self._last, t1)
    self._previous = _kept_dataset(data, self.field, self.copy)
    self._last = t1
    return self._particles(data)

tack.dataset.Pathlines

Pathlines(
    field,
    seeds=None,
    integrator="rk4",
    step_size=None,
    copy=True,
)

Each particle's path through the steps: Advection (its arguments) and extend_paths, each update(data) returning the current pathlines -- a polyline per particle (seed and time on the points), a particle's first point alone until it has moved. inject releases particles as Advection.inject does.

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def __init__(self, field, seeds=None, integrator="rk4", step_size=None, copy=True):
    self.advection = Advection(field, seeds, integrator, step_size, copy)
    self.paths = None

tack.dataset.Streaklines

Streaklines(
    field,
    seeds,
    releases=1,
    integrator="rk4",
    step_size=None,
    copy=True,
)

The lines of particles released at seeds (world positions) through the steps: every update releases releases particles at each seed, at times spread evenly over the interval since the previous step (one at each step's time by default), carries them all (Advection), and returns a polyline per seed (source on the cells) through the particles it released still in the mesh, newest (at the seed) first, with source, release and time on the points; a seed with one particle so far keeps its point, in no line.

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def __init__(self, field, seeds, releases=1, integrator="rk4", step_size=None,
             copy=True):
    if int(releases) < 1:
        raise ValueError("releases is at least 1")
    self.seeds = np.asarray(seeds, np.float64).reshape(-1, 3)
    self.releases = int(releases)
    self.advection = Advection(field, None, integrator, step_size, copy)

tack.dataset.extend_paths

extend_paths(paths, particles)

paths (None, or what this returned before) with particles -- an Advection's output -- appended: each particle's points in time (seed and time on the points), and a polyline through those of each particle with two or more. Each step adds the particles that moved since (a stopped one adds nothing). It stands alone, so particles written step by step can be joined into paths afterwards.

Source code in packages/tack-vis/src/tack/dataset/temporal.py

def extend_paths(paths, particles):
    """``paths`` (``None``, or what this returned before) with ``particles`` -- an
    ``Advection``'s output -- appended: each particle's points in time (``seed``
    and ``time`` on the points), and a polyline through those of each particle
    with two or more. Each step adds the particles that moved since (a
    stopped one adds nothing). It stands alone, so particles written step by
    step can be joined into paths afterwards."""
    x, rows, pieces = _rows(particles)
    seed, time = rows["seed"].astype(np.int64), rows["time"].astype(np.float64)
    if paths is not None:
        old_seed = paths.fields["seed"].values.numpy().astype(np.int64)
        old_time = paths.fields["time"].values.numpy().astype(np.float64)
        # A particle's point is new when it is later than its path's last.
        order = np.lexsort((old_time, old_seed))
        last_of = dict(zip(old_seed[order].tolist(), old_time[order].tolist()))
        new = np.array([t_ > last_of.get(s_, -np.inf)
                        for s_, t_ in zip(seed.tolist(), time.tolist())], bool)
        x = np.concatenate([paths.world_positions(), x[new]])
        seed = np.concatenate([old_seed, seed[new]])
        time = np.concatenate([old_time, time[new]])
    order = np.lexsort((time, seed))
    like = pieces[0]
    unit = like.fields["time"].unit
    return _lines(seed[order], x[order], like._frame(), like,
                  {"seed": (seed[order], None), "time": (time[order], unit)})

Kernels over meshes

tack.dataset.for_each

for_each(kernel, data, domain, *args)

Run kernel(view, *args) once per shape in domain of data.

view is the domain's view for that shape (for i in view); every Field among args becomes its view for the same shape, so u.value(i, pc) and u.at(i) read the entity i the loop is at. Other arguments pass through.

domain may also be a Space: the kernel then runs over the space's nodes (DataSet.node_views), taking arrays, not fields.

Source code in packages/tack-vis/src/tack/dataset/dataset.py

@on_inputs_device
def for_each(kernel, data, domain, *args):
    """Run ``kernel(view, *args)`` once per shape in ``domain`` of ``data``.

    ``view`` is the domain's view for that shape (``for i in view``); every
    ``Field`` among ``args`` becomes its view for the same shape, so
    ``u.value(i, pc)`` and ``u.at(i)`` read the entity ``i`` the loop is at.
    Other arguments pass through.

    ``domain`` may also be a ``Space``: the kernel then runs over the space's
    nodes (``DataSet.node_views``), taking arrays, not fields.
    """
    if isinstance(domain, Space):
        if any(isinstance(arg, Field) for arg in args):
            raise TypeError("a kernel over a space's nodes takes arrays, not fields: "
                            "pass field.values")
        for view in data.node_views(domain):
            kernel(view, *args)
        return
    kind = data._kind(domain)
    for arg in args:
        if isinstance(arg, Field):
            if arg.space.topology is not data.topology:
                raise ValueError(f"a field on {arg.space!r} of another topology")
            _check_domain(arg.space, kind)
    for group in data.launch_groups(domain, [a for a in args if isinstance(a, Field)]):
        if not group.count:
            continue
        view = data.domain_view(domain, group)
        kernel(view, *(a.view(group) if isinstance(a, Field) else a for a in args))

Filters of your own

What the built-in filters are written with. A filter that makes a new dataset states how its output's points, cells, faces and edges come from the input's, and carry takes the fields, sets, ghost marks, global ids, origin, unit and time along by those maps. A filter that builds its DataSet some other way passes its input's marks on with keep_marks (the same cells and points) or records their loss with drop_marks, so that reductions over pieces refuse what they would count twice. Filters that make cells start from visible(data), which leaves out hidden cells, and a reduction over a piece takes its owned_values. A selection made by a kernel becomes a set with Set.from_flags, and extract_cells keeps it.

tack.dataset.carry

How a filter's output takes its input's fields.

A filter that makes a new dataset says how each kind of output entity comes from the input -- the same entity (Same), entity ids[i] of the input (Take), a piece of input cell ids[i] (Pieces: a contour's triangles, a surface's faces), or a point between two input points (Interpolate) -- and carry applies one set of rules to every field:

• point data (H1, values on points) goes through the point map. An order-2 field brings its corner values, as the filters read it. A field is interpolated only if its values are floating point; an integer field meets an Interpolate map and is left behind.
• cell data (Constant, values on cells) goes through the cell map and keeps its kind.
• L2 (DG) data is copied cell by cell, keeping each cell's order, when the cell map keeps or takes whole cells (not pieces of them) and the output has reference cells to lay it out on.
• values on faces go through the face map -- negated where the map turns a face around, if they are oriented -- or become values on the output's cells when those cells are input faces (faces_to_cells).
• values on edges go through the edge map, which only a filter whose edges are the input's gives (as_polyhedra), or become values on the output's cells when those cells are input edges (edges_to_cells).
• traces, and fields on any entity without a map, are left behind: the output's own edges and faces are derived afresh, so the old numbering means nothing there.

The geometry follows the same rules unless the filter computed positions itself. fields selects which fields come along: None (all), a list of names, or an empty list (none); naming a field the input does not have raises KeyError.

carry

carry(
    data,
    out,
    *,
    points=None,
    cells=None,
    faces=None,
    faces_to_cells=None,
    edges=None,
    edges_to_cells=None,
    geometry=None,
    fields=None,
    sets=None,
)

A dataset on topology out with data's fields, by the maps given.

points, cells, faces and edges say how out's come from data's (Same, Take, Pieces for cells, Interpolate for points); faces_to_cells (edges_to_cells) that out's cells are data's faces (edges). A kind of entity with no map takes no fields. geometry is the output's geometry field when the filter computed it, otherwise the input's is carried like any point (or L2) field.

Ghost marks go by the same maps (ghosts.carry_marks): what the output makes takes the owner of what it came from. Global ids go with the entities that are input entities (global_ids.carry_ids).

Sets go by the same maps: an output entity is in a set when the input entity it comes from is, and a face or edge set becomes a cell set where faces or edges become cells (faces_to_cells, edges_to_cells). A set whose kind has no map, or one that makes new entities (Interpolate), is left behind; sets, if given, replaces them.

Source code in packages/tack-vis/src/tack/dataset/carry.py

@on_inputs_device
def carry(data, out, *, points=None, cells=None, faces=None, faces_to_cells=None,
          edges=None, edges_to_cells=None, geometry=None, fields=None, sets=None):
    """A dataset on topology ``out`` with ``data``'s fields, by the maps given.

    ``points``, ``cells``, ``faces`` and ``edges`` say how ``out``'s come from
    ``data``'s (``Same``, ``Take``, ``Pieces`` for cells, ``Interpolate`` for
    points); ``faces_to_cells`` (``edges_to_cells``) that ``out``'s cells are
    ``data``'s faces (edges). A
    kind of entity with no map takes no fields. ``geometry`` is the output's
    geometry field when the filter computed it, otherwise the input's is
    carried like any point (or L2) field.

    Ghost marks go by the same maps (``ghosts.carry_marks``): what the output makes
    takes the owner of what it came from. Global ids go with the entities that are
    input entities (``global_ids.carry_ids``).

    Sets go by the same maps: an output entity is in a set when the input entity
    it comes from is, and a face or edge set becomes a cell set where faces or
    edges become cells (``faces_to_cells``, ``edges_to_cells``). A set whose kind has no map, or one that makes new
    entities (``Interpolate``), is left behind; ``sets``, if given, replaces
    them.
    """
    if geometry is None:
        geometry = _carry_field(data.geometry, data, out, points, cells, faces,
                                faces_to_cells)
        if geometry is None:
            raise ValueError("the geometry cannot be carried by these maps; pass it")
    carried = {}
    for name in selected(data, fields):
        field = _carry_field(data.fields[name], data, out, points, cells, faces,
                             faces_to_cells, edges, edges_to_cells)
        if field is not None:
            carried[name] = field
    if sets is None:
        sets = {}
        for name, member in data.sets.items():
            moved = _carry_set(member, out, points, cells, faces, faces_to_cells, edges,
                               edges_to_cells)
            if moved is not None:
                sets[name] = moved
    from tack.dataset.ghosts import carry_marks
    from tack.dataset.global_ids import carry_ids

    # A geometry the filter computed is in the input's spatial unit.
    made = DataSet(out, geometry, fields=carried, sets=sets, **data._frame())
    carry_ids(data, made, points=points, cells=None if faces_to_cells is not None or
              edges_to_cells is not None else cells)
    return carry_marks(data, made, points=points, cells=cells, faces_to_cells=faces_to_cells,
                       edges_to_cells=edges_to_cells)

tack.dataset.Same

Output entity i is input entity i.

tack.dataset.Take

Take(ids, turned=None)

Output entity i is input entity ids[i]. turned, for faces: an i32 array, 1 where the output face is the input face turned around, so its oriented values change sign.

Source code in packages/tack-vis/src/tack/dataset/carry.py

def __init__(self, ids, turned=None):
    self.ids = _id_field(ids)
    self.turned = turned if turned is None or isinstance(turned, _DeviceArray) \
        else tack.array(np.asarray(turned, np.int32))

tack.dataset.Pieces

Pieces(ids)

Bases: Take

Output cell i lies inside input cell ids[i]: it takes that cell's values, but not its L2 data, which belongs to the whole cell's shape.

Source code in packages/tack-vis/src/tack/dataset/carry.py

def __init__(self, ids):
    super().__init__(ids)

tack.dataset.Interpolate

Interpolate(
    ends, weights, averages=None, average_cells=None
)

Output point i lies between input points ends[i] (a 2-vector of ids) at weights[i] from the first. averages, (offsets, members): further output points after those, each the average of earlier output points members[offsets[k]:offsets[k + 1]] (clip's centroid points), and inside input cell average_cells[k], whose ghost mark it takes.

Source code in packages/tack-vis/src/tack/dataset/carry.py

def __init__(self, ends, weights, averages=None, average_cells=None):
    self.ends = ends
    self.weights = weights
    self.averages = averages
    self.average_cells = average_cells

tack.dataset.keep_marks

keep_marks(source, out)

out, a dataset on source's cells and points (the same topology), with source's ghost marks -- and its record of having lost them -- and its origin, spatial unit and time.

Source code in packages/tack-vis/src/tack/dataset/ghosts.py

def keep_marks(source, out):
    """``out``, a dataset on ``source``'s cells and points (the same topology), with
    ``source``'s ghost marks -- and its record of having lost them -- and its
    origin, spatial unit and time."""
    out.ghosts = source.ghosts
    out._ghosts_lost = source._ghosts_lost
    out.global_ids = source.global_ids
    out.origin = source.origin
    return _keep_frame(source, out)

tack.dataset.drop_marks

drop_marks(source, out, reason=True)

out, made from source's cells and points by a filter that cannot give it ghost marks: if source had marks, or had lost them, record that out lost them (reason, a sentence saying why, or True), so reductions refuse it rather than count shared cells and points once per piece.

Source code in packages/tack-vis/src/tack/dataset/ghosts.py

def drop_marks(source, out, reason=True):
    """``out``, made from ``source``'s cells and points by a filter that cannot give
    it ghost marks: if ``source`` had marks, or had lost them, record that
    ``out`` lost them (``reason``, a sentence saying why, or ``True``), so
    reductions refuse it rather than count shared cells and points once per
    piece."""
    if source.ghosts is not None or source._ghosts_lost:
        out._ghosts_lost = source._ghosts_lost or reason
    out.origin = source.origin
    return _keep_frame(source, out)

tack.dataset.visible

visible(data, *arguments)

data without its hidden cells (HIDDENCELL: data not valid there) and refined cells (REFINEDCELL: an AMR cell a finer level covers), which filters that make cells never process: data itself when it has none, else its other cells as extract_cells keeps them, marks and fields carried. Given field arguments too (Fields of data or names; anything else passes through), returns (data, *arguments), each Field carried onto the cells kept.

Source code in packages/tack-vis/src/tack/dataset/ghosts.py

@on_inputs_device
def visible(data, *arguments):
    """``data`` without its hidden cells (``HIDDENCELL``: data not valid there) and
    refined cells (``REFINEDCELL``: an AMR cell a finer level covers), which
    filters that make cells never process: ``data`` itself when it has none,
    else its other cells as ``extract_cells`` keeps them, marks and fields
    carried. Given field ``arguments`` too (``Field``s of
    ``data`` or names; anything else passes through), returns ``(data,
    *arguments)``, each ``Field`` carried onto the cells kept."""
    from tack.dataset.dataset import DataSet, Field

    ghosts = data.ghosts
    if ghosts is None or not _has(ghosts.cells, Ghosts.MASKED_CELLS):
        return (data, *arguments) if arguments else data
    from tack.dataset.algorithms._extract import _keep_cells

    names = {f"__tack_argument{k}__": a for k, a in enumerate(arguments)
             if isinstance(a, Field) and a.space.topology is data.topology}
    source = data
    if names:
        source = DataSet(data.topology, data.geometry, fields={**data.fields, **names},
                         sets=dict(data.sets), ghosts=data.ghosts, **data._frame())
    keep = tack.empty((data.num_cells,), tack.i32)
    keep.fill(1)
    out = _keep_cells(source, keep)
    if not arguments:
        return out
    carried = {id(a): out.fields.pop(name) for name, a in names.items() if name in out.fields}
    return (out, *(carried.get(id(a), a) for a in arguments))

tack.dataset.owned_values

owned_values(data, field)

The values of field (a Field of data, or its name) that data owns: all of them when it has no marks for where they live, else those of cells and points not marked a duplicate or hidden, materialized. Fields on faces or edges, and order-2 fields, have no marks of their own yet; on a dataset with ghost marks they are refused rather than counted twice.

Source code in packages/tack-vis/src/tack/dataset/ghosts.py

@on_inputs_device
def owned_values(data, field):
    """The values of ``field`` (a ``Field`` of ``data``, or its name) that ``data`` owns: all of
    them when it has no marks for where they live, else those of cells and
    points not marked a duplicate or hidden, materialized. Fields on faces or
    edges, and order-2 fields, have no marks of their own yet; on a dataset
    with ghost marks they are refused rather than counted twice."""
    from tack.dataset.carry import _take
    from tack.dataset.spaces import H1, L2, Constant, Values

    if isinstance(field, str):
        field = data.fields[field]
    if data._ghosts_lost:
        why = data._ghosts_lost if isinstance(data._ghosts_lost, str) else (
            "the filter that made it does not carry the marks")
        raise NotImplementedError(
            f"this dataset is a filter's output from data with ghost marks, and {why}: "
            "its shared cells and points would count once per piece. Reduce the input "
            "instead, or set the output's ghosts")
    values = materialize(field.values)
    ghosts = data.ghosts
    if ghosts is None or (ghosts.cells is None and ghosts.points is None):
        return values
    space = field.space
    n = size_of(values)
    on_points = (isinstance(space, H1) and space.order == 1) or (
        isinstance(space, Values) and space.on == "points")
    on_cells = isinstance(space, (Constant, L2)) or (
        isinstance(space, Values) and space.on == "cells")
    if on_points:
        marks, skipped = ghosts.points, Ghosts.SKIPPED_POINTS
    elif on_cells:
        marks, skipped = ghosts.cells, Ghosts.SKIPPED_CELLS
    else:
        raise NotImplementedError(f"which values of a field on {space!r} a piece owns is "
                                  "not decided yet; this dataset has ghost marks")
    if marks is None:
        return values
    flags = tack.empty((n,), tack.i32)
    if isinstance(space, L2):
        if marks.shape[0]:
            _keep_blocks(marks, skipped, space.offsets, flags)
    elif n:
        _keep(marks, skipped, flags)
    idt = data.id_dtype
    slots = tack.empty((n,), idt)
    count = exclusive_scan(flags, slots, n) if n else 0
    if count == n:
        return values
    ids = tack.empty((count,), idt)
    if count:
        _kept(flags, slots, ids)
    return _take(values, ids)

Units

A filter's parameters take quantities ("300 K", (1, "cm")) where the built-ins do: units.magnitude(value, field.unit, "the bounds") gives a number in a field's unit, units.magnitude(value, data.unit, "the radius") one in the spatial unit, and a plain number passes through. A field a filter computes gets the unit of what it is: units.delta for a difference (a difference of degrees Celsius is in kelvin), units.derived for a product of powers, units.named for the SI name. Unknown in, unknown out: None, never a guess.

tack.dataset.units

Physical units of fields: plain strings, checked and combined by unyt.

A field's unit (Field.unit) is meaning, not storage: arrays and kernels never see one. None is unknown and "dimensionless" ("1" on the way in) is a pure number; nothing derived from an unknown unit is guessed. A dataset's spatial unit (DataSet.unit) is its geometry field's, and its origin is in it.

Units are stored as unyt's canonical spelling, so "m/s", "m*s**-1" and "m/s^1" are one unit: ^ is read as a power, as LLNL units and VTK write it. Code units are scaled units ("3.0857e21*cm"), which need no registry of their own. Tack's registry is unyt's plus the US survey foot, ftUS (1200/3937 m), which State Plane coordinates and older LiDAR use and which differs from the international foot by 2 parts per million.

Derived units (derived) are products of powers of units; None in gives None out. Units with an offset zero (degrees Celsius and Fahrenheit) are not multiplied: a difference of them, as a gradient, a jump or a variance takes, is in kelvin (Rankine), the same size of degree with zero at zero (delta). named shows a derived unit by its SI name (Pa*m**3 as J) where a filter names what it computes.

Parameters take quantities (quantity, magnitude): a string such as "10 cm", a pair (value, unit) (the value a number or an array), or a unyt or pint quantity. A plain number is in the unit the parameter is in (the geometry's for a length, the field's for an isovalue), as before.

unyt (with sympy) is imported on first use, so a program that gives no unit pays nothing. Tack hands out strings only, never unyt objects.

canonical

canonical(unit)

unit (a string, or None for unknown) in unyt's canonical spelling; a string that is no unit raises ValueError.

Source code in packages/tack-vis/src/tack/dataset/units.py

def canonical(unit):
    """``unit`` (a string, or ``None`` for unknown) in unyt's canonical spelling;
    a string that is no unit raises ``ValueError``."""
    if unit is None:
        return None
    if not isinstance(unit, str):
        raise TypeError(f"a unit is a string such as 'm/s', or None, not {unit!r}")
    if not unit.strip():
        raise ValueError("an empty unit: give 'dimensionless' for a pure number, or None "
                         "for unknown")
    return str(_parsed(unit))

conversion

conversion(source, target)

(factor, offset): a value in source is factor * value + offset in target. Both must be known and of the same dimensions.

Source code in packages/tack-vis/src/tack/dataset/units.py

def conversion(source, target):
    """``(factor, offset)``: a value in ``source`` is ``factor * value + offset``
    in ``target``. Both must be known and of the same dimensions."""
    if source is None:
        raise ValueError("the unit is unknown: give the field a unit before converting it")
    if target is None:
        raise ValueError("convert to which unit? None is unknown")
    a, b = _parsed(source), _parsed(target)
    if not a.same_dimensions_as(b):
        raise ValueError(f"cannot convert {source!r} to {target!r}: {a.dimensions} and "
                         f"{b.dimensions} differ")
    factor, offset = a.get_conversion_factor(b)
    # unyt's offset is subtracted after scaling.
    return float(factor), 0.0 if offset is None else -float(offset)

same

same(a, b)

Whether units a and b (canonical strings or None) are one unit: equal, or both known and one a relabelling of the other (factor 1, no offset). Unknown is the same only as unknown.

Source code in packages/tack-vis/src/tack/dataset/units.py

def same(a, b):
    """Whether units ``a`` and ``b`` (canonical strings or ``None``) are one unit:
    equal, or both known and one a relabelling of the other (factor 1, no
    offset). Unknown is the same only as unknown."""
    if a == b:
        return True
    if a is None or b is None:
        return False
    a, b = _parsed(a), _parsed(b)
    if not a.same_dimensions_as(b):
        return False
    factor, offset = a.get_conversion_factor(b)
    return factor == 1.0 and offset is None

compatible

compatible(a, b)

Whether units a and b may meet: one unit, or either unknown (which cannot be checked, so it is not refused).

Source code in packages/tack-vis/src/tack/dataset/units.py

def compatible(a, b):
    """Whether units ``a`` and ``b`` may meet: one unit, or either unknown (which
    cannot be checked, so it is not refused)."""
    return a is None or b is None or same(a, b)

common

common(what, units)

The one unit of units (canonical strings or None), or ValueError naming what ("the pieces' geometry", say) when two known units differ. None if any is unknown.

Source code in packages/tack-vis/src/tack/dataset/units.py

def common(what, units):
    """The one unit of ``units`` (canonical strings or ``None``), or ``ValueError``
    naming ``what`` (``"the pieces' geometry"``, say) when two known units
    differ. ``None`` if any is unknown."""
    units = list(units)
    known = [unit for unit in units if unit is not None]
    for unit in known[1:]:
        if not same(known[0], unit):
            raise ValueError(f"{what} is in {known[0]!r} and {unit!r}: convert one first "
                             "(tack.dataset.convert)")
    return known[0] if known and len(known) == len(units) else None

check

check(what, a, b)

ValueError naming what unless units a and b may meet (compatible).

Source code in packages/tack-vis/src/tack/dataset/units.py

def check(what, a, b):
    """``ValueError`` naming ``what`` unless units ``a`` and ``b`` may meet
    (``compatible``)."""
    if not compatible(a, b):
        raise ValueError(f"{what}: {a!r} and {b!r} differ; convert one first "
                         "(tack.dataset.convert)")

delta

delta(unit)

The unit of a difference of values in unit: itself, or for a unit with an offset zero, the absolute unit of the same size (kelvin for degrees Celsius, Rankine for Fahrenheit). None stays None.

Source code in packages/tack-vis/src/tack/dataset/units.py

def delta(unit):
    """The unit of a difference of values in ``unit``: itself, or for a unit with an
    offset zero, the absolute unit of the same size (kelvin for degrees Celsius,
    Rankine for Fahrenheit). ``None`` stays ``None``."""
    if unit is None:
        return None
    parsed = _parsed(unit)
    if not parsed.base_offset:
        return unit
    for absolute in ("K", "R"):
        if math.isclose(parsed.base_value, _parsed(absolute).base_value, rel_tol=1e-12):
            return absolute
    return canonical(f"{parsed.base_value!r}*K")

derived

derived(*terms)

The unit of a product of powers: terms are (unit, exponent) pairs (a bare unit is to the first power). None if any unit is unknown; a unit with an offset zero is refused (take its delta first, or convert).

Source code in packages/tack-vis/src/tack/dataset/units.py

def derived(*terms):
    """The unit of a product of powers: ``terms`` are ``(unit, exponent)`` pairs (a
    bare unit is to the first power). ``None`` if any unit is unknown; a unit
    with an offset zero is refused (take its ``delta`` first, or convert)."""
    terms = [(t, 1) if not isinstance(t, tuple) else t for t in terms]
    if any(unit is None for unit, _ in terms):
        return None
    result = None
    for unit, exponent in terms:
        parsed = _parsed(unit)
        if parsed.base_offset:
            raise ValueError(f"{unit!r} has an offset zero, so its values do not multiply: "
                             "convert to an absolute unit (K) first")
        if parsed.is_dimensionless and parsed.base_value == 1.0:
            continue
        power = parsed if exponent == 1 else parsed ** exponent
        result = power if result is None else result * power
    if result is None or (result.is_dimensionless and result.base_value == 1.0):
        return "dimensionless"
    return canonical(str(result))

named

named(unit)

unit shown by the SI derived unit of its dimensions and size, if one is ("Pa*m**3" as "J", "kg*m/s**2" as "N"), else unchanged. Only the dimensions and the size decide, never how the unit was written.

Source code in packages/tack-vis/src/tack/dataset/units.py

def named(unit):
    """``unit`` shown by the SI derived unit of its dimensions and size, if one is
    (``"Pa*m**3"`` as ``"J"``, ``"kg*m/s**2"`` as ``"N"``), else unchanged. Only
    the dimensions and the size decide, never how the unit was written."""
    if unit is None:
        return None
    parsed = _parsed(unit)
    for name in _NAMED:
        candidate = _parsed(name)
        if parsed.dimensions == candidate.dimensions and \
                math.isclose(parsed.base_value, candidate.base_value, rel_tol=1e-12):
            return name
    return unit

dimensionless_or_none

dimensionless_or_none(what, unit)

"dimensionless" for a pure number, None for unknown, and ValueError naming what for anything with dimensions or a scale (a logarithm of metres, or of kilometres per metre, depends on the unit).

Source code in packages/tack-vis/src/tack/dataset/units.py

def dimensionless_or_none(what, unit):
    """``"dimensionless"`` for a pure number, ``None`` for unknown, and
    ``ValueError`` naming ``what`` for anything with dimensions or a scale (a
    logarithm of metres, or of kilometres per metre, depends on the unit)."""
    if unit is None:
        return None
    parsed = _parsed(unit)
    if parsed.is_dimensionless and parsed.base_value == 1.0:
        return "dimensionless"
    raise ValueError(f"{what} takes a pure number, not values in {unit!r}: convert to "
                     "'dimensionless' (divide by a reference) first")

quantity

quantity(value)

(magnitude, unit) of a quantity -- a string "10 cm", a pair (value, unit) with unit a string, or a unyt or pint quantity (.units and .value or .magnitude) -- or None for anything else (a plain number, an array, a field's name is not parsed here). A string with a number and no unit is a pure number.

Source code in packages/tack-vis/src/tack/dataset/units.py

def quantity(value):
    """``(magnitude, unit)`` of a quantity -- a string ``"10 cm"``, a pair ``(value,
    unit)`` with ``unit`` a string, or a unyt or pint quantity (``.units`` and
    ``.value`` or ``.magnitude``) -- or ``None`` for anything else (a plain
    number, an array, a field's name is not parsed here). A string with a number
    and no unit is a pure number."""
    if isinstance(value, str):
        match = _NUMBER.match(value)
        if match is None:
            return None
        number, unit = match.groups()
        return float(number), canonical(unit or "dimensionless")
    if isinstance(value, (tuple, list)) and len(value) == 2 and isinstance(value[1], str) \
            and not isinstance(value[0], str):
        return value[0], canonical(value[1])
    units_of = getattr(value, "units", None)
    if units_of is not None:
        number = getattr(value, "value", None)
        if number is None:
            number = getattr(value, "magnitude", None)
        if number is not None:
            if np.ndim(number) == 0:
                number = float(number)
            return number, canonical(str(units_of))
    return None

magnitude

magnitude(value, unit, what)

value as a number (or numbers) in unit: a quantity (quantity) converted, scale and offset; a sequence of quantities each; anything else as given, taken to be in unit already. Converting a quantity needs unit known; what names the parameter in the refusal.

Source code in packages/tack-vis/src/tack/dataset/units.py

def magnitude(value, unit, what):
    """``value`` as a number (or numbers) in ``unit``: a quantity (``quantity``)
    converted, scale and offset; a sequence of quantities each; anything else
    as given, taken to be in ``unit`` already. Converting a quantity needs
    ``unit`` known; ``what`` names the parameter in the refusal."""
    found = quantity(value)
    if found is None:
        if isinstance(value, (tuple, list)) and any(quantity(v) is not None
                                                    for v in value):
            return type(value)(magnitude(v, unit, what) for v in value)
        return value
    number, given = found
    if unit is None:
        raise ValueError(f"{what} is given in {given!r}, but what it is compared with has "
                         "no unit: give the dataset (or field) a unit, or pass a plain "
                         "number")
    factor, offset = conversion(given, unit)
    if factor == 1.0 and offset == 0.0:
        return number
    if isinstance(number, (int, float)):
        return factor * number + offset
    return np.asarray(number, np.float64) * factor + offset

number_and_unit

number_and_unit(value)

(number, unit) of a parameter that may bring its own unit (a mass, a velocity): a quantity's magnitude and unit, else the value as given and None.

Source code in packages/tack-vis/src/tack/dataset/units.py

def number_and_unit(value):
    """``(number, unit)`` of a parameter that may bring its own unit (a mass, a
    velocity): a quantity's magnitude and unit, else the value as given and
    ``None``."""
    found = quantity(value)
    return (value, None) if found is None else found

time_unit

time_unit(unit)

unit (a string or None) canonical, refused unless it is a unit of time.

Source code in packages/tack-vis/src/tack/dataset/units.py

def time_unit(unit):
    """``unit`` (a string or ``None``) canonical, refused unless it is a unit of
    time."""
    unit = canonical(unit)
    if unit is not None:
        unyt, _ = _unyt()
        if _parsed(unit).dimensions != unyt.dimensions.time:
            raise ValueError(f"{unit!r} is not a unit of time")
    return unit

Interoperability

tack.interop.vtk

tack.interop.vtk — Zero-copy interop between VTK arrays and Tack arrays.

from tack.interop.vtk import vtk_to_array, array_to_vtk

field = vtk_to_array(vtk_array)   # VTK  -> Tack
array = array_to_vtk(field)       # Tack -> VTK

Both directions share memory rather than copying, on the host and on the device, and each side keeps the other alive for as long as it is needed.

Shape

A vtkDataArray is tuples x components, and a Tack array is n-dimensional, so the two line up directly: a 2-D array of shape (n, 3) is n tuples of 3 components, and field.shape[1] is the component count. Nothing needs to be declared -- the shape already says it.

(1000, 3) field  <->  1000 tuples x 3 components
(1000,)   field  <->  1000 tuples x 1 component

The one wrinkle is that Tack's own visualization algorithms work on flat interleaved arrays -- flying_edges returns points as (n*3,), and compute_normals indexes points[i * 3 + c]. For those, ask for the flat form explicitly:

points = vtk_to_array(vtk_points, flatten=True)   # (n*3,)
array  = array_to_vtk(points, n_components=3)     # back to n x 3

Requirements

VTK with vtkmodules.util.dlpack_support. Device arrays additionally need VTK built with Viskores; host arrays work with any VTK that has the DLPack module.

Level Zero

CUDA and HIP pointers identify themselves: the runtime keeps one context per device for the whole process. A Level Zero pointer does not -- it means something only inside the context it was allocated from, and DLPack has no field for a context. So the two libraries have to be in the same one before any array is exchanged:

from tack.interop.vtk import init_level_zero
init_level_zero()        # instead of tack.set_default_device("level_zero")

That starts Tack's Level Zero backend inside the context VTK's Viskores device already uses. It needs VTK built with Viskores on Kokkos' SYCL backend. Fields made before the call belong to another context and cannot be shared.

init_level_zero

init_level_zero()

Make a Level Zero device in the context VTK's device memory uses the default device, and return it.

Raises RuntimeError if this VTK cannot report one: it predates dlpack_support.level_zero_handles, was built without Viskores on Kokkos' SYCL backend, or is not running on a Level Zero device.

Source code in packages/tack-vis/src/tack/interop/vtk.py

def init_level_zero():
    """Make a Level Zero device in the context VTK's device memory uses the
    default device, and return it.

    Raises RuntimeError if this VTK cannot report one: it predates
    ``dlpack_support.level_zero_handles``, was built without Viskores on
    Kokkos' SYCL backend, or is not running on a Level Zero device.
    """
    dlpack_support = _dlpack_support()
    handles = getattr(dlpack_support, "level_zero_handles", lambda: None)()
    if not handles:
        raise RuntimeError(
            "this VTK has no Level Zero context to share. Device interop on "
            "Level Zero needs VTK built with Viskores on Kokkos' SYCL "
            "backend, running on a Level Zero device.")
    return tack.set_default_device("level_zero", external_context=handles)

vtk_to_array

vtk_to_array(vtk_array, flatten=False)

Wrap a vtkDataArray as a Tack array, without copying.

Parameters:

Name	Type	Description	Default
vtk_array		any vtkDataArray, on the host or on a device.	required
flatten		return a 1-D array of tuples * components values instead of a 2-D one. Tack's visualization algorithms take their point and vector inputs in this interleaved form.	False

Returns:

Type	Description
	A Tack array sharing the VTK array's memory, shaped (tuples, components) -- or (tuples,) when there is one component, since a scalar array is naturally 1-D. The VTK array is held for as long as the array lives, so it may be dropped immediately.

The array is on the default device when that can address the memory, else on the memory's own device: a CUDA array with the CPU as default is a CUDA array.

Raises:

Type	Description
RuntimeError	for Level Zero memory from another context than Tack's (see init_level_zero).
ValueError	for an array with no memory to share (implicit and computed arrays have none) or a non-contiguous layout.

Source code in packages/tack-vis/src/tack/interop/vtk.py

def vtk_to_array(vtk_array, flatten=False):
    """Wrap a vtkDataArray as a Tack array, without copying.

    Args:
        vtk_array: any vtkDataArray, on the host or on a device.
        flatten: return a 1-D array of ``tuples * components`` values
            instead of a 2-D one. Tack's visualization algorithms take
            their point and vector inputs in this interleaved form.

    Returns:
        A Tack array sharing the VTK array's memory, shaped
        ``(tuples, components)`` -- or ``(tuples,)`` when there is one
        component, since a scalar array is naturally 1-D. The VTK array is
        held for as long as the array lives, so it may be dropped
        immediately.

    The array is on the default device when that can address the memory,
    else on the memory's own device: a CUDA array with the CPU as default is
    a CUDA array.

    Raises:
        RuntimeError: for Level Zero memory from another context than
            Tack's (see `init_level_zero`).
        ValueError: for an array with no memory to share (implicit and
            computed arrays have none) or a non-contiguous layout.
    """
    dlpack_support = _dlpack_support()

    # VTK always exports 2-D, (tuples, components). The array lands on the
    # default device if that can wrap the memory, else on the memory's own
    # device, whose context is the one to check.
    array = tack.from_dlpack(dlpack_support.vtk_to_dlpack(vtk_array))
    _require_shared_level_zero_context(dlpack_support, array._backend)

    if flatten:
        return array.reshape((array.size,))
    if len(array.shape) == 2 and array.shape[1] == 1:
        return array.reshape((array.shape[0],))
    return array

array_to_vtk

array_to_vtk(array, n_components=None, name=None)

Wrap a Tack array as a vtkDataArray, without copying.

Parameters:

Name	Type	Description	Default
array		the Tack array to share. A 2-D array maps straight over -- (1000, 3) becomes 1000 tuples of 3.	required
n_components		only needed for a flat array holding interleaved data, which is the form Tack's own algorithms produce. A (3000,) array with n_components=3 becomes 1000 tuples of 3. Leave it unset for a 2-D array, whose shape already says.	None
name		array name, as VTK uses to identify it in a dataset.	None

Returns:

Type	Description
	A vtkDataArray sharing the array's memory. The array is held for as long as the array lives.

Raises:

Type	Description
ValueError	if n_components contradicts the array's shape, does not divide its size, or the array has more than 2 dimensions.

Source code in packages/tack-vis/src/tack/interop/vtk.py

def array_to_vtk(array, n_components=None, name=None):
    """Wrap a Tack array as a vtkDataArray, without copying.

    Args:
        array: the Tack array to share. A 2-D array maps straight over --
            ``(1000, 3)`` becomes 1000 tuples of 3.
        n_components: only needed for a *flat* array holding interleaved
            data, which is the form Tack's own algorithms produce. A
            ``(3000,)`` array with n_components=3 becomes 1000 tuples of 3.
            Leave it unset for a 2-D array, whose shape already says.
        name: array name, as VTK uses to identify it in a dataset.

    Returns:
        A vtkDataArray sharing the array's memory. The array is held for as
        long as the array lives.

    Raises:
        ValueError: if n_components contradicts the array's shape, does not
            divide its size, or the array has more than 2 dimensions.
    """
    # Work out the layout before reaching for VTK, so a shape mistake
    # reports itself rather than being masked by a missing-VTK error.
    if getattr(array, 'layout', 'aos') == 'soa':
        raise TypeError(
            "array_to_vtk() shares one buffer with VTK, and an SOA array has one per "
            "component; share each of its components, or copy it to one buffer with "
            "tack.array(a, layout='aos')")
    if hasattr(array, '_vector_n'):
        width = array._vector_n
        if n_components is not None and n_components != width:
            raise ValueError(
                f"n_components={n_components} contradicts the element width {width}")
        if len(array.shape) > 1:
            raise ValueError("reshape a typed array to one array dimension for VTK tuples")
        array = array._scalar_view().reshape((array.size, width))
        n_components = width
    ndim = len(array.shape)

    if n_components is None:
        if ndim == 2:
            n_components = array.shape[1]
        elif ndim <= 1:
            n_components = 1
        else:
            raise ValueError(
                f"a {ndim}-dimensional array has no obvious VTK layout; "
                f"reshape it to (tuples, components), or pass n_components "
                f"if it holds flat interleaved data")
    else:
        if n_components < 1:
            raise ValueError(f"n_components must be at least 1, got {n_components}")
        if ndim == 2 and array.shape[1] != n_components:
            raise ValueError(
                f"n_components={n_components} contradicts the array's shape "
                f"{array.shape}, which already says {array.shape[1]}")
        if array.size % n_components:
            raise ValueError(
                f"an array of {array.size} values does not divide into tuples "
                f"of {n_components}")

    dlpack_support = _dlpack_support()
    _require_shared_level_zero_context(dlpack_support, array._backend)

    # VTK reads the tensor as (tuples, components), so give it that shape.
    shaped = array
    if array.shape != (array.size // n_components, n_components):
        shaped = array.reshape((array.size // n_components, n_components))

    return dlpack_support.dlpack_to_vtk(shaped, name=name)

vtk_to_dataset

vtk_to_dataset(
    grid,
    dtype=tack.f32,
    polyhedral=None,
    orient=False,
    id_dtype=None,
    origin=None,
)

Copy a vtkUnstructuredGrid, vtkPolyData, vtkStructuredGrid, vtkRectilinearGrid or vtkImageData into a tack.dataset.DataSet: an image's points as UniformCoordinates (with its direction matrix), a rectilinear grid's as a CartesianProduct, a structured grid's as its explicit points.

Points and floating-point arrays become fields of dtype; integer arrays keep their type. Named point data become H1 fields and named cell data Constant fields (a cell array named as a point array is "<name> (cells)", with a warning); a vtkGhostType array becomes the ghost marks (data.ghosts), not a field; a constant or affine array (VTK's implicit arrays) stays implicit, as a ConstantArray or CountingArray. An unstructured grid's cells must all be linear shapes -- unless it holds polyhedra (or polyhedral is true), when it becomes a PolyhedralTopology, every cell a polyhedron of its faces, polylines, when it becomes a PolylineTopology (its lines two-point polylines, and no other cells allowed), or polygons (VTK_POLYGON), when it becomes a PolygonalTopology (its triangles and quads polygons too, and no other cells allowed). Its types, offsets and connectivity stay implicit where VTK's are (a cell type given as a number, affine offsets). A vtkPolyData is taken as the unstructured grid of its cells in VTK's order -- verts, lines, polygons (of three or four points triangles and quads), and triangle strips, each strip as its triangles, which take its cell data -- so the same rules hold: a surface of triangles and quads is an UnstructuredTopology, one with larger polygons a PolygonalTopology, lines a PolylineTopology, and polylines or polygons beside other cells are refused. With orient, inconsistently wound polyhedra are repaired first (PolyhedralTopology.from_cell_faces); without it they are refused. VTK's ids are 64-bit; the topology keeps them in id_dtype, by default i32 unless the grid is too large for it (tack.dataset.ids).

origin (three numbers) becomes the dataset's origin: VTK's points, an image's origin and a rectilinear grid's lines are taken relative to it, subtracted in double precision before they become dtype, so points far from zero keep their precision.

A vtkOverlappingAMR becomes an AMRDataSet: each block (an image or a rectilinear grid) converted as above, placed in its level by its vtkAMRBox, the levels' refinement ratios kept; the refined cells are marked afresh (VTK's vtkGhostType bits come in and are recomputed). Blocks with ghost layers are refused.

A vtkPartitionedDataSet becomes a PartitionedDataSet, each piece converted as above (null pieces skipped), and a vtkPartitionedDataSetCollection a dict of them, keyed by each partitioned dataset's vtkCompositeDataSet.NAME() ("block<i>" where it has none, "<name>#<i>" for a name already taken).

VTK's DATA_TIME_STEP becomes the step's time (time), its unit unknown.

Source code in packages/tack-vis/src/tack/interop/vtk.py

def vtk_to_dataset(grid, dtype=tack.f32, polyhedral=None, orient=False, id_dtype=None,
                   origin=None):
    """Copy a ``vtkUnstructuredGrid``, ``vtkPolyData``, ``vtkStructuredGrid``,
    ``vtkRectilinearGrid`` or ``vtkImageData`` into a ``tack.dataset.DataSet``: an image's points as
    ``UniformCoordinates`` (with its direction matrix), a rectilinear grid's as a
    ``CartesianProduct``, a structured grid's as its explicit points.

    Points and floating-point arrays become fields of ``dtype``; integer
    arrays keep their type. Named point data become ``H1`` fields and named
    cell data ``Constant`` fields (a cell array named as a point array is
    ``"<name> (cells)"``, with a warning); a ``vtkGhostType`` array becomes the ghost marks
    (``data.ghosts``), not a field; a constant or affine array (VTK's implicit
    arrays) stays implicit, as a ``ConstantArray`` or ``CountingArray``. An
    unstructured grid's cells must all be linear shapes -- unless it holds
    polyhedra (or ``polyhedral`` is true), when it becomes a
    ``PolyhedralTopology``, every cell a polyhedron of its faces, polylines,
    when it becomes a ``PolylineTopology`` (its lines two-point polylines, and
    no other cells allowed), or polygons (``VTK_POLYGON``), when it becomes a
    ``PolygonalTopology`` (its triangles and quads polygons too, and no other
    cells allowed). Its types, offsets and connectivity stay implicit where
    VTK's are (a cell type given as a number, affine offsets).
    A ``vtkPolyData`` is taken as the unstructured grid of its cells in VTK's
    order -- verts, lines, polygons (of three or four points triangles and
    quads), and triangle strips, each strip as its triangles, which take its
    cell data -- so the same rules hold: a surface of triangles and quads is
    an ``UnstructuredTopology``, one with larger polygons a ``PolygonalTopology``,
    lines a ``PolylineTopology``, and polylines or polygons beside other cells
    are refused.
    With ``orient``, inconsistently wound polyhedra are repaired first
    (``PolyhedralTopology.from_cell_faces``); without it they are refused.
    VTK's ids are 64-bit; the topology keeps them in ``id_dtype``, by default
    ``i32`` unless the grid is too large for it (``tack.dataset.ids``).

    ``origin`` (three numbers) becomes the dataset's ``origin``: VTK's points,
    an image's origin and a rectilinear grid's lines are taken relative to it,
    subtracted in double precision before they become ``dtype``, so points far
    from zero keep their precision.

    A ``vtkOverlappingAMR`` becomes an ``AMRDataSet``: each block (an image or a
    rectilinear grid) converted as above, placed in its level by its
    ``vtkAMRBox``, the levels' refinement ratios kept; the refined cells are
    marked afresh (VTK's ``vtkGhostType`` bits come in and are recomputed).
    Blocks with ghost layers are refused.

    A ``vtkPartitionedDataSet`` becomes a ``PartitionedDataSet``, each piece
    converted as above (null pieces skipped), and a
    ``vtkPartitionedDataSetCollection`` a ``dict`` of them, keyed by each
    partitioned dataset's ``vtkCompositeDataSet.NAME()`` (``"block<i>"`` where
    it has none, ``"<name>#<i>"`` for a name already taken).

    VTK's ``DATA_TIME_STEP`` becomes the step's time (``time``), its unit
    unknown."""
    out = _vtk_to_dataset(grid, dtype=dtype, polyhedral=polyhedral, orient=orient,
                          id_dtype=id_dtype, origin=origin)
    if grid.GetInformation().Has(_time_key()):
        _stamp(out, grid.GetInformation().Get(_time_key()))
    return out

dataset_to_vtk

dataset_to_vtk(data)

Copy a tack.dataset.DataSet into a new VTK dataset; an AMRDataSet of images and rectilinear grids, its ratios the same along every axis, into a vtkOverlappingAMR (its marks as vtkGhostType arrays), any other into a vtkPartitionedDataSetCollection of a partitioned dataset per level ("level<l>"), VTK's AMR holding no curvilinear grids; a PartitionedDataSet into a vtkPartitionedDataSet of its pieces, and a dict of names to PartitionedDataSets into a vtkPartitionedDataSetCollection whose partitioned datasets carry their names (vtkCompositeDataSet.NAME()).

A structured topology over uniform coordinates gives a vtkImageData (with the direction matrix), over rectilinear ones a vtkRectilinearGrid, over explicit ones (a warped or transformed grid) a vtkStructuredGrid; an unstructured one a vtkUnstructuredGrid, a polygonal or polyline one a vtkUnstructuredGrid of VTK_POLYGON or VTK_POLY_LINE cells, and a polyhedral one a vtkUnstructuredGrid of VTK_POLYHEDRON cells, each keeping its own copy of its faces, outward for it (VTK's layout has no orientation bit, so a shared face cannot be stored once). H1 fields and values on points become point data -- an order-2 field its values at the points, the grid being linear -- and Constant fields and values on cells cell data. L2 fields have no VTK array to go to without duplicating points, and fields on faces or edges none at all; they are left out. The geometry (data.geometry) becomes the points, so it must be H1. Implicit types, offsets and connectivity, and constant and counting field values, become VTK's constant and affine arrays. Ghost marks become vtkGhostType arrays.

The step's time (data.time) becomes VTK's DATA_TIME_STEP; VTK has no unit for it.

Source code in packages/tack-vis/src/tack/interop/vtk.py

def dataset_to_vtk(data):
    """Copy a ``tack.dataset.DataSet`` into a new VTK dataset; an ``AMRDataSet`` of
    images and rectilinear grids, its ratios the same along every axis, into a
    ``vtkOverlappingAMR`` (its marks as ``vtkGhostType`` arrays), any other into
    a ``vtkPartitionedDataSetCollection`` of a partitioned dataset per level
    (``"level<l>"``), VTK's AMR holding no curvilinear grids; a
    ``PartitionedDataSet`` into a ``vtkPartitionedDataSet`` of its pieces, and a
    ``dict`` of names to
    ``PartitionedDataSet``s into a ``vtkPartitionedDataSetCollection`` whose
    partitioned datasets carry their names (``vtkCompositeDataSet.NAME()``).

    A structured topology over uniform coordinates gives a ``vtkImageData`` (with
    the direction matrix), over rectilinear ones a ``vtkRectilinearGrid``, over
    explicit ones (a warped or transformed grid) a ``vtkStructuredGrid``; an
    unstructured one a ``vtkUnstructuredGrid``, a polygonal or polyline one a
    ``vtkUnstructuredGrid`` of ``VTK_POLYGON`` or ``VTK_POLY_LINE`` cells, and a
    polyhedral one a ``vtkUnstructuredGrid`` of ``VTK_POLYHEDRON`` cells, each
    keeping its own copy of its faces, outward for it (VTK's layout has no
    orientation bit, so a shared face cannot be stored once).
    ``H1`` fields and values on points become point data -- an order-2 field
    its values at the points, the grid being linear -- and ``Constant`` fields
    and values on cells cell data. ``L2`` fields have no VTK array to go
    to without duplicating points, and fields on faces or edges none at
    all; they are left out. The geometry (``data.geometry``) becomes the
    points, so it must be ``H1``. Implicit types, offsets and connectivity,
    and constant and counting field values, become VTK's constant and affine
    arrays. Ghost marks become ``vtkGhostType`` arrays.

    The step's time (``data.time``) becomes VTK's ``DATA_TIME_STEP``; VTK has
    no unit for it."""
    out = _dataset_to_vtk(data)
    time = _time_of(data)
    if time is not None:
        from vtkmodules.vtkCommonDataModel import vtkDataObject

        out.GetInformation().Set(vtkDataObject.DATA_TIME_STEP(), time)
    return out

tack.interop.mfem

MFEM meshes and grid functions as tack.dataset datasets, through PyMFEM.

mfem_to_dataset(mesh, fields) copies an mfem.ser.Mesh and named GridFunctions on it. The mesh's elements are the topology, and its element attributes an i32 Constant field, "attribute".

MFEM lists every element's corners in VTK's order, the prism included: both put the triangle (0, 1, 2) with its normal toward (3, 4, 5), and a positively oriented MFEM prism has a positive Jacobian as a VTK wedge, as VTK's own wedges do. (MFEM's VTK writer swaps a prism's corners 1 and 2, and 4 and 5 -- PrismMap in mesh/vtk.cpp -- which turns such a prism inside out against VTK's documented and implemented convention.)

Each grid function becomes a field in the space its collection matches:

• H1 order 1 (H1_*_P1): an H1 field, one value per vertex, which is how MFEM numbers an order-1 H1 space's DOFs;
• L2 order 0: a Constant field;
• L2 order 1: an L2 field of each element's values at its corners, where MFEM evaluates its own basis (its DOFs sit at interior Gauss points, not corners). On tetrahedra, hexahedra and prisms the two spaces are the same, so this is exact; on pyramids MFEM's order-1 L2 space is not the collapsed trilinear one, and the corners interpolate it;
• H1 order 2: an H1 order-2 field, MFEM's function evaluated at each element's quadratic nodes (corners, edge middles, quad-face middles, a hexahedron's center). MFEM's order-2 nodes are these same points, so it is exact; pyramids, whose MFEM basis is another, are refused.

A curved mesh -- nodes of order 2 -- gets an order-2 geometry the same way. Vector grid functions (vdim > 1) become fields of vectors, in either ordering. Higher orders are refused until their bases exist. PyMFEM (pip install mfem) is imported on first use.

mfem_field

mfem_field(data, mesh, gf, dtype=tack.f32)

gf, an MFEM grid function on mesh, as a field of data (the dataset mfem_to_dataset made from mesh).

Source code in packages/tack-vis/src/tack/interop/mfem.py

def mfem_field(data, mesh, gf, dtype=tack.f32):
    """``gf``, an MFEM grid function on ``mesh``, as a field of ``data`` (the dataset
    ``mfem_to_dataset`` made from ``mesh``)."""
    from tack.dataset import H1, L2, Constant, Field

    fes = gf.FESpace()
    family, order = _collection(fes)
    if family == "H1" and order == 1:
        if fes.GetNDofs() != mesh.GetNV():
            raise ValueError("an order-1 H1 space has a DOF per vertex")
        return Field(H1(data), _field_values(_by_dof(fes, gf.GetDataArray()), dtype))
    if family == "L2" and order == 0:
        return Field(Constant(data), _field_values(_by_dof(fes, gf.GetDataArray()), dtype))
    if family == "L2" and order == 1:
        return Field(L2(data), _field_values(_corner_values(mesh, gf), dtype))
    if family == "H1" and order == 2:
        space, values = _quadratic_values(data, mesh, gf)
        return Field(space, _field_values(values, dtype))
    raise NotImplementedError(f"MFEM {family} order {order} has no space here yet")

mfem_to_dataset

mfem_to_dataset(
    mesh, fields=None, dtype=tack.f32, id_dtype=None
)

A tack.dataset.DataSet of an MFEM mesh and named grid functions on it.

The geometry is the vertex positions (z = 0 for a 2D mesh), or for a curved mesh -- nodes of order 2 -- an H1 order-2 field of positions; nodes of higher order are refused. Element attributes become fields["attribute"]. id_dtype is the topology's (tack.dataset.ids).

Source code in packages/tack-vis/src/tack/interop/mfem.py

def mfem_to_dataset(mesh, fields=None, dtype=tack.f32, id_dtype=None):
    """A ``tack.dataset.DataSet`` of an MFEM mesh and named grid functions on it.

    The geometry is the vertex positions (z = 0 for a 2D mesh), or for a curved
    mesh -- nodes of order 2 -- an H1 order-2 field of positions; nodes of
    higher order are refused. Element attributes become ``fields["attribute"]``.
    ``id_dtype`` is the topology's (``tack.dataset.ids``).
    """
    from tack.dataset import Constant, DataSet, Field

    nodes = mesh.GetNodes()
    curved = nodes is not None and _collection(nodes.FESpace())[1] > 1
    if curved and _collection(nodes.FESpace())[1] > 2:
        raise NotImplementedError("curved MFEM meshes of order above 2 are not supported "
                                  "yet")
    topology = _topology(mesh, id_dtype)
    if curved:
        positions = np.asarray(mesh.GetVertexArray(), dtype=float).reshape(
            mesh.GetNV(), mesh.SpaceDimension())
    elif nodes is not None:
        # Order-1 nodes are the vertices' positions, and may have moved since
        # the vertex array was set (Mesh.Transform moves the nodes).
        positions = _by_dof(nodes.FESpace(), nodes.GetDataArray())
    else:
        positions = np.asarray(mesh.GetVertexArray(), dtype=float).reshape(
            mesh.GetNV(), mesh.SpaceDimension())
    if positions.shape[1] < 3:
        positions = np.hstack([positions, np.zeros((positions.shape[0],
                                                    3 - positions.shape[1]))])
    data = DataSet(topology, positions, dtype=dtype)
    if curved:
        # The corners first, to have a dataset to number the quadratic nodes on;
        # then the curved geometry, MFEM's nodes at those nodes.
        space, values = _quadratic_values(data, mesh, nodes)
        data.geometry = Field(space, _field_values(values, dtype))
    attributes = np.array([mesh.GetAttribute(e) for e in range(mesh.GetNE())], np.int32)
    data.fields["attribute"] = Field(Constant(data), _field_values(attributes, tack.i32))
    for name, gf in (fields or {}).items():
        data.fields[name] = mfem_field(data, mesh, gf, dtype=dtype)
    return data

tack.interop.warp

Arrays and datasets between Tack and NVIDIA Warp, without copies where they agree.

With the Warp device (tack.device("warp")) a Tack array's memory is a wp.array on Warp's device, and Tack launches on that device's current Warp stream, the stream a caller's wp.launch uses. So a program can run its own Warp kernels and Tack's filters on the same arrays:

tack.set_default_device("warp")
wp.launch(my_kernel, n, inputs=[f])
data.fields["f"] = td.Field(td.H1(data), from_warp(f))
iso = alg.contour(data, "f", 0.5)
points = to_warp(iso.geometry.values)          # a wp.array of wp.vec3

from_warp and to_warp map Warp's vector and matrix types to Tack's (wp.vec3 and tack.vec(3, tack.f32), wp.mat33 and tack.mat(3, 3, tack.f32)), keep the other side's memory alive for as long as the result lives, and do not mark the memory as handed outside Tack: on one device and stream, Tack's launches need not complete before returning. Under wp.ScopedStream Tack follows the scope's stream too.

Warp arrays have no read-only flag: a read-only Tack array given to Warp can be written there.

to_simdata and from_simdata exchange datasets with warp-simdata (NVIDIA's simulation-data library on Warp, an optional import): image data, structured grids, polydata and unstructured grids (polyhedra included) against Tack's topologies, NODE and ELEMENT fields against point and cell fields, collections against partitioned datasets. Most of warp-simdata's other models come across as they are, sharing their arrays (see from_simdata); the rest through its to_vtk_data_model.

from_warp

from_warp(array, *, writable=True)

array, a contiguous wp.array, as a Tack array on the Warp device, sharing its memory: a write through either is seen by the other.

The element type follows Warp's (wp.vec3 gives tack.vec(3, tack.f32)) and the shape is the array's. The Tack array keeps array alive. array must be on the device Tack's Warp device drives (Warp's default device).

Source code in packages/tack-vis/src/tack/interop/warp.py

def from_warp(array, *, writable=True):
    """``array``, a contiguous ``wp.array``, as a Tack array on the Warp device,
    sharing its memory: a write through either is seen by the other.

    The element type follows Warp's (``wp.vec3`` gives ``tack.vec(3,
    tack.f32)``) and the shape is the array's. The Tack array keeps
    ``array`` alive. ``array`` must be on the device Tack's Warp device
    drives (Warp's default device).
    """
    from tack.runtime.warp_backend import WarpBuffer
    wp = _warp()
    if not isinstance(array, wp.array):
        raise TypeError(f"from_warp takes a wp.array, not {type(array).__name__}")
    if not array.is_contiguous:
        raise ValueError("from_warp wraps a contiguous wp.array; copy it with "
                         "wp.clone first")
    backend = _warp_backend()
    if str(array.device) != str(backend.wp_device):
        raise ValueError(f"the array is on Warp's {array.device}, but Tack's Warp device "
                         f"drives {backend.wp_device} (Warp's default device)")
    dtype = _tack_dtype(array.dtype)
    shape = tuple(array.shape)
    scalar = dtype.scalar_dtype if isinstance(dtype, ElementType) else dtype
    count = math.prod(shape) * (dtype.size if isinstance(dtype, ElementType) else 1)
    storage = (count,) if isinstance(dtype, ElementType) else shape
    if count:
        flat = wp.array(ptr=array.ptr, dtype=_warp_dtype(scalar), shape=(count,),
                        capacity=count * np.dtype(scalar.numpy_dtype).itemsize,
                        device=array.device)
    else:
        flat = None                     # nothing to share; WarpBuffer makes its own
    buffer = WarpBuffer(backend, scalar, storage, flat)
    buffer._owner = array               # the memory lives as long as the Tack array
    return _array_view(dtype, shape, buffer, writable=writable, device=backend)

to_warp

to_warp(array)

array as a wp.array sharing its memory, with Warp's element type (tack.vec(3, tack.f32) gives wp.vec3) at its logical shape.

An array on the Warp device is its own memory, kept alive by the result. An array on another Tack device goes through DLPack, which completes its queued launches. An implicit array (tack.arrays: ConstantArray, UniformCoordinates, ...) has no memory, and is materialized on its device first. Structure-of-arrays storage has no single wp.array; give its components one at a time.

Source code in packages/tack-vis/src/tack/interop/warp.py

def to_warp(array):
    """``array`` as a ``wp.array`` sharing its memory, with Warp's element type
    (``tack.vec(3, tack.f32)`` gives ``wp.vec3``) at its logical shape.

    An array on the Warp device is its own memory, kept alive by the result.
    An array on another Tack device goes through DLPack, which completes its
    queued launches. An implicit array (``tack.arrays``: ``ConstantArray``,
    ``UniformCoordinates``, ...) has no memory, and is materialized on its
    device first. Structure-of-arrays storage has no single ``wp.array``;
    give its components one at a time.
    """
    wp = _warp()
    if isinstance(array, _Implicit):
        array = materialize(array)
    if isinstance(array, SoaArray):
        raise TypeError("a structure-of-arrays array has no single wp.array; "
                        "give its components (array.components) one at a time")
    if not isinstance(array, Array):
        raise TypeError(f"to_warp takes a Tack array, not {type(array).__name__}")
    shape, dtype, nbytes = tuple(array.shape), array.dtype, array._buffer.nbytes
    backend = array._device
    if backend.name == "warp":
        address, device = array._buffer.address, backend.wp_device
    else:
        source = wp.from_dlpack(array)          # Tack's export completes its launches
        address, device = source.ptr, source.device
        array = (array, source)                 # keep both alive
    if not math.prod(shape):
        return wp.empty(shape, dtype=_warp_dtype(dtype), device=device)

    def keep(_ptr, _size, _held=array):
        """Warp calls this when the result is freed; it holds the Tack array till then."""

    return wp.array(ptr=address, dtype=_warp_dtype(dtype), shape=shape, capacity=nbytes,
                    device=device, deleter=keep)

to_simdata

to_simdata(data)

A Tack DataSet (or PartitionedDataSet) as a warp-simdata dataset (or DatasetCollection) on Warp's device.

A uniform grid becomes image data (an origin, a spacing and an extent; a rotated one becomes a structured grid), any other structured grid a structured grid, polygons and polylines polydata, and every other topology an unstructured grid, polyhedra with their faces. Point fields become NODE fields and cell fields ELEMENT ones. Arrays already f32 positions, i32 ids and u8 types on the Warp device are shared; others are copied, positions as f32 world positions (Tack's origin added, so far coordinates lose the precision the origin kept). Ghost marks, global ids, sets, units and time have no counterpart there and are left behind.

Source code in packages/tack-vis/src/tack/interop/warp.py

def to_simdata(data):
    """A Tack ``DataSet`` (or ``PartitionedDataSet``) as a warp-simdata dataset (or
    ``DatasetCollection``) on Warp's device.

    A uniform grid becomes image data (an origin, a spacing and an extent; a
    rotated one becomes a structured grid), any other structured grid a
    structured grid, polygons and polylines polydata, and every other topology
    an unstructured grid, polyhedra with their faces. Point fields become NODE
    fields and cell fields ELEMENT ones. Arrays already f32 positions, i32 ids
    and u8 types on the Warp device are shared; others are copied, positions
    as f32 world positions (Tack's origin added, so far coordinates lose the
    precision the origin kept). Ghost marks, global ids, sets, units and time
    have no counterpart there and are left behind.
    """
    sd = _simdata()
    from warp_simdata.data_models.vtk import (
        image_data,
        polydata,
        structured_grid,
        unstructured_grid,
    )

    from tack.arrays import UniformCoordinates
    from tack.dataset import (
        PartitionedDataSet,
        PolygonalTopology,
        PolyhedralTopology,
        PolylineTopology,
        StructuredTopology,
    )

    wp = _warp()
    if isinstance(data, PartitionedDataSet):
        from warp_simdata.data_models.vtk.conversion import to_common_vtk_data_model
        return sd.DatasetCollection.from_datasets(
            to_common_vtk_data_model([to_simdata(piece) for piece in data]))
    device = _wp_device(data)
    topology, coordinates = data.topology, data.geometry.values
    if isinstance(topology, StructuredTopology):
        dims = list(topology.point_dims) + [1] * (3 - len(topology.point_dims))
        lo = list(topology.offset) + [0] * (3 - len(topology.offset))
        hi = [o + d - 1 for o, d in zip(lo, dims)]
        if isinstance(coordinates, UniformCoordinates) and np.allclose(
                coordinates._matrix, np.eye(3)):
            spacing = np.asarray(coordinates._spacing, np.float64)
            start = np.asarray(coordinates._origin, np.float64) + np.asarray(data.origin)
            out = image_data.create_dataset(
                origin=wp.vec3f(*(start - np.asarray(lo) * spacing)),   # node ijk at 0, as VTK
                spacing=wp.vec3f(*spacing), extent_min=wp.vec3i(*lo),
                extent_max=wp.vec3i(*hi), device=str(device))
        else:
            out = structured_grid.create_dataset(_points(data, device), wp.vec3i(*lo),
                                                 wp.vec3i(*hi))
    elif isinstance(topology, (PolygonalTopology, PolylineTopology)):
        polygons = isinstance(topology, PolygonalTopology)
        offsets, points = topology.loops() if polygons else topology.lines()
        rings = {"polys_offsets" if polygons else "lines_offsets":
                 _shared_or_copied(offsets, wp.int32, device),
                 "polys_connectivity" if polygons else "lines_connectivity":
                 _shared_or_copied(points, wp.int32, device)}
        out = polydata.create_dataset(_points(data, device), **rings)
    elif isinstance(topology, PolyhedralTopology):
        # Faces shared between cells: warp-simdata finds each cell's outward
        # winding itself (its face orientations).
        point_offsets, point_ids = topology.cell_points()
        out = unstructured_grid.create_dataset(
            _points(data, device),
            wp.full(topology.num_cells, 42, dtype=wp.uint8, device=device),
            _shared_or_copied(point_offsets, wp.int32, device),
            _shared_or_copied(point_ids, wp.int32, device),
            faces_offsets=_shared_or_copied(topology.face_offsets, wp.int32, device),
            faces_connectivity=_shared_or_copied(topology.face_points, wp.int32, device),
            face_locations_offsets=_shared_or_copied(topology.cell_offsets, wp.int32, device),
            face_locations_connectivity=_shared_or_copied(topology.cell_faces, wp.int32,
                                                          device))
    else:
        out = unstructured_grid.create_dataset(
            _points(data, device),
            _shared_or_copied(topology.types, wp.uint8, device),
            _shared_or_copied(topology.offsets, wp.int32, device),
            _shared_or_copied(topology.connectivity, wp.int32, device))
    _add_fields(data, out, device)
    return out

from_simdata

from_simdata(ds)

A warp-simdata dataset (or DatasetCollection) as a Tack DataSet (or PartitionedDataSet) on the Warp device, sharing its arrays.

Point clouds, curves, surface meshes, OpenFOAM polyMesh and boundary patches, CGNS sections of linear elements (one shape, MIXED, NGON_n and NFACE_n), EnSight parts of one kind, corner-point grids (inactive cells hidden) and FLASH axisymmetric data (an AMRDataSet of revolved boxes) are taken as they are: their arrays become Tack's, CGNS's 1-based and MIXED connectivities read lazily (see "models imported as they are"). Every other model goes through warp-simdata's to_vtk_data_model first. Image data becomes a uniform grid (offset from its extent), a structured grid a structured topology over its points, polydata of polygons or polylines a polygonal or polyline topology (vertices alone become vertex cells; a mix of kinds is refused: split it first), and an unstructured grid an unstructured topology, or a polyhedral one when it has polyhedra. NODE fields become H1 fields and ELEMENT fields Constant ones. The dataset must be on the device Tack's Warp device drives.

Source code in packages/tack-vis/src/tack/interop/warp.py

def from_simdata(ds):
    """A warp-simdata dataset (or ``DatasetCollection``) as a Tack ``DataSet`` (or
    ``PartitionedDataSet``) on the Warp device, sharing its arrays.

    Point clouds, curves, surface meshes, OpenFOAM polyMesh and boundary
    patches, CGNS sections of linear elements (one shape, MIXED, NGON_n and
    NFACE_n), EnSight parts of one kind, corner-point grids (inactive cells
    hidden) and FLASH axisymmetric data (an ``AMRDataSet`` of revolved boxes)
    are taken as they are: their arrays become Tack's, CGNS's 1-based and MIXED
    connectivities read lazily (see "models imported as they are"). Every
    other model goes through warp-simdata's ``to_vtk_data_model`` first. Image data becomes a uniform grid (``offset`` from its extent), a structured
    grid a structured topology over its points, polydata of polygons or
    polylines a polygonal or polyline topology (vertices alone become vertex
    cells; a mix of kinds is refused: split it first), and an unstructured grid
    an unstructured topology, or a polyhedral one when it has polyhedra. NODE
    fields become H1 fields and ELEMENT fields Constant ones. The dataset must
    be on the device Tack's Warp device drives.
    """
    sd = _simdata()
    from warp_simdata.data_models.vtk import (
        image_data,
        polydata,
        structured_grid,
        unstructured_grid,
    )
    from warp_simdata.data_models.vtk.conversion import to_vtk_data_model

    from tack.arrays import ConstantArray, UniformCoordinates
    from tack.dataset import (
        PartitionedDataSet,
        PolygonalTopology,
        PolyhedralTopology,
        PolylineTopology,
        UnstructuredTopology,
    )
    from tack.dataset.dataset import grid_topology

    if isinstance(ds, sd.DatasetCollection):
        pieces = ds.datasets
        if any(not piece.has_field(name) for name in ds.get_field_names() for piece in pieces):
            # Fields of the collection itself: its VTK conversion hands them out.
            pieces = to_vtk_data_model(ds).datasets
        return PartitionedDataSet(from_simdata(piece) for piece in pieces)
    device = tack.device("warp")
    from warp_simdata.data_models.flash import axisymmetric as flash_model
    if ds.data_model is flash_model.DataModel:
        return _flash_amr(ds, device)
    with device:
        direct = _direct_topology(ds, device)
        if direct is None:
            more = _more_direct(ds, device)
            if more is not None:
                topology, geometry, ghosts = more
                return _with_fields(ds, topology, geometry, device, ghosts)
    if direct is not None:
        return _with_fields(ds, *direct, device)
    ds = to_vtk_data_model(ds)
    handle, kind = ds.handle, ds.data_model.DatasetHandle
    if kind is image_data.DatasetHandle:
        lo = np.asarray(handle.extent_min, np.int64)
        dims = np.asarray(handle.extent_max, np.int64) - lo + 1
        spacing = np.asarray(handle.spacing, np.float64)
        start = np.asarray(handle.origin, np.float64) + lo * spacing
        topology = grid_topology(dims.tolist(), offset=lo.tolist(), device=device)
        geometry = UniformCoordinates(dims.tolist(), start, spacing, np.eye(3), device=device)
    elif kind is structured_grid.DatasetHandle:
        lo = np.asarray(handle.extent_min, np.int64)
        dims = np.asarray(handle.extent_max, np.int64) - lo + 1
        topology = grid_topology(dims.tolist(), offset=lo.tolist(), device=device)
        geometry = from_warp(handle.points)
    elif kind is polydata.DatasetHandle:
        def present(offsets):
            return offsets is not None and offsets.shape[0] > 1
        kinds = [k for k in ("verts", "lines", "polys") if present(getattr(handle, f"{k}_offsets"))]
        if len(kinds) > 1:
            raise ValueError(f"polydata of {' and '.join(kinds)} together has no one Tack "
                             "topology: split it by kind first")
        n = ds.get_num_nodes()
        kind = kinds[0] if kinds else "verts"
        offsets = getattr(handle, f"{kind}_offsets")
        points = getattr(handle, f"{kind}_connectivity")
        offsets = from_warp(offsets) if present(offsets) else np.zeros(1, np.int32)
        points = from_warp(points) if present(getattr(handle, f"{kind}_offsets")) \
            else np.zeros(0, np.int32)
        with device:                        # host arrays among them go to the Warp device
            if kind == "polys":
                topology = PolygonalTopology(offsets, points, num_points=n)
            elif kind == "lines":
                topology = PolylineTopology(offsets, points, num_points=n)
            else:
                count = (len(offsets) if isinstance(offsets, np.ndarray) else offsets.shape[0]) - 1
                if points.shape[0] == count:
                    # One point each: every cell a VTK_VERTEX, an implicit array stored
                    # nowhere, as Tack's own point clouds' (uniform_cells).
                    types = ConstantArray(1, count, dtype=tack.u8)
                else:
                    # Poly-vertices among them (VTK_POLY_VERTEX, 2): a type per cell.
                    counts = np.diff(offsets.numpy() if isinstance(offsets, Array) else offsets)
                    types = np.where(counts == 1, 1, 2).astype(np.uint8)
                topology = UnstructuredTopology(types, offsets, points, num_points=n)
        geometry = from_warp(handle.points)
    elif kind is unstructured_grid.DatasetHandle:
        n = ds.get_num_nodes()
        types = from_warp(handle.cell_types)
        present = _present_types(types)
        if 42 in present:
            with device:                    # its new arrays on the Warp device
                topology = _polyhedra_from_simdata(handle, PolyhedralTopology, n)
        else:
            topology = UnstructuredTopology(types, from_warp(handle.cell_offsets),
                                            from_warp(handle.cell_connectivity), num_points=n)
        geometry = from_warp(handle.points)
    else:
        raise TypeError(f"no Tack topology for warp-simdata's {kind.__name__}")
    return _with_fields(ds, topology, geometry, device)