Skip to content

Algorithms

Every algorithm of tack.dataset, by what it does, with its parameters. Each group opens with a table: what each algorithm gives, what it takes (a DataSet; pieces, a PartitionedDataSet, over MPI ranks or not; an AMRDataSet; an OctreeDataSet), and its nearest counterpart in VTK and in Viskores. Below the table, each algorithm's signature and description.

import tack.dataset as td

hot = td.threshold(data, "temperature", 300.0, 400.0)
grad = td.gradients(data, "pressure", at="points")
s = td.statistics(data, "pressure")

Every algorithm is in tack.dataset.algorithms, one module each, and tack.dataset re-exports them all. They share conventions, set out in the Dataset API reference's Algorithms and filters:

  • the dataset comes first, and a field argument is a Field or the name of one of the dataset's fields;
  • an algorithm that makes a new dataset takes fields=: which fields to carry (all by default, a name, a list of names, or [] for none), carried by the rules in how fields carry through filters;
  • an algorithm that computes a field returns a Field, several fields a named tuple of them, a summary Python numbers, a selection a Set;
  • nothing changes its input, except where an entry says "in place";
  • every algorithm runs, and allocates, on its dataset's device; algorithms over pieces run each piece on its own device;
  • hidden cells are never processed, and ghost cells and points are left out of summaries;
  • a parameter that is a length, a time or a field's value (a radius, a step size, a threshold's range) also takes a quantity ("10 cm", (10, "cm")), converted to the dataset's or the field's unit; a plain number is taken as it is (Units).

Which algorithms work on which topologies (polyhedra, polygons, polylines) is in the Dataset API reference's support table. The classes the algorithms take and give -- datasets, fields, spaces, topologies, sets, implicit functions, locators -- are in Data Model.

Contents

Cells kept whole

New datasets of some of the input's cells, unchanged, with only the points they use.

Algorithm Gives Takes VTK Viskores
threshold cells whose field lies in a range dataset vtkThreshold Threshold
extract_cells cells by id, array or cell set dataset vtkExtractCells —
extract_geometry cells inside (or outside, or cut by) an implicit function dataset vtkExtractGeometry ExtractGeometry
extract_extent a box of a structured grid, every n-th point dataset vtkExtractGrid ExtractStructured
mask every n-th cell dataset — Mask

threshold 

threshold(
    data, field, lower, upper, all_points=True, fields=None
)

The cells whose field lies in [lower, upper], with only the points they use.

field is a scalar field or its name. A cell field (Constant, values on cells) keeps a cell by its value; a field with a basis (H1, L2 of any orders) or values on points by its values at the cell's corners: all of them in range, or with all_points=False any (VTK's AllScalars). The result's cells keep their shapes and order; its points are the used ones, renumbered in order. Every field it can carry comes along: point fields gathered to the kept points, cell fields to the kept cells, L2 fields (and an L2 geometry) block by block, keeping each cell's order. Fields on faces and edges, which the new topology derives afresh, are left behind.

On a polyhedral topology the kept cells keep their faces, stored once and numbered in their old order, so values on faces come along. A face whose side-0 cell is dropped is turned around for its remaining cell, and its oriented values (Values(..., "faces", oriented=True)) are negated. A polygonal topology's kept polygons keep their loops, and a polyline topology's polylines their chains; their facets' values, whose numbering is derived afresh, are left behind.

Source code in packages/tack-vis/src/tack/dataset/algorithms/threshold.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
@on_inputs_device
def threshold(data, field, lower, upper, all_points=True, fields=None):
    """The cells whose ``field`` lies in ``[lower, upper]``, with only the points they use.

    ``field`` is a scalar field or its name. A cell field (``Constant``, values
    on cells) keeps a cell by its value; a field with a basis (``H1``, ``L2``
    of any orders) or values on points by its values at the cell's corners:
    all of them in range, or with ``all_points=False`` any (VTK's
    ``AllScalars``). The result's cells keep their shapes and order; its points
    are the used ones, renumbered in order. Every field it can carry comes
    along: point fields gathered to the kept points, cell fields to the kept
    cells, ``L2`` fields (and an ``L2`` geometry) block by block, keeping each
    cell's order. Fields on faces and edges, which the new topology derives
    afresh, are left behind.

    On a polyhedral topology the kept cells keep their faces, stored once and
    numbered in their old order, so values on faces come along. A face whose
    side-0 cell is dropped is turned around for its remaining cell, and its
    oriented values (``Values(..., "faces", oriented=True)``) are negated. A
    polygonal topology's kept polygons keep their loops, and a polyline
    topology's polylines their chains; their facets' values, whose numbering is
    derived afresh, are left behind.
    """
    lower, upper = (_in_unit_of(_field(data, field), v, "threshold's bounds")
                    for v in (lower, upper))
    field = _field(data, field)
    if arrays.width_of(field.values):
        raise TypeError("threshold needs a scalar field")
    field = _with_basis(data, field)
    space = field.space
    if not space.interpolated or space.on not in ("cells", "points"):
        raise TypeError(f"threshold needs a field on the points or cells, not {space!r}")
    return _keep_cells(data, _threshold_flags(data, field, lower, upper, all_points), fields)

extract_cells 

extract_cells(data, cells, fields=None)

The cells with the given ids (a host array, an integer field, or a cell Set or its name in data.sets), with only the points they use; cells keep their order, not the order given.

Source code in packages/tack-vis/src/tack/dataset/algorithms/extract_cells.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
@on_inputs_device
def extract_cells(data, cells, fields=None):
    """The cells with the given ids (a host array, an integer field, or a cell
    ``Set`` or its name in ``data.sets``), with only the points they use; cells
    keep their order, not the order given."""
    from tack.dataset.sets import Set

    if isinstance(cells, str):
        cells = data.sets[cells]
    if isinstance(cells, Set):
        if cells.on != "cells" or cells.topology is not data.topology:
            raise ValueError(f"extract_cells takes a set of this dataset's cells, not a "
                             f"set of {cells.on}" + ("" if cells.topology is data.topology
                                                    else " of another topology"))
        cells = cells.ids
    if not hasattr(cells, "numpy"):
        values = np.asarray(cells).reshape(-1)
        if values.size and not np.issubdtype(values.dtype, np.integer):
            raise TypeError("cell ids must be integers")
        cells = values
    chosen = ids.as_ids(cells, data.id_dtype)
    keep = tack.zeros((data.num_cells,), tack.i32)
    if chosen.shape[0]:
        ids.check_bounds(chosen, data.num_cells, "cell")
        _flag_ids(chosen, keep)
    return _keep_cells(data, keep, fields)

extract_geometry 

extract_geometry(
    data, function, inside=True, boundary=False, fields=None
)

The whole cells inside function's region (all their points at or below zero), or, with inside=False, outside it (all above); with boundary, also the cells it cuts (some point on the kept side). Cells are not split -- that is clip. Points are compacted and fields carried, as by threshold.

Source code in packages/tack-vis/src/tack/dataset/algorithms/extract_geometry.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
@on_inputs_device
def extract_geometry(data, function, inside=True, boundary=False, fields=None):
    """The whole cells inside ``function``'s region (all their points at or below zero),
    or, with ``inside=False``, outside it (all above); with ``boundary``, also the
    cells it cuts (some point on the kept side). Cells are not split -- that is
    clip. Points are compacted and fields carried, as by ``threshold``."""

    values = implicit_values(data, function)
    if inside:
        return threshold(data, values, -np.inf, 0.0, all_points=not boundary, fields=fields)
    # Outside is what is not inside, exactly: all points above zero is no point at
    # or below it, and some point above is not all at or below.
    keep = _threshold_flags(data, values, -np.inf, 0.0, all_points=boundary)
    _complement(keep)
    return _keep_cells(data, keep, fields)

extract_extent 

extract_extent(data, extent, rate=1, fields=None)

The points of data, a structured grid, in the box extent -- one (low, high) pair of point indices per axis, both included, or VTK's flat (i0, i1, j0, j1, k0, k1), in the whole grid's indices (as VTK's extents count them, from the topology's offset), clipped to the grid -- every rate-th along each axis (a number, or one per axis): a structured grid of the same kind, an image's coordinates still uniform, a rectilinear grid's axes cut, explicit points taken. VTK's vtkExtractGrid.

An axis of one point is dropped: a k-plane of a 3D grid is a 2D grid, whose cells are faces of the 3D grid's. Point fields come along at the points kept; each cell takes the cell fields, marks and ids of the input cell at its first corner -- the cell beyond a plane, or before it at the grid's far side -- as VTK takes them. The output's offset is its first point's index, so its extent is the box's (with rate 1).

Source code in packages/tack-vis/src/tack/dataset/algorithms/extract_extent.py 
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
@on_inputs_device
def extract_extent(data, extent, rate=1, fields=None):
    """The points of ``data``, a structured grid, in the box ``extent`` -- one
    (low, high) pair of point indices per axis, both included, or VTK's flat
    (i0, i1, j0, j1, k0, k1), in the whole grid's indices (as VTK's extents
    count them, from the topology's ``offset``), clipped to the grid -- every
    ``rate``-th along each axis (a number, or one per axis): a structured grid
    of the same kind, an image's coordinates still uniform, a rectilinear
    grid's axes cut, explicit points taken. VTK's ``vtkExtractGrid``.

    An axis of one point is dropped: a k-plane of a 3D grid is a 2D grid,
    whose cells are faces of the 3D grid's. Point fields come along at the
    points kept; each cell takes the cell fields, marks and ids of the input
    cell at its first corner -- the cell beyond a plane, or before it at the
    grid's far side -- as VTK takes them. The output's ``offset`` is its first
    point's index, so its extent is the box's (with ``rate`` 1)."""
    topology = data.topology
    if not isinstance(topology, StructuredTopology):
        raise TypeError("extract_extent takes a structured grid (an image, a rectilinear "
                        "or a curvilinear grid); extract_cells selects any mesh's cells")
    dims = list(topology.point_dims)
    n = len(dims)
    rates = [int(rate)] * n if np.ndim(rate) == 0 else [int(r) for r in rate]
    if len(rates) != n or min(rates) < 1:
        raise ValueError(f"rate is a step of at least 1, or one per axis, not {rate!r}")
    offset = list(topology.offset)
    starts, counts = [], []
    for (lo, hi), d, o, r in zip(_ranges(extent, n), dims, offset, rates, strict=True):
        lo, hi = max(lo, o), min(hi, o + d - 1)
        if hi < lo:
            raise ValueError(f"the extent {extent!r} misses the grid's points "
                             f"{[(o, o + d - 1) for o, d in zip(offset, dims)]}")
        starts.append(lo - o)
        counts.append((hi - lo) // r + 1)
    kept = [a for a in range(n) if counts[a] > 1]
    if not kept:
        raise ValueError(f"the extent {extent!r} holds one point of the grid: no cells")
    out_dims = [counts[a] for a in kept]
    device = data.device
    out = StructuredTopology(out_dims, id_dtype=topology.id_dtype, device=device,
                             offset=[offset[a] + starts[a] for a in kept])
    # Each output cell's first corner, in the input's cells; along a dropped
    # axis, the cell beyond the point, or before it at the far side.
    cells = [d - 1 for d in dims]
    cell_starts = [min(s, c - 1) if counts[a] == 1 else s
                   for a, (s, c) in enumerate(zip(starts, cells))]
    cell_counts = [max(c - 1, 1) for c in counts]
    point_ids = _ids(dims, starts, counts, rates, topology.id_dtype)
    cell_ids = (_ids(cells, cell_starts, cell_counts, rates, topology.id_dtype)
                if out_dims and min(cells) > 0 else tack.empty((0,), topology.id_dtype))
    coordinates = data.geometry.values
    geometry = None
    if isinstance(data.geometry.space, H1) and isinstance(coordinates,
                                                          (UniformCoordinates, CartesianProduct)):
        geometry = Field(H1(out), _strided_geometry(coordinates, starts, counts, rates, device),
                         data.unit)
    return carry(data, out, points=Take(point_ids), cells=Take(cell_ids), geometry=geometry,
                 fields=fields)

mask 

mask(data, stride, fields=None)

Every stride-th cell (0, stride, 2 * stride, ...), as Viskores' Mask, with only the points they use: fields carry as threshold carries them. mask_points is its point twin.

Source code in packages/tack-vis/src/tack/dataset/algorithms/mask.py 
10
11
12
13
14
15
16
17
18
19
20
@on_inputs_device
def mask(data, stride, fields=None):
    """Every ``stride``-th cell (0, stride, 2 * stride, ...), as Viskores' ``Mask``,
    with only the points they use: fields carry as ``threshold`` carries them.
    ``mask_points`` is its point twin."""
    if int(stride) < 1:
        raise ValueError(f"stride must be at least 1, not {stride!r}")
    keep = tack.empty((data.num_cells,), tack.i32)
    if data.num_cells:
        _every(keep, int(stride))
    return _keep_cells(data, keep, fields)

Cutting

New datasets of pieces of cells, with points made on their edges.

Algorithm Gives Takes VTK Viskores
contour isosurfaces, isolines or points of a field, one level or several dataset vtkContourFilter Contour
slice where an implicit function is zero (or a level) dataset vtkCutter Slice
slice_plane a plane's cut, or parallel planes dataset vtkCutter Slice
clip the part above (or below) a field's level or inside an implicit function dataset vtkTableBasedClipDataSet ClipWithField, ClipWithImplicitFunction

contour 

contour(
    data,
    field,
    value,
    merge_points=True,
    fields=None,
    dimension=None,
)

Where field crosses value: a surface of a 3D mesh, as triangles (as polygons on polyhedra); lines of a 2D mesh.

value may be a sequence: the surfaces of every value, in one dataset made in one pass over the cells (each cell compared with every value), with each output cell's value in a cell field contour_value. Points are not shared between values.

field is a scalar field with a basis, or its name: H1 of order 1 (point data), or L2 (DG, an order per cell too). An order-2 H1 field is refused: its cells would need subdividing, and their corner values alone are another, linear field. Each 3D cell contributes the triangles of VTK's case tables for its corner values, interpolated along its edges; in a mesh with 3D cells, cells of lower dimension contribute nothing, and a cell constant has no level set inside its cells.

dimension (3, 2 or 1) chooses which cells to contour, by default the mesh's highest: the cells of that dimension alone, as a mesh of their own on the points they use (no cells of it, no contour), so a mesh of tetrahedra and the triangles of its boundary gives the tetrahedra's surface, or with dimension=2 the triangles' lines. Cells of other dimensions are never mixed in.

A 2D mesh -- triangles, quads, pixels or polygons -- gives line cells from its point values: each cell's crossings, found going round its edges, are joined in pairs, each line with the higher values on its right as the cell winds. Lines share their points through the mesh's faces (its edges), whatever merge_points says -- unless the field or the geometry is L2: then each cell's lines have their own points, from its own values and corners. No cells give no lines.

A 1D mesh -- lines or polylines -- gives a vertex cell at each crossing of its segments, from its point values, interpolated from the lower point id; a crossing at a point whose value is the level is that point, and with merge_points every crossing is one vertex, each taking the cell fields of the first cell to cross there.

With merge_points, triangles share their points: by global edge when the field and the geometry are both continuous, so the surface is watertight; otherwise only within each cell, so a DG field's surface keeps its jumps and an L2 geometry's its gaps. H1 fields carry onto the surface's points, interpolated along the same edges; cell fields (Constant, values on cells) onto its triangles, from the cell each came from.

Source code in packages/tack-vis/src/tack/dataset/algorithms/contour.py 
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
@on_inputs_device
def contour(data, field, value, merge_points=True, fields=None, dimension=None):
    """Where ``field`` crosses ``value``: a surface of a 3D mesh, as triangles
    (as polygons on polyhedra); lines of a 2D mesh.

    ``value`` may be a sequence: the surfaces of every value, in one dataset
    made in one pass over the cells (each cell compared with every value),
    with each output cell's value in a cell field ``contour_value``. Points
    are not shared between values.

    ``field`` is a scalar field with a basis, or its name: ``H1`` of order 1
    (point data), or ``L2`` (DG, an order per cell too). An order-2 ``H1``
    field is refused: its cells would need subdividing, and their corner
    values alone are another, linear field. Each 3D cell contributes the
    triangles of VTK's case tables for its corner values, interpolated along
    its edges; in a mesh with 3D cells, cells of lower dimension contribute
    nothing, and a cell constant has no level set inside its cells.

    ``dimension`` (3, 2 or 1) chooses which cells to contour, by default the
    mesh's highest: the cells of that dimension alone, as a mesh of their own on
    the points they use (no cells of it, no contour), so a mesh of tetrahedra
    and the triangles of its boundary gives the tetrahedra's surface, or with
    ``dimension=2`` the triangles' lines. Cells of other dimensions are never
    mixed in.

    A 2D mesh -- triangles, quads, pixels or polygons -- gives line cells from
    its point values: each cell's crossings, found going round its edges, are
    joined in pairs, each line with the higher values on its right as the
    cell winds. Lines share their points through the mesh's faces (its edges),
    whatever ``merge_points`` says -- unless the field or the geometry is
    ``L2``: then each cell's lines have their own points, from its own values
    and corners. No cells give no lines.

    A 1D mesh -- lines or polylines -- gives a vertex cell at each crossing of
    its segments, from its point values, interpolated from the lower point id;
    a crossing at a point whose value is the level is that point, and with
    ``merge_points`` every crossing is one vertex, each taking the cell fields
    of the first cell to cross there.

    With ``merge_points``, triangles share their points: by global edge when
    the field and the geometry are both continuous, so the surface is
    watertight; otherwise only within each cell, so a DG field's surface
    keeps its jumps and an ``L2`` geometry's its gaps. ``H1`` fields carry
    onto the surface's points, interpolated along the same edges; cell
    fields (``Constant``, values on cells) onto its triangles, from the cell
    each came from.
    """
    value = _in_unit_of(_field(data, field), value, "contour's value")
    from tack.dataset.ghosts import visible

    data, field = visible(data, field)
    field = _field(data, field)
    space = field.space
    if isinstance(space, H1) and space.order != 1:
        raise NotImplementedError(
            f"contour of an order-{space.order} H1 field is not built: it would subdivide "
            "each cell, and its corner values alone give the contour of a linear field")
    _scalar_floating("contour", field)
    levels, nlevels, several = _levels(value, field)
    top = data.topology.dimension
    if dimension is None:
        dimension = top if top >= 1 else 3       # no cells (or vertices): no surface
    if dimension not in (1, 2, 3):
        raise ValueError(f"contour's dimension is 3, 2 or 1 (the cells to contour), "
                         f"not {dimension!r}")
    mixed = any(g.count and g.shape.DIMENSION != dimension for g in data.topology.groups())
    if mixed or dimension > top:
        data, field = _cells_of_dimension(data, field, dimension)
    if dimension == 2:
        return _contour_lines(data, field, levels, nlevels, several, fields)
    if dimension == 1:
        return _contour_vertices(data, field, levels, nlevels, several, merge_points, fields)
    if not getattr(data.topology, "reference_cells", True):
        return _polyhedral_contour(data, field, levels, nlevels, several, fields)
    field = _with_basis(data, field)
    space = field.space
    if not space.interpolated or space.on not in ("cells", "points"):
        raise TypeError(f"contour needs a field with a basis on the cells, not {space!r}")
    per_cell = 0 if isinstance(space, H1) and isinstance(data.geometry.space, H1) else 1
    groups = [g for g in data.launch_groups("cells", [field])
              if g.count and g.shape.CONTOUR_TRIANGLES]
    n_cells = data.num_cells
    idt = data.id_dtype
    odt = ids.for_output(data, 15 * n_cells * nlevels)     # up to five triangles a cell
    # Keys hold the level beside a point or cell id.
    kdt = ids.at_least(idt, (n_cells if per_cell else data.num_points) * nlevels)
    counts = tack.zeros((n_cells,), tack.i32)
    for group in groups:
        _contour_counts(data.domain_view("cells", group), field.view(group), levels, nlevels,
                        counts)
    starts, triangles = _scanned(counts, n_cells, odt)

    n = 3 * triangles
    points = tack.empty((n,), tack.vec(3, data.dtype))
    edges = tack.empty((n,), tack.vec(2, idt))
    weights = tack.empty((n,), arrays.dtype_of(field.values))
    sources = tack.empty((n,), idt)
    keys = tack.empty((n,), tack.vec(2, kdt))
    which = tack.empty((n,), tack.i32)
    for group in groups:
        _contour_points(data.domain_view("cells", group), field.view(group), levels, nlevels,
                        per_cell, starts, points, edges, weights, sources, keys, which)
    if triangles:
        keep = tack.empty((triangles,), tack.i32)
        _distinct_triangles(keys, keep)
        slots = tack.empty((triangles,), odt)
        kept = exclusive_scan(keep, slots, triangles)
        if kept < triangles:
            arrays_in = (points, edges, weights, sources, keys, which)
            arrays_out = tuple(_like(a, 3 * kept) for a in arrays_in)
            _kept_triangles(keep, slots, *arrays_in, *arrays_out)
            points, edges, weights, sources, keys, which = arrays_out
            triangles, n = kept, 3 * kept
    triangle_cells = tack.empty((triangles,), idt)
    triangle_levels = tack.empty((triangles,), tack.i32)
    if triangles:
        _triangle_cells(sources, which, triangle_cells, triangle_levels)

    if merge_points and n:
        # A key is the edge's (low, high) points, or (within cells) (cell, local edge),
        # its first part times the levels plus the level.
        order = key_order(keys, (n_cells if per_cell else data.num_points) * nlevels,
                                odt)
        offsets, count = run_offsets(keys, order)
        merged = tack.empty((count,), tack.vec(3, data.dtype))
        merged_edges = tack.empty((count,), tack.vec(2, idt))
        merged_weights = tack.empty((count,), weights.dtype)
        _first_of_runs(order, offsets, points, edges, weights, merged, merged_edges,
                       merged_weights, count)
        rows = tack.empty((n,), odt)
        _run_of_each(order, offsets, rows, count)
        points, edges, weights = merged, merged_edges, merged_weights
    else:
        rows, count = CountingArray(n, dtype=odt), n          # each triangle its own points

    # Every cell a triangle, three points each: the types and offsets are
    # implicit, and the connectivity is the triangles' rows of merged points.
    types, tri_offsets = uniform_cells(shapes.TRIANGLE, triangles, odt)
    surface = UnstructuredTopology(types, tri_offsets, rows, num_points=count)
    return _with_levels(carry(data, surface, points=Interpolate(edges, weights),
                              cells=Pieces(triangle_cells), geometry=Field(H1(surface), points),
                              fields=fields), levels, triangle_levels, several,
                        field.unit)

slice 

slice(
    data,
    function,
    merge_points=True,
    fields=None,
    value=0.0,
)

Where function (tack.dataset.implicit) is value (zero by default): a contour of its values on the geometry, carrying fields as contour does (merge_points and fields are contour's). value may be a sequence, as contour's: every level in one pass, each output cell's in contour_value. slice_plane is this with a Plane. On a curved geometry (H1 above order 1) the cells are cut linearized at their corners, as contour would need its cells subdivided for the curve inside them.

Source code in packages/tack-vis/src/tack/dataset/algorithms/slice.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
@on_inputs_device
def slice(data, function, merge_points=True, fields=None, value=0.0):
    """Where ``function`` (``tack.dataset.implicit``) is ``value`` (zero by default):
    a contour of its values on the geometry, carrying fields as ``contour``
    does (``merge_points`` and ``fields`` are ``contour``'s). ``value`` may be
    a sequence, as ``contour``'s: every level in one
    pass, each output cell's in ``contour_value``. ``slice_plane`` is this
    with a ``Plane``. On a curved geometry (``H1`` above order 1) the cells are cut
    linearized at their corners, as ``contour`` would need its cells subdivided
    for the curve inside them."""
    values = implicit_values(data, function)
    if isinstance(values.space, H1) and values.space.order != 1:
        # The function's values at the points: projected into order-1 H1, a
        # continuous field keeps them. contour refuses order 2, so slicing a
        # curved dataset raised.
        values = project(data, values, H1(data))
    return contour(data, values, value, merge_points=merge_points, fields=fields)

slice_plane 

slice_plane(
    data,
    origin,
    normal,
    merge_points=True,
    fields=None,
    value=0.0,
)

The cut of the 3D cells by the plane through origin with normal (any length), as triangles: contour at zero of each position's signed distance to the plane along normal. The distance is computed in the geometry's own space -- per point, or per cell corner for an L2 geometry -- so the cut is the cells' own. value, a signed distance along the unit normal, or a sequence of them, gives parallel planes, all in one pass. slice with a Plane; merge_points and fields are contour's.

Source code in packages/tack-vis/src/tack/dataset/algorithms/slice.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
@on_inputs_device
def slice_plane(data, origin, normal, merge_points=True, fields=None, value=0.0):
    """The cut of the 3D cells by the plane through ``origin`` with ``normal`` (any
    length), as triangles: ``contour`` at zero of each position's signed distance
    to the plane along ``normal``. The distance is computed in the geometry's own
    space -- per point, or per cell corner for an ``L2`` geometry -- so the cut is
    the cells' own. ``value``, a signed distance along the unit normal, or a
    sequence of them, gives parallel planes, all in one pass. ``slice`` with a
    ``Plane``; ``merge_points`` and ``fields`` are ``contour``'s."""
    origin = _length(data, origin, "slice_plane's origin")
    value = _length(data, value, "slice_plane's distance")
    return slice(data, Plane(origin, normal), merge_points=merge_points, fields=fields,
                 value=value)

clip 

clip(data, by, value=0.0, inside=False, fields=None)

The part of data where by -- point data (a field or its name) or an implicit function (tack.dataset.implicit) -- is at or above value, or with inside below it: Viskores' ClipWithField and (value 0) ClipWithImplicitFunction, by its case tables. Cells wholly kept keep their shape; cut cells become the pieces their case lists. Points are the kept corners, the crossings on cut edges (merged across cells) and centroid points; point data is interpolated onto them, cell data goes to each piece.

A polyhedral topology is clipped cell by cell into polyhedra: each cut cell becomes the polyhedron of its kept part, its faces clipped (each face once, so faces two cells share stay one face) and its iso-polygons -- as contour traces them, by López's method -- as caps wound out of it; the output is a PolyhedralTopology. A point at exactly value is kept on either side (inside keeps at or below). Polygons and polylines are not clipped yet.

Source code in packages/tack-vis/src/tack/dataset/algorithms/clip.py 
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
@on_inputs_device
def clip(data, by, value=0.0, inside=False, fields=None):
    """The part of ``data`` where ``by`` -- point data (a field or its name) or an
    implicit function (``tack.dataset.implicit``) -- is at or above ``value``, or
    with ``inside`` below it: Viskores' ClipWithField and (``value`` 0)
    ClipWithImplicitFunction, by its case tables. Cells wholly kept keep their
    shape; cut cells become the pieces their case lists. Points are the kept
    corners, the crossings on cut edges (merged across cells) and centroid
    points; point data is interpolated onto them, cell data goes to each piece.

    A polyhedral topology is clipped cell by cell into polyhedra: each cut cell
    becomes the polyhedron of its kept part, its faces clipped (each face once,
    so faces two cells share stay one face) and its iso-polygons -- as
    ``contour`` traces them, by López's method -- as caps wound out of it; the
    output is a ``PolyhedralTopology``. A point at exactly ``value`` is kept on
    either side (``inside`` keeps at or below). Polygons and polylines are not
    clipped yet."""
    from tack.dataset.ghosts import visible

    data, by = visible(data, by)
    from tack.dataset.polyhedra import PolygonalTopology, PolyhedralTopology

    polyhedral = isinstance(data.topology, PolyhedralTopology) and \
        not isinstance(data.topology, PolygonalTopology)
    if not polyhedral and not getattr(data.topology, "reference_cells", True):
        raise NotImplementedError("clip takes cells of the linear shapes or polyhedra, not "
                                  "polygons or polylines")
    if not isinstance(data.geometry.space, H1):
        # Clip's points are merged across cells, so each needs one position.
        raise TypeError("clip needs an H1 geometry: one position per point")
    if hasattr(by, "value") and not isinstance(by, (str, Field)):
        field = implicit_values(data, by)
    else:
        field = _field(data, by)
        space = field.space
        if not (isinstance(space, H1) or (isinstance(space, Values) and space.on == "points")):
            raise TypeError(f"clip takes point data or an implicit function, not {space!r}")
    # In the field's unit; an implicit function's values have none.
    value = _in_unit_of(field, value, "clip's value")
    if arrays.width_of(field.values):
        raise TypeError("clip needs a scalar field")
    precision = _floating_dtype(field.values)
    values = _converted(_point_values(field.values, data.num_points), precision)
    if polyhedral:
        from tack.dataset.algorithms._clip_polyhedra import clip_polyhedra

        return clip_polyhedra(data, arrays.materialize(values), float(value), inside, fields)
    invert = 1 if inside else 0           # the tables' name: keep what is below
    n = data.num_cells

    groups = []
    for group in data.topology.groups():
        if group.count:
            t = _shape_tables(int(group.shape.ID), bool(invert))
            if t is None:
                raise NotImplementedError(f"clip has no table for {group.shape.__name__} cells")
            groups.append((group, t))
    counts = tack.empty((n,), tack.vec(_KINDS, tack.i32))
    if n:
        counts.fill(0)
    for group, (table, cases, _, remap, _, _) in groups:
        _clip_counts(group.view(), values, float(value), invert, table, cases, remap, counts)
    idt = data.id_dtype
    odt = ids.for_output(data, _MOST * n)
    starts = tack.empty((n,), tack.vec(_KINDS, odt))
    if n <= _SERIAL_KINDS:
        # A small mesh: every kind's scan in one thread, one launch and one read,
        # where the parallel scans below are eleven launches.
        sums = tack.empty((1,), tack.vec(_KINDS, odt))
        if n:
            _scan_kinds(counts, starts, sums, n)
        totals = [int(v) for v in sums.numpy().reshape(-1)] if n else [0] * _KINDS
    else:
        totals, parts = [], []
        for kind in range(_KINDS):
            part = tack.empty((n,), tack.i32)
            kind_starts = tack.empty((n,), odt)
            _split_counts(counts, kind, part)
            totals.append(exclusive_scan(part, kind_starts, n))
            parts.append(kind_starts)
        _join_starts(*parts, starts)
    n_cells, n_conn, n_records, n_centroids, n_members = totals

    types = tack.empty((n_cells,), tack.u8)
    offsets = tack.zeros((n_cells + 1,), odt)
    connectivity = tack.empty((n_conn,), odt)
    sources = tack.empty((n_cells,), idt)
    keys = tack.empty((n_records,), tack.vec(2, idt))
    slots = tack.empty((n_records,), odt)
    centroid_offsets = tack.empty((n_centroids + 1,), odt)
    centroid_cells = tack.empty((n_centroids,), idt)
    for group, (table, cases, edges, remap, whole_type, whole_order) in groups:
        _clip_emit(group.view(), values, float(value), invert, table, cases, edges, remap,
                   whole_type, whole_order, starts, types, offsets, connectivity, sources,
                   keys, slots, centroid_offsets, centroid_cells)
    _close(offsets, n_cells, n_conn)
    _close(centroid_offsets, n_centroids, n_members)

    order = key_order(keys, data.num_points, odt)
    point_runs, count = run_offsets(keys, order)
    members = tack.empty((n_members,), odt)
    ends = tack.empty((count,), tack.vec(2, idt))
    # Weights are fractions even when the values are integers (Viskores interpolates
    # an integer field's crossings the same way).
    weights = tack.empty((count,), precision)
    if count:
        _resolve_points(keys, order, point_runs, slots, values, float(value), connectivity,
                        members, ends, weights, count)
    if n_conn:
        _centroid_entries(connectivity, count)
    out = UnstructuredTopology(types, offsets, connectivity, num_points=count + n_centroids)
    points = Interpolate(ends, weights, averages=(centroid_offsets, members),
                         average_cells=centroid_cells)
    return carry(data, out, points=points, cells=Pieces(sources), fields=fields)

Surfaces and edges

A mesh's faces and edges, as sets or datasets on the same points.

Algorithm Gives Takes VTK Viskores
boundary_faces the boundary, as a face set dataset — —
extract_surface faces (the boundary by default) as cells dataset vtkDataSetSurfaceFilter ExternalFaces
extract_edges edges as line cells dataset vtkExtractEdges —
feature_edges a surface's boundary, feature, non-manifold and manifold edges dataset vtkFeatureEdges —
surface_normals cell and point normals of a surface, wound consistently dataset vtkPolyDataNormals SurfaceNormals

boundary_faces 

boundary_faces(data)

The faces with one side -- the boundary -- as a face Set: a 3D mesh's boundary faces, a 2D mesh's boundary edges, a polyline mesh's free ends, numbered as data.topology.faces() numbers them. Name it in data.sets (data.sets["boundary"] = boundary_faces(data)) to iterate over it or carry it through filters; extract_surface makes a dataset of it.

Source code in packages/tack-vis/src/tack/dataset/algorithms/boundary_faces.py 
 9
10
11
12
13
14
15
16
@on_inputs_device
def boundary_faces(data):
    """The faces with one side -- the boundary -- as a face ``Set``: a 3D mesh's
    boundary faces, a 2D mesh's boundary edges, a polyline mesh's free ends,
    numbered as ``data.topology.faces()`` numbers them. Name it in ``data.sets``
    (``data.sets["boundary"] = boundary_faces(data)``) to iterate over it or
    carry it through filters; ``extract_surface`` makes a dataset of it."""
    return Set._from_sorted_ids(data.topology, "faces", data.topology.faces().boundary())

extract_surface 

extract_surface(data, faces=None, fields=None)

The faces -- a face Set (as boundary_faces returns), a set's name in data.sets, or face ids (a field or an array); by default the boundary (boundary_faces), VTK's external faces -- as a dataset on the same points: a 3D mesh's faces as triangles, quads or polygons, each in its side 0's outward order; a 2D mesh's edges as line cells; a polyline mesh's ends as vertex cells. By default a 2D mesh is its own surface: its cells pass through (its hidden ones left out), as VTK's geometry filter passes polygons; its boundary edges are extract_surface(data, boundary_faces(data)). Fields on faces become the surface's cell fields, and so do cell fields (Constant, values on cells), each face taking its side 0 cell's value; point fields (H1, values on points) stay on the same points. A 3D mesh's order-2 H1 fields and geometry stay order 2 on the surface, the trace of the solid's: each face's corners, edge middles and (a quad's) center are the surface cell's, so their values are gathered, not evaluated. A 2D mesh's boundary lines take an order-2 field's values at the points. The geometry must be H1: an L2 one has no shared points for the faces to stand on.

Source code in packages/tack-vis/src/tack/dataset/algorithms/extract_surface.py 
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
@on_inputs_device
def extract_surface(data, faces=None, fields=None):
    """The ``faces`` -- a face ``Set`` (as ``boundary_faces`` returns), a set's name
    in ``data.sets``, or face ids (a field or an array); by default the boundary
    (``boundary_faces``), VTK's external faces -- as a dataset on the same points: a 3D
    mesh's faces as triangles, quads or polygons, each in its side 0's outward
    order; a 2D mesh's edges as line cells; a polyline mesh's ends as vertex
    cells. By default a 2D mesh is its own surface: its cells pass through (its
    hidden ones left out), as VTK's geometry filter passes polygons; its
    boundary edges are ``extract_surface(data, boundary_faces(data))``. Fields on faces become the surface's
    cell fields, and so do cell fields (``Constant``, values on cells), each face
    taking its side 0 cell's value; point fields (``H1``, values on points) stay on
    the same points. A 3D mesh's order-2 ``H1`` fields and geometry stay order 2
    on the surface, the trace of the solid's: each face's corners, edge middles
    and (a quad's) center are the surface cell's, so their values are gathered,
    not evaluated. A 2D mesh's boundary lines take an order-2 field's values at
    the points. The geometry must be ``H1``: an ``L2`` one has no shared points
    for the faces to stand on."""
    from tack.dataset.ghosts import refuse_hidden, visible

    if faces is None:
        data = visible(data)
    else:
        refuse_hidden(data, "extract_surface of chosen faces")
    if not isinstance(data.geometry.space, H1):
        raise TypeError("extract_surface keeps the points, so needs an H1 geometry")
    if faces is None and data.topology.dimension == 2:
        # A surface mesh is its own surface: its polygons pass through, as VTK's
        # geometry filter passes them.
        return carry(data, data.topology, points=Same(), cells=Same(), fields=fields)
    if faces is None:
        fresh = _fresh_boundary(data, fields)
        if fresh is not None:
            return fresh
        from tack.dataset.algorithms.boundary_faces import boundary_faces
        faces = boundary_faces(data)
    ids = _face_ids(data, faces)
    if not getattr(data.topology, "reference_cells", True):
        return _polygon_surface(data, ids, fields)
    faces = data.topology.faces()
    idt = data.id_dtype
    n = ids.shape[0]
    # Up to 4 points a face: past the topology's own bound (12 a cell) for a
    # surface of faces that share no points.
    odt = id_types.for_output(data, 4 * n)
    sizes = tack.empty((n,), tack.i32)
    starts = tack.empty((n,), odt)
    if n:
        _surface_sizes(ids, faces.kinds, sizes)
    length = exclusive_scan(sizes, starts, n) if n else 0
    types = tack.empty((n,), tack.u8)
    offsets = tack.zeros((n + 1,), odt)
    connectivity = tack.empty((length,), odt)
    if n:
        _surface_cells(ids, faces.kinds, faces.rows, starts, types, offsets, connectivity,
                       length, n)
    surface = UnstructuredTopology(types, offsets, connectivity, num_points=data.num_points)
    cells = tack.empty((n,), idt)
    if n:
        _face_cells(ids, faces.sides, cells)
    out = carry(data, surface, points=Same(), cells=Pieces(cells), faces_to_cells=Take(ids),
                fields=fields)
    if data.topology.dimension == 3:
        _keep_order_two(data, out, ids, fields)
    return out

extract_edges 

extract_edges(data, edges=None, fields=None)

The edges -- an edge Set, a set's name in data.sets, or edge ids (a field or an array); by default every edge of the topology (edges()) -- as a dataset of line cells on the same points, cell i running from edge i's lower point id to its higher (edges().rows). Values on edges (Values(data, "edges")) become the output's cell values, and edge sets its cell sets; point fields (H1, values on points) and point sets stay on the same points -- an order-2 field's or geometry's values at the points, the lines being linear. Cell and face fields are left behind: an edge lies on every cell and face around it. The geometry must be H1: an L2 one has no shared points for the lines to stand on.

Source code in packages/tack-vis/src/tack/dataset/algorithms/extract_edges.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
@on_inputs_device
def extract_edges(data, edges=None, fields=None):
    """The ``edges`` -- an edge ``Set``, a set's name in ``data.sets``, or edge ids (a
    field or an array); by default every edge of the topology (``edges()``) --
    as a dataset of line cells on the same points, cell ``i`` running from edge
    ``i``'s lower point id to its higher (``edges().rows``). Values on edges
    (``Values(data, "edges")``) become the output's cell values, and edge sets
    its cell sets; point fields (``H1``, values on points) and point sets stay
    on the same points -- an order-2 field's or geometry's values at the
    points, the lines being linear. Cell and face fields are left behind: an
    edge lies on every cell and face around it. The geometry must be ``H1``:
    an ``L2`` one has no shared points for the lines to stand on."""
    from tack.dataset.ghosts import refuse_hidden, visible

    if edges is None:
        data = visible(data)
    else:
        refuse_hidden(data, "extract_edges of chosen edges")
    if not isinstance(data.geometry.space, H1):
        raise TypeError("extract_edges keeps the points, so needs an H1 geometry")
    topology = data.topology
    found = topology.edges()
    idt = data.id_dtype
    if edges is None:
        n = found.num_edges
        chosen = CountingArray(n, dtype=idt)              # every edge, in order
        edges_to_cells = Same()
    else:
        chosen = _edge_ids(data, edges, found.num_edges)
        n = chosen.shape[0]
        edges_to_cells = Take(chosen)
    types, offsets = uniform_cells(shapes.LINE, n, idt)       # all lines: implicit
    connectivity = tack.empty((2 * n,), idt)
    if n:
        _edge_cells(chosen, found.rows, connectivity, n)
    lines = UnstructuredTopology(types, offsets, connectivity, num_points=data.num_points)
    return carry(data, lines, points=Same(), edges_to_cells=edges_to_cells, fields=fields)

feature_edges 

feature_edges(
    data,
    angle=30.0,
    boundary=True,
    feature=True,
    non_manifold=True,
    manifold=False,
    fields=None,
)

A surface's edges of the kinds asked for, as line cells on the same points, each from its lower point id to its higher, in order of those ids, with kind (EdgeKind) on the cells -- as VTK's vtkFeatureEdges:

  • boundary: edges of one cell;
  • feature: edges of two cells whose normals are at least angle degrees apart (two cells wound opposite ways meet at almost 180);
  • non_manifold: edges of three or more cells;
  • manifold: the other edges of two cells -- all of them when feature edges are not asked for, as in VTK (off by default).

The surface is the 2D cells -- triangles, quads, pixels, polygons; a polygon's normal is Newell's. Other cells are left out; a volume mesh's surface comes from extract_surface first. Point fields stay on the same points (fields chooses which, as extract_edges); cell fields are left behind, an edge belonging to every cell around it. Hidden cells are left out. The geometry must be H1.

Source code in packages/tack-vis/src/tack/dataset/algorithms/feature_edges.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
@on_inputs_device
def feature_edges(data, angle=30.0, boundary=True, feature=True, non_manifold=True,
                  manifold=False, fields=None):
    """A surface's edges of the kinds asked for, as line cells on the same points,
    each from its lower point id to its higher, in order of those ids, with
    ``kind`` (``EdgeKind``) on the cells -- as VTK's ``vtkFeatureEdges``:

    - ``boundary``: edges of one cell;
    - ``feature``: edges of two cells whose normals are at least ``angle``
      degrees apart (two cells wound opposite ways meet at almost 180);
    - ``non_manifold``: edges of three or more cells;
    - ``manifold``: the other edges of two cells -- all of them when feature
      edges are not asked for, as in VTK (off by default).

    The surface is the 2D cells -- triangles, quads, pixels, polygons; a
    polygon's normal is Newell's. Other cells are left out; a volume mesh's
    surface comes from ``extract_surface`` first. Point fields stay on the same
    points (``fields`` chooses which, as ``extract_edges``); cell fields are
    left behind, an edge belonging to every cell around it. Hidden cells are
    left out. The geometry must be ``H1``."""
    from tack.dataset.ghosts import visible
    from tack.dataset.polyhedra import PolygonalTopology

    data = visible(data)
    if not isinstance(data.geometry.space, H1):
        raise TypeError("feature_edges keeps the points, so needs an H1 geometry")
    topology = data.topology
    if topology.dimension == 3:
        raise ValueError("feature_edges takes a surface; take a volume mesh's with "
                         "extract_surface first")
    if not 0.0 <= angle <= 180.0:
        raise ValueError(f"the feature angle is in degrees, 0 to 180, not {angle!r}")
    idt = data.id_dtype
    n_cells = data.num_cells
    # Each 2D cell's loop of point ids.
    if isinstance(topology, PolygonalTopology):
        offsets, loops = (ids.as_ids(a, idt) for a in topology.loops())
    else:
        sizes = tack.zeros((n_cells,), tack.i64)
        groups = [g for g in data.launch_groups("cells", []) if g.count]
        for group in groups:
            if group.shape in _ROUND:
                _loop_sizes(data.domain_view("cells", group), len(_ROUND[group.shape]), sizes)
        offsets = tack.empty((n_cells + 1,), tack.i64)
        total = int(exclusive_scan(sizes, offsets, n_cells, close=True))
        loops = tack.empty((total,), idt)
        for group in groups:
            if group.shape in _ROUND:
                order = _table(("round", group.shape), _ROUND[group.shape])
                _shape_loops(data.domain_view("cells", group), order, offsets, loops)
    total = int(loops.shape[0])
    positions = materialize(_point_values(data.geometry.values, data.num_points))
    keys = tack.empty((total,), tack.vec(2, idt))
    owners = tack.empty((total,), tack.i64)
    normals = _vectors(n_cells, data.dtype)
    if n_cells:
        _sides(loops, offsets, positions, n_cells, keys, owners, normals)
    count = 0
    if total:
        order = key_order(keys, data.num_points, idt)
        runs, count = run_offsets(keys, order)
    wanted = tack.array(np.array([0, boundary, feature, non_manifold, manifold], np.int32))
    kinds = tack.empty((count,), tack.i32)
    flags = tack.empty((count,), tack.i32)
    rows = tack.empty((count,), tack.vec(2, idt))
    chosen = 0
    slots = tack.empty((count,), tack.i64)
    if count:
        _classify(order, runs, keys, owners, normals, math.cos(math.radians(angle)), wanted,
                  kinds, flags, rows, count)
        chosen = int(exclusive_scan(flags, slots, count))
    connectivity = tack.empty((2 * chosen,), idt)
    kind_out = tack.empty((chosen,), tack.i32)
    if chosen:
        _chosen(flags, slots, kinds, rows, connectivity, kind_out, count)
    types, line_offsets = uniform_cells(shapes.LINE, chosen, idt)
    lines = UnstructuredTopology(types, line_offsets, connectivity, num_points=data.num_points)
    out = carry(data, lines, points=Same(), fields=fields)
    out.fields["kind"] = Field(Values(out, "cells"), kind_out)
    return out

EdgeKind

Bases: IntEnum

The values of feature_edges' kind cell field.

surface_normals 

surface_normals(data, consistent=True)

A surface -- triangles, quads, pixels or polygons -- with its normals: the unit normal of each cell (Newell's, round its corners) as the cell field cell_normals, and at each point the cells' around it summed and made unit as the point field normals (zero at a point of no cells), as VTK's vtkPolyDataNormals computes them without splitting.

With consistent, the polygons are first wound alike: each connected region's as its first (lowest-numbered) cell winds, by union-find over the edges on the device, the cells that wind against it turned round (their corners reversed from the first), so their normals agree -- VTK's consistency, which is on by default there too. A surface no winding makes consistent (a Moebius strip) is refused, as is an edge of three or more polygons (the mesh's faces refuse it). A polygonal topology is consistent already (it refuses to be made otherwise), so it is left as it is. The result keeps every field; the normals are dimensionless.

Source code in packages/tack-vis/src/tack/dataset/algorithms/surface_normals.py 
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
@on_inputs_device
def surface_normals(data, consistent=True):
    """A surface -- triangles, quads, pixels or polygons -- with its normals: the
    unit normal of each cell (Newell's, round its corners) as the cell field
    ``cell_normals``, and at each point the cells' around it summed and made unit
    as the point field ``normals`` (zero at a point of no cells), as VTK's
    ``vtkPolyDataNormals`` computes them without splitting.

    With ``consistent``, the polygons are first wound alike: each connected
    region's as its first (lowest-numbered) cell winds, by union-find over the
    edges on the device, the cells that wind against it turned round (their
    corners reversed from the first), so their normals agree -- VTK's
    consistency, which is on by default there too. A surface no winding makes
    consistent (a Moebius strip) is refused, as is an edge of three or more
    polygons (the mesh's faces refuse it). A polygonal topology is consistent
    already (it refuses to be made otherwise), so it is left as it is. The
    result keeps every field; the normals are dimensionless."""
    from tack.dataset.ghosts import visible

    data = visible(data)
    if data.num_cells and data.topology.dimension != 2:
        raise ValueError(f"surface_normals takes 2D cells (a surface), not a mesh of "
                         f"dimension {data.topology.dimension}: extract_surface a solid "
                         "first")
    if consistent and isinstance(data.topology, UnstructuredTopology) and data.num_cells:
        data = _consistent(data)
    else:
        data = carry(data, data.topology, points=Same(), cells=Same())    # not the input
    n = data.num_cells
    cell_normals = _vectors(n, data.dtype)
    if not getattr(data.topology, "reference_cells", True):
        data.topology.faces()                      # a polygon's walk
    for group in data.launch_groups("cells", []):
        if not group.count:
            continue
        view = data.domain_view("cells", group)
        if group.shape in _ROUND:
            _shape_normals(view, _table(("round", group.shape), _ROUND[group.shape]),
                           cell_normals)
        else:
            _polygon_normals(view, cell_normals)
    points = tack.empty((data.num_points,), tack.vec(3, data.dtype))
    points.fill(0)
    # Each cell's normal at each of its corners, as point_links lays them out
    # (by cell and corner on shapes, by cell_points' entries elsewhere).
    if getattr(data.topology, "reference_cells", True):
        start_of, total = corner_layout(data.topology)
        corners = _vectors(total, data.dtype)
        for group in data.launch_groups("cells", []):
            if group.count:
                _corner_normals(group.view(), cell_normals, start_of[group.root], corners)
    else:
        point_offsets, point_ids = data.topology.cell_points()
        total = point_ids.shape[0]
        owners = tack.empty((total,), data.id_dtype)
        if total:
            _entry_cells(point_offsets, owners)
        corners = _take(cell_normals, owners)
    if total:
        offsets, links = data.topology.point_links()
        _average_links(offsets, links, corners, points)
        _unit(points)
    data.fields["normals"] = Field(Values(data, "points"), points, "dimensionless")
    data.fields["cell_normals"] = Field(Values(data, "cells"), cell_normals, "dimensionless")
    return data

Connectivity and winding

Which cells belong together, which are wound wrong, and which points are one.

Algorithm Gives Takes VTK Viskores
connected_components regions of cells connected through points (or within a radius) dataset, pieces vtkConnectivityFilter, vtkEuclideanClusterExtraction CellSetConnectivity
check_winding cells whose faces do not wind consistently, as a cell set dataset — —
merge_points coincident points (or points within a tolerance) made one dataset vtkStaticCleanUnstructuredGrid CleanGrid

connected_components 

connected_components(
    data,
    radius=None,
    min_size=None,
    periodic=None,
    locator=None,
)

Which cells are connected, through shared points: data with region, a field on its cells and one on its points, the regions numbered from 0 in order of their lowest cell id, as VTK's vtkConnectivityFilter (RegionId, all regions) and Viskores' CellSetConnectivity: region on the cells, point_region on the points. A point takes its cells' region; a point no cell uses, and a hidden cell, has -1. The largest region, or any, is a threshold of region away.

Computed on the device: union-find over the cells and their points, each cell-point pair hooking the larger root under the smaller with an atomic minimum, then every node jumping to its root, until nothing changes.

Over a PartitionedDataSet, regions are numbered across every piece of every rank: each piece labels its own, and pieces join the regions that meet at the points they share. With global ids on every piece, regions are numbered in order of their lowest global cell id -- the whole's numbering, for pieces carrying the whole's ids -- else of their lowest (piece, cell). With radius, the pieces are point clouds: each labels itself with copies of the other pieces' points within radius (ghost_points), and a copy joins its region to its owner's. min_size counts each point once, where it is owned.

With radius, points within that distance of each other (inclusive) are connected too: on a point cloud, whose cells are its points, that is Euclidean clustering (VTK's vtkEuclideanClusterExtraction, PCL's) and, with the linking length, friends-of-friends halo finding (yt's FoF). min_size gives regions of fewer points -1, the rest numbered again in the same order. periodic (a box, as PointLocator takes it) joins points across its faces, as a cosmology box needs. locator reuses a PointLocator of data.

Source code in packages/tack-vis/src/tack/dataset/algorithms/connected_components.py 
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
@on_inputs_device
def connected_components(data, radius=None, min_size=None, periodic=None, locator=None):
    """Which cells are connected, through shared points: ``data`` with ``region``, a
    field on its cells and one on its points, the regions numbered from 0 in
    order of their lowest cell id, as VTK's vtkConnectivityFilter (``RegionId``,
    all regions) and Viskores' CellSetConnectivity: ``region`` on the cells,
    ``point_region`` on the points. A point takes its cells' region; a point no
    cell uses, and a hidden cell, has -1. The largest region, or any, is a
    ``threshold`` of ``region`` away.

    Computed on the device: union-find over the cells and their points, each
    cell-point pair hooking the larger root under the smaller with an atomic
    minimum, then every node jumping to its root, until nothing changes.

    Over a ``PartitionedDataSet``, regions are numbered across every piece of
    every rank: each piece labels its own, and pieces join the regions that
    meet at the points they share. With global ids on
    every piece, regions are numbered in order of their lowest global cell id
    -- the whole's numbering, for pieces carrying the whole's ids -- else of
    their lowest (piece, cell). With ``radius``, the pieces are point clouds:
    each labels itself with copies of the other pieces' points within
    ``radius`` (``ghost_points``), and a copy joins its region to its owner's.
    ``min_size`` counts each point once, where it is owned.

    With ``radius``, points within that distance of each other (inclusive) are
    connected too: on a point cloud, whose cells are its points, that is
    Euclidean clustering (VTK's vtkEuclideanClusterExtraction, PCL's) and,
    with the linking length, friends-of-friends halo finding (yt's FoF).
    ``min_size`` gives regions of fewer points -1, the rest numbered again in
    the same order. ``periodic`` (a box, as ``PointLocator`` takes it) joins
    points across its faces, as a cosmology box needs. ``locator`` reuses a
    ``PointLocator`` of ``data``."""
    radius = _length(data, radius, "connected_components's radius")
    periodic = _length(data, periodic, "connected_components's periodic box")
    if radius is not None and not radius >= 0:
        raise ValueError(f"radius must be at least zero, not {radius!r}")
    if isinstance(data, PartitionedDataSet):
        if locator is not None:
            raise ValueError("a PointLocator finds one piece's points; pieces take none")
        if radius is not None:
            return _components_near(data, float(radius), min_size, periodic)
        if periodic is not None:
            raise ValueError("periodic joins points within a radius across the box's faces")
        return _components_pieces(data, min_size)
    cell_region, point_region, count = _labels(data, radius, locator, periodic)
    if min_size is not None and count:
        _smallest_dropped(cell_region, point_region, count, int(min_size))
    return _with_regions(data, cell_region, point_region)

check_winding 

check_winding(data)

The cells whose faces do not wind consistently, as a cell Set.

Walking every face outward (side_point), a closed, consistently wound cell walks each of its edges exactly twice, once each way. A cell that does not -- a face with the wrong side, an open cell, an edge of three faces -- is listed. Combinatorial: geometry plays no part. For polygons: the cells whose edges, walked their way, do not close into a loop. Polylines have no winding, and shape-based cells take their faces from their shape's table, so are refused.

Source code in packages/tack-vis/src/tack/dataset/algorithms/check_winding.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
@on_inputs_device
def check_winding(data):
    """The cells whose faces do not wind consistently, as a cell ``Set``.

    Walking every face outward (``side_point``), a closed, consistently wound cell
    walks each of its edges exactly twice, once each way. A cell that does not --
    a face with the wrong side, an open cell, an edge of three faces -- is listed.
    Combinatorial: geometry plays no part. For polygons: the cells whose edges,
    walked their way, do not close into a loop. Polylines have no winding, and
    shape-based cells take their faces from their shape's table, so are refused.
    """
    t = getattr(data, "topology", data)
    if getattr(t, "reference_cells", True):
        raise TypeError("check_winding checks cells that list their faces (polyhedra and "
                        "polygons); a shape's faces come from its table, wound one way")
    if t.dimension == 1:
        raise ValueError("polylines have no winding to check")
    if t.dimension == 2:
        return _check_loops(t)
    edges = t.edges()
    n = t.num_cells
    sizes = tack.empty((n,), tack.i32)
    starts = tack.empty((n,), t.id_dtype)
    if n:
        _cell_sizes(t.cell_offsets, t.cell_faces, t.face_offsets, sizes, n)
    total = exclusive_scan(sizes, starts, n) if n else 0
    bad = tack.zeros((n,), tack.i32)
    if total:
        keys = tack.empty((total,), tack.vec(2, t.id_dtype))
        directions = tack.empty((total,), tack.i32)
        _directed_edges(t.cell_offsets, t.cell_faces, t.cell_face_sides, t.face_offsets,
                        edges.face_edge, edges.face_edge_sign, starts, keys, directions, n)
        order = key_order(keys, n, t.id_dtype)
        offsets, runs = run_offsets(keys, order)
        _unbalanced(keys, directions, order, offsets, bad, runs)
    return Set._from_flags(t, "cells", bad)

merge_points 

merge_points(data, tolerance=0.0, fields=None)

data with each set of coincident points made one: points at most tolerance apart (inclusive; 0, the default, merges only points of equal coordinates) join, and so do chains of them, so a set may span more than tolerance. A set becomes the point of its smallest id, which keeps that point's position, point fields, ghost marks and global id; the merged points keep the order of those ids, and the cells are the same cells on them, a cell's corners that merged kept as they are. Cell fields are unchanged.

For meshes written without merging their points (each piece of a simulation's output storing a point once for every cell or block that uses it), which nothing can stitch by matching points across cells or pieces until each point is one. tolerance is in the dataset's spatial unit, or a quantity ("1e-9 m"); it should be far below the distance between points that are not copies, since every point compares itself with every point that close. Takes an unstructured topology; a structured grid's points are its lattice, and polyhedral, polygonal and polyline topologies are not handled yet. Hidden cells are dropped first, as every filter that makes cells drops them.

Source code in packages/tack-vis/src/tack/dataset/algorithms/merge_points.py 
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
@on_inputs_device
def merge_points(data, tolerance=0.0, fields=None):
    """``data`` with each set of coincident points made one: points at most
    ``tolerance`` apart (inclusive; 0, the default, merges only points of equal
    coordinates) join, and so do chains of them, so a set may span more than
    ``tolerance``. A set becomes the point of its smallest id, which keeps that
    point's position, point fields, ghost marks and global id; the merged points
    keep the order of those ids, and the cells are the same cells on them, a
    cell's corners that merged kept as they are. Cell fields are unchanged.

    For meshes written without merging their points (each piece of a
    simulation's output storing a point once for every cell or block that
    uses it), which nothing can stitch by matching points across cells or
    pieces until each point is one. ``tolerance`` is in the dataset's spatial
    unit, or a quantity (``"1e-9 m"``); it should be far below the distance
    between points that are not copies, since every point compares itself with
    every point that close. Takes an unstructured topology; a structured grid's
    points are its lattice, and polyhedral, polygonal and polyline topologies
    are not handled yet. Hidden cells are dropped first, as every filter that
    makes cells drops them."""
    tolerance = _length(data, tolerance, "merge_points's tolerance")
    if not tolerance >= 0:
        raise ValueError(f"tolerance must be at least zero, not {tolerance!r}")
    if isinstance(data.topology, StructuredTopology):
        raise TypeError("merge_points takes an unstructured topology: a structured "
                        "grid's points are its lattice")
    if not isinstance(data.topology, UnstructuredTopology):
        raise NotImplementedError(f"merge_points of a {type(data.topology).__name__} is "
                                  f"not supported yet")
    if not isinstance(data.geometry.space, H1) or data.geometry.space.order != 1:
        raise TypeError("merge_points needs an order-1 H1 geometry: one position per point")
    data = visible(data)
    topology = data.topology
    n = data.num_points
    dtype = ids.at_least(tack.i32, n + 1)
    if dtype not in get_backend().supported_atomic_dtypes:
        raise NotImplementedError(f"merging {n} points needs {dtype} atomics, which the "
                                  f"{get_backend().name} backend does not have")
    # Union-find over the points within the tolerance, through the locator's
    # bins: each pair hooks the larger root under the smaller (an atomic
    # minimum), until nothing changes, so a set's root is its smallest point
    # whatever the order of the hooks.
    labels = tack.arange(n, dtype=dtype)
    if n:
        grid = cached_point_locator(data).grid()
        every = tack.empty((n,), tack.u8)
        every.fill(1)
        changed = tack.zeros((1,), tack.i32)
        while True:
            changed.fill(0)
            _hook_near(grid, 0, every, float(tolerance), labels, changed)
            _jump(labels)
            if not int(changed.numpy()[0]):
                break
    flags = tack.empty((max(n, 1),), dtype)
    index = tack.empty((max(n, 1),), dtype)
    count = 0
    if n:
        _roots(labels, flags)
        count = int(exclusive_scan(flags, index, n))
    new_id = tack.empty((max(n, 1),), dtype)
    kept = tack.empty((max(count, 1),), dtype)
    if n:
        _numbered(labels, index, new_id, kept)
    length = arrays.size_of(topology.connectivity)
    connectivity = tack.empty((length,), data.id_dtype)
    if length:
        _renumbered(topology.connectivity, new_id, connectivity)
    merged = UnstructuredTopology(topology.types, topology.offsets, connectivity,
                                  num_points=count)
    kept = kept if count else tack.empty((0,), dtype)
    return carry(data, merged, points=Take(kept), cells=Same(), fields=fields)

Changing cell types

The same points and region, in other cells.

Algorithm Gives Takes VTK Viskores
tetrahedralize 3D cells as tetrahedra that meet face to face dataset vtkDataSetTriangleFilter Tetrahedralize
triangulate 2D cells as triangles dataset vtkTriangleFilter Triangulate
as_polyhedra a shape-based mesh as a polyhedral one dataset — —
as_polygons 2D cells as a polygonal mesh dataset — —
shrink every cell shrunk toward its centroid dataset vtkShrinkFilter Shrink

tetrahedralize 

tetrahedralize(data, fields=None)

Every 3D cell as tetrahedra that meet face to face: two cells sharing a face split it along one diagonal, as VTK's vtkDataSetTriangleFilter does. A structured grid's hexahedra take 5 each, the two mirror splits alternating with the parity of (i, j, k) in the whole grid (the topology's offset added); any other cell is the cone from its corner of smallest point id -- global id, where data has global ids -- over the faces without it, each quadrilateral face split from its own smallest corner (a hexahedron or voxel in 6, a wedge in 3, a pyramid in 2), which fills a convex cell. So the pieces of a partitioned dataset split the faces they share alike too. Other cells are left out. The points are the same; cell fields go to each cell's tetrahedra.

Source code in packages/tack-vis/src/tack/dataset/algorithms/tetrahedralize.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
@on_inputs_device
def tetrahedralize(data, fields=None):
    """Every 3D cell as tetrahedra that meet face to face: two cells sharing a face
    split it along one diagonal, as VTK's vtkDataSetTriangleFilter does. A
    structured grid's hexahedra take 5 each, the two mirror splits alternating
    with the parity of (i, j, k) in the whole grid (the topology's ``offset``
    added); any other cell is the cone from its corner of smallest point id --
    global id, where ``data`` has global ids -- over the faces without it, each
    quadrilateral face split from its own smallest corner (a hexahedron or
    voxel in 6, a wedge in 3, a pyramid in 2), which fills a convex cell. So
    the pieces of a partitioned dataset split the faces they share alike too. Other cells are left out. The
    points are the same; cell fields go to each cell's tetrahedra."""
    from tack.dataset.ghosts import visible

    data = visible(data)
    if not getattr(data.topology, "reference_cells", True):
        raise NotImplementedError("tetrahedralize takes cells of the linear shapes")
    return _split(data, "tetrahedra", fields)

triangulate 

triangulate(data, fields=None)

Every 2D cell as triangles: a quad or pixel in 2, a polygon in a fan from its first point; other cells are left out. The points are the same; cell fields go to each cell's triangles.

Source code in packages/tack-vis/src/tack/dataset/algorithms/triangulate.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
@on_inputs_device
def triangulate(data, fields=None):
    """Every 2D cell as triangles: a quad or pixel in 2, a polygon in a fan from its
    first point; other cells are left out. The points are the same; cell fields
    go to each cell's triangles."""
    from tack.dataset.ghosts import visible

    data = visible(data)
    t = data.topology
    if getattr(t, "reference_cells", True):
        return _split(data, "triangles", fields)
    if t.dimension != 2:
        raise NotImplementedError(f"triangulate takes 2D cells, not {t.dimension}D ones")
    loop_offsets, loop_points = t.loops()
    n = data.num_cells
    odt = ids.for_output(data, 3 * size_of(loop_points))
    counts = tack.empty((n,), tack.i32)
    starts = tack.empty((n,), odt)
    if n:
        _fan_counts(loop_offsets, counts)
    count = exclusive_scan(counts, starts, n) if n else 0
    connectivity = tack.empty((3 * count,), odt)
    sources = tack.empty((count,), data.id_dtype)
    if count:
        _fans(loop_offsets, loop_points, starts, connectivity, sources)
    return _pieces_topology(data, connectivity, sources, count, shapes.TRIANGLE, 3, fields)

as_polyhedra 

as_polyhedra(data, fields=None)

A shape-based dataset as a polyhedral one, every cell a polyhedron.

Its derived faces already hold every face once, in side 0's outward order (a voxel's pixel faces as quads), with the side of each cell's reference, so they become the polyhedral faces as they stand: face ids, and so fields on faces, are unchanged. Edges come out numbered as before. Cells must be 3D, and neighbours wound alike: an inside-out cell beside an upright one is refused, since their shared face would have them on one side. Constant and H1 order-1 fields, values on points, cells, faces and edges, and sets carry over; fields needing a reference element are dropped.

Source code in packages/tack-vis/src/tack/dataset/algorithms/as_polyhedra.py 
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
@on_inputs_device
def as_polyhedra(data, fields=None):
    """A shape-based dataset as a polyhedral one, every cell a polyhedron.

    Its derived faces already hold every face once, in side 0's outward order
    (a voxel's pixel faces as quads), with the side of each cell's reference, so
    they become the polyhedral faces as they stand: face ids, and so fields on
    faces, are unchanged. Edges come out numbered as before. Cells must be 3D,
    and neighbours wound alike: an inside-out cell beside an upright one is
    refused, since their shared face would have them on one side.
    ``Constant`` and ``H1`` order-1 fields, values on points, cells, faces and
    edges, and sets carry over; fields needing a reference element are dropped.
    """
    from tack.dataset.spaces import H1

    if not isinstance(data.geometry.space, H1) or data.geometry.space.order != 1:
        raise TypeError("as_polyhedra needs an order-1 H1 geometry: positions per point")
    topology = data.topology
    for group in topology.groups():
        if group.count and group.shape.DIMENSION != 3:
            raise ValueError(f"{group.shape.__name__} cells are not 3D: a polyhedral "
                             "topology is of 3D cells (as_polygons takes 2D ones)")
    faces = topology.faces()
    same_way = tack.zeros((1,), tack.i32)
    for group in topology.groups():
        if group.count:
            _same_way(data.domain_view("cells", group), same_way)
    if same_way[0]:
        raise ValueError(f"{same_way[0]} faces are wound the same way by both their cells "
                         "(an inside-out cell beside an upright one, as cell_geometry's "
                         "negative measures show); a polyhedral topology needs each face's "
                         "cells on opposite sides of it")
    nf = faces.num_faces
    idt = topology.id_dtype
    # The faces' points: 4 a face, past 12 a cell when faces are not shared (a
    # hexahedron has 24), and the polyhedral topology holds twice that.
    fdt = ids.at_least(idt, 8 * nf)
    sizes = tack.empty((nf,), tack.i32)
    starts = tack.empty((nf + 1,), fdt)
    if nf:
        _face_sizes(faces.kinds, sizes)
    total = exclusive_scan(sizes, starts, nf, close=True)
    face_points = tack.empty((total,), fdt)
    if nf:
        _face_rows(faces.rows, starts, sizes, face_points)

    n = data.num_cells
    counts = tack.zeros((n,), tack.i32)
    for group in topology.groups():
        if group.count:
            _cell_face_counts(group.view(), counts)
    cell_offsets = tack.empty((n + 1,), idt)
    entries = exclusive_scan(counts, cell_offsets, n, close=True)
    cell_faces = tack.empty((entries,), idt)
    cell_face_sides = tack.empty((entries,), tack.u8)
    for group in topology.groups():
        if group.count:
            _cell_face_entries(data.domain_view("cells", group), cell_offsets, cell_faces,
                               cell_face_sides)
    out = PolyhedralTopology(starts, face_points, cell_offsets, cell_faces, cell_face_sides,
                             num_points=data.num_points)

    # The same points, cells and faces (in the shape path's numbering), so the
    # sets still apply; edges too, both paths numbering them by (low, high)
    # point ids.
    return carry(data, out, points=Same(), cells=Same(), faces=Same(), edges=Same(),
                 fields=fields)

as_polygons 

as_polygons(data, fields=None)

A shape-based dataset of 2D cells -- triangles, quads, pixels -- as a polygonal one: each cell's points in order become its loop. Fields carry as as_polyhedra's do, but values on faces and edges, whose numbering the polygonal topology derives afresh, are dropped.

Source code in packages/tack-vis/src/tack/dataset/algorithms/as_polygons.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
@on_inputs_device
def as_polygons(data, fields=None):
    """A shape-based dataset of 2D cells -- triangles, quads, pixels -- as a polygonal
    one: each cell's points in order become its loop. Fields carry as
    ``as_polyhedra``'s do, but values on faces and edges, whose numbering the
    polygonal topology derives afresh, are dropped."""

    topology = data.topology
    for group in topology.groups():
        if group.count and group.shape.DIMENSION != 2:
            raise ValueError(f"{group.shape.__name__} cells are not 2D")
    n = data.num_cells
    idt = topology.id_dtype
    sizes = tack.zeros((n,), tack.i32)
    for group in topology.groups():
        if group.count:
            _polygon_sizes(group.view(), sizes)
    offsets = tack.empty((n + 1,), idt)
    total = exclusive_scan(sizes, offsets, n, close=True)
    loops = tack.empty((total,), idt)
    for group in topology.groups():
        if group.count:
            _polygon_loops(group.view(), offsets, loops)
    out = PolygonalTopology(offsets, loops, num_points=data.num_points)
    # The same points and cells; edge numbering is derived afresh.
    return carry(data, out, points=Same(), cells=Same(), fields=fields)

shrink 

shrink(data, factor=0.5, fields=None)

Every cell shrunk toward its centroid (the mean of its corners) by factor, as Viskores' Shrink and VTK's vtkShrinkFilter. Viskores gives each cell its own copies of its points; here that is an L2 geometry on the same topology, and point data becomes L2 the same way (each cell's own copy of its corners' values). Everything else is kept, the topology's faces and edges with it.

Source code in packages/tack-vis/src/tack/dataset/algorithms/shrink.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
@on_inputs_device
def shrink(data, factor=0.5, fields=None):
    """Every cell shrunk toward its centroid (the mean of its corners) by ``factor``,
    as Viskores' Shrink and VTK's vtkShrinkFilter. Viskores gives each cell its
    own copies of its points; here that is an ``L2`` geometry on the same
    topology, and point data becomes ``L2`` the same way (each cell's own copy
    of its corners' values). Everything else is kept, the topology's faces and
    edges with it."""

    if not 0 <= factor <= 1:
        raise ValueError(f"factor must lie in [0, 1], not {factor!r}")
    corners = project(data, data.geometry, L2(data))
    if data.num_cells:
        _shrink_corners(corners.values, corners.space.offsets, float(factor), data.num_cells)
    from tack.dataset.carry import selected

    out = {}
    for name in selected(data, fields):
        f = data.fields[name]
        space = f.space
        if isinstance(space, Values) and space.on == "points":
            f, space = Field(H1(data), f.values, f.unit), H1(data)
        out[name] = project(data, f, L2(data)) if isinstance(space, H1) else f
    from tack.dataset.ghosts import keep_marks

    return keep_marks(data, DataSet(data.topology, corners, fields=out, sets=dict(data.sets)))

Moving the geometry

New positions for the same topology.

Algorithm Gives Takes VTK Viskores
transform an affine map (4x4 or 3x4 matrix) dataset vtkTransformFilter PointTransform
warp positions moved along a vector field or vector dataset vtkWarpVector, vtkWarpScalar Warp
cylindrical Cartesian to cylindrical coordinates, or back dataset — CylindricalCoordinateTransform
spherical Cartesian to spherical coordinates, or back dataset — SphericalCoordinateTransform

transform 

transform(data, matrix, fields=None)

The geometry under an affine map: a 4x4 matrix (its last row 0 0 0 1) or 3x4, applied to each position as a column (Viskores' PointTransform, VTK's vtkTransformFilter). On a dataset with an origin the map moves the origin, in double precision, and only its linear part acts on the geometry's values, so a rotation far from zero loses nothing.

Source code in packages/tack-vis/src/tack/dataset/algorithms/transform.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
@on_inputs_device
def transform(data, matrix, fields=None):
    """The geometry under an affine map: a 4x4 matrix (its last row 0 0 0 1) or 3x4,
    applied to each position as a column (Viskores' PointTransform, VTK's
    vtkTransformFilter). On a dataset with an origin the map moves the origin,
    in double precision, and only its linear part acts on the geometry's values,
    so a rotation far from zero loses nothing."""
    m = np.asarray(matrix, float)
    if m.shape == (4, 4):
        if not np.allclose(m[3], [0, 0, 0, 1]):
            raise ValueError("a projective matrix (last row not 0 0 0 1) is not affine")
        m = m[:3]
    if m.shape != (3, 4):
        raise ValueError(f"transform takes a 4x4 or 3x4 matrix, not {m.shape}")
    origin = None
    if any(data.origin):
        origin = m[:, :3] @ np.asarray(data.origin) + m[:, 3]
        m = np.hstack([m[:, :3], np.zeros((3, 1))])
    n = size_of(data.geometry.values)
    out = _vectors(n, data.dtype)
    if n:
        _affine(materialize(data.geometry.values), *map(float, m.reshape(-1)), out)
    out = _moved(data, out, fields)
    if origin is not None:
        out.origin = origin
    return out

warp 

warp(
    data, direction, scale=1.0, scale_by=None, fields=None
)

The geometry moved along direction -- a 3-vector field in the geometry's space (or its name), or one 3-vector for every point -- times scale and, if given, a scalar field scale_by in the same space (Viskores' Warp; VTK's vtkWarpVector and vtkWarpScalar).

Source code in packages/tack-vis/src/tack/dataset/algorithms/warp.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
@on_inputs_device
def warp(data, direction, scale=1.0, scale_by=None, fields=None):
    """The geometry moved along ``direction`` -- a 3-vector field in the geometry's
    space (or its name), or one 3-vector for every point -- times ``scale``
    and, if given, a scalar field ``scale_by`` in the same space (Viskores'
    Warp; VTK's vtkWarpVector and vtkWarpScalar)."""
    geometry = data.geometry
    n = size_of(geometry.values)
    dtype = data.dtype
    factors = (scale_by if not isinstance(scale_by, str) else data.fields[scale_by])
    if factors is None:
        factors = tack.full((n,), 1.0, dtype) if n else tack.empty((0,), dtype)
    else:
        if factors.space is not geometry.space or width_of(factors.values):
            raise ValueError("scale_by must be a scalar field in the geometry's space")
        factors = materialize(factors.values)
    out = _vectors(n, dtype)
    if isinstance(direction, (str, Field)):
        direction = data.fields[direction] if isinstance(direction, str) else direction
        if direction.space is not geometry.space or width_of(direction.values) != 3:
            raise ValueError("direction must be a 3-vector field in the geometry's space")
        if n:
            _warp_by_field(materialize(geometry.values), materialize(direction.values),
                           factors, float(scale), out)
    else:
        d = np.asarray(direction, float).reshape(-1)
        if d.size != 3:
            raise ValueError("a constant direction is three numbers")
        if n:
            _warp_constant(materialize(geometry.values), *map(float, d), factors,
                           float(scale), out)
    return _moved(data, out, fields)

cylindrical 

cylindrical(data, inverse=False, fields=None)

The geometry's Cartesian positions as cylindrical (r, theta, z), or back with inverse -- Viskores' CylindricalCoordinateTransform: theta from asin(y / r), in [-pi/2, 3pi/2). The positions are the world's (the origin added), and the result's origin is zero and its spatial unit unknown, an angle and lengths being no one unit. Fields come along.

Source code in packages/tack-vis/src/tack/dataset/algorithms/cylindrical.py 
 9
10
11
12
13
14
15
16
@on_inputs_device
def cylindrical(data, inverse=False, fields=None):
    """The geometry's Cartesian positions as cylindrical ``(r, theta, z)``, or back
    with ``inverse`` -- Viskores' CylindricalCoordinateTransform: ``theta`` from
    ``asin(y / r)``, in ``[-pi/2, 3pi/2)``. The positions are the world's
    (the origin added), and the result's origin is zero and its spatial unit
    unknown, an angle and lengths being no one unit. Fields come along."""
    return _recoordinate(data, _from_cylindrical if inverse else _to_cylindrical, fields)

spherical 

spherical(data, inverse=False, fields=None)

The geometry's Cartesian positions as spherical (r, theta, phi) -- theta from the z axis, phi in [0, 2pi) -- or back with inverse, as Viskores' SphericalCoordinateTransform. The positions are the world's (the origin added), and the result's origin is zero and its spatial unit unknown, an angle and lengths being no one unit. Fields come along.

Source code in packages/tack-vis/src/tack/dataset/algorithms/spherical.py 
 9
10
11
12
13
14
15
16
@on_inputs_device
def spherical(data, inverse=False, fields=None):
    """The geometry's Cartesian positions as spherical ``(r, theta, phi)`` -- ``theta``
    from the z axis, ``phi`` in ``[0, 2pi)`` -- or back with ``inverse``, as
    Viskores' SphericalCoordinateTransform. The positions are the world's (the
    origin added), and the result's origin is zero and its spatial unit
    unknown, an angle and lengths being no one unit. Fields come along."""
    return _recoordinate(data, _from_spherical if inverse else _to_spherical, fields)

Measures and quality

Fields describing the mesh itself.

Algorithm Gives Takes VTK Viskores
cell_geometry each cell's volume (area, length) and centroid dataset vtkCellSizeFilter CellMeasures
face_geometry each face's normal and area dataset — —
cell_centers each cell's parametric center (or cell point) dataset vtkCellCenters —
face_centers each face's area centroid dataset — —
edge_lengths each edge's length dataset — —
cell_quality scaled Jacobian, Jacobian ratio, condition number, edge ratio, edge lengths dataset vtkMeshQuality MeshQuality
point_elevation each position along a line, mapped to a range dataset vtkElevationFilter PointElevation
implicit_values an implicit function at the geometry's values dataset vtkSampleImplicitFunctionFilter —
cell_ids each cell's id (provenance through filters) dataset vtkGenerateIds —
point_ids each point's id (provenance through filters) dataset vtkGenerateIds —

cell_geometry 

cell_geometry(data)

CellGeometry(measure, centroid) of each cell: fields on cells. A 3D cell's volume, by the divergence theorem over its faces -- polyhedral or of a shape -- with faces taken as fans from their points' means, so a cell with non-planar faces is measured as that polyhedron (a hexahedron with warped faces is not quite its trilinear volume); its faces must be wound out of it. A 2D cell's area -- triangle, quad, pixel or polygon -- along its own normal (Newell's), from its edges going round it. A line's or a polyline's length, and the centroid of its segments; a vertex's point, with measure 0.

Source code in packages/tack-vis/src/tack/dataset/algorithms/cell_geometry.py 
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
@on_inputs_device
def cell_geometry(data):
    """``CellGeometry(measure, centroid)`` of each cell: fields on cells. A 3D cell's volume, by the
    divergence theorem over its faces -- polyhedral or of a shape -- with faces
    taken as fans from their points' means, so a cell with non-planar faces is
    measured as that polyhedron (a hexahedron with warped faces is not quite its
    trilinear volume); its faces must be wound out of it. A 2D cell's area --
    triangle, quad, pixel or polygon -- along its own normal (Newell's), from
    its edges going round it. A line's or a polyline's length, and the
    centroid of its segments; a vertex's point, with measure 0."""
    if data.topology.dimension > 1 and not getattr(data.topology, "reference_cells", True):
        data.topology.faces()                      # polyhedra and polygons walk theirs
    n = data.num_cells
    volumes = tack.empty((n,), data.dtype)
    centroids = _vectors(n, data.dtype)
    for group in data.launch_groups("cells", []):
        if group.count and group.shape.DIMENSION == 3:
            _cell_geometry(data.domain_view("cells", group), volumes, centroids)
        elif group.count and group.shape.DIMENSION < data.topology.dimension and \
                group.shape in _ROUND:
            _corner_polygon_geometry(data.domain_view("cells", group),
                                     _table(("round", group.shape), _ROUND[group.shape]),
                                     volumes, centroids)
        elif group.count and group.shape.DIMENSION == 2:
            _polygon_geometry(data.domain_view("cells", group), volumes, centroids)
        elif group.count and group.shape is shapes.PolyLine:
            _line_geometry(data.domain_view("cells", group), volumes, centroids)
        elif group.count:
            _corner_geometry(data.domain_view("cells", group), volumes, centroids)
    return CellGeometry(Field(Values(data, "cells"), volumes, _measure_unit(data)),
                        Field(Values(data, "cells"), centroids, data.unit))

CellGeometry

Bases: NamedTuple

cell_geometry's fields, on the cells.

face_geometry 

face_geometry(data)

FaceGeometry(normal, area): fields on faces. A face's normal points out of its side 0. A 3D mesh's faces are polygons of any size, the normal the mean one on a non-planar face; a 2D mesh's faces are edges, the normal in the plane of side 0's cell and the "area" the edge's length; a polyline mesh's faces are the polylines' ends, the normal along the line out of side 0 and the "area" 1.

Source code in packages/tack-vis/src/tack/dataset/algorithms/face_geometry.py 
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
@on_inputs_device
def face_geometry(data):
    """``FaceGeometry(normal, area)``: fields on faces. A face's normal points out of
    its side 0. A 3D mesh's faces are polygons of any size, the normal the mean
    one on a non-planar face; a 2D mesh's faces are edges, the normal in the
    plane of side 0's cell and the "area" the edge's length; a polyline mesh's
    faces are the polylines' ends, the normal along the line out of side 0 and
    the "area" 1."""
    faces = data.topology.faces()
    dtype = data.dtype
    normals = _vectors(faces.num_faces, dtype)
    areas = tack.empty((faces.num_faces,), dtype)
    if data.topology.dimension == 2:
        for_each(_edge_geometry, data, "cells", normals, areas)
    elif data.topology.dimension == 1:
        for_each(_end_geometry, data, "cells", normals, areas)
    else:
        for_each(_face_geometry, data, "faces", normals, areas)
    # A face of a k-dimensional mesh measures k - 1 dimensions; a polyline end's
    # "area" is 1.
    area = units.derived((data.unit, data.topology.dimension - 1))
    return FaceGeometry(Field(Values(data, "faces"), normals, "dimensionless"),
                        Field(Values(data, "faces"), areas, area))

FaceGeometry

Bases: NamedTuple

face_geometry's fields, on the faces.

cell_centers 

cell_centers(data)

Each cell's center, where the geometry maps its parametric center: a field on cells. That is values_at_centers(data, data.geometry), and not the centroid on a cell that is not affine; cell_geometry(data).centroid is the centroid, on every topology, polyhedra included.

A polyhedron, polygon or polyline has no parametric center: its center is its cell point, the mean of its distinct points, where its interpolation fans its faces from. That is not the centroid either, unless the points are spread evenly.

Source code in packages/tack-vis/src/tack/dataset/algorithms/cell_centers.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
@on_inputs_device
def cell_centers(data):
    """Each cell's center, where the geometry maps its parametric center: a field on cells.
    That is ``values_at_centers(data, data.geometry)``, and not the centroid on a
    cell that is not affine; ``cell_geometry(data).centroid`` is the centroid, on
    every topology, polyhedra included.

    A polyhedron, polygon or polyline has no parametric center: its center is its
    *cell point*, the mean of its distinct points, where its interpolation fans
    its faces from. That is not the centroid either, unless the points are
    spread evenly."""
    out = _vectors(data.num_cells, data.dtype)
    if not getattr(data.topology, "reference_cells", True):
        from tack.dataset.algorithms.project import _distinct_point_means

        for_each(_distinct_point_means, data, "cells", materialize(data.geometry.values), out)
    else:
        for_each(_centers, data, "cells", out)
    return Field(Values(data, "cells"), out, data.unit)

face_centers 

face_centers(data)

Each face's area centroid: a field on faces. The face is fanned from its points' mean, so both its cells -- which store it once -- see the same fan. A 2D mesh's face, an edge, has its midpoint, and a polyline's end its point.

Source code in packages/tack-vis/src/tack/dataset/algorithms/face_centers.py 
37
38
39
40
41
42
43
44
@on_inputs_device
def face_centers(data):
    """Each face's area centroid: a field on faces. The face is fanned from its
    points' mean, so both its cells -- which store it once -- see the same fan. A
    2D mesh's face, an edge, has its midpoint, and a polyline's end its point."""
    out = _vectors(data.topology.faces().num_faces, data.dtype)
    for_each(_face_centers, data, "faces", out)
    return Field(Values(data, "faces"), out, data.unit)

edge_lengths 

edge_lengths(data)

Each edge's length, the distance between its ends in the spatial unit: a field on edges (Values(data, "edges")), numbered as data.topology.edges() numbers them. An order-2 geometry's edges are measured straight. Polylines have no derived edges.

Source code in packages/tack-vis/src/tack/dataset/algorithms/edge_lengths.py 
17
18
19
20
21
22
23
24
25
@on_inputs_device
def edge_lengths(data):
    """Each edge's length, the distance between its ends in the spatial unit: a
    field on edges (``Values(data, "edges")``), numbered as
    ``data.topology.edges()`` numbers them. An order-2 geometry's edges are
    measured straight. Polylines have no derived edges."""
    out = tack.empty((data.topology.edges().num_edges,), data.dtype)
    for_each(_edge_lengths, data, "edges", out)
    return Field(Values(data, "edges"), out, data.unit)

cell_quality 

cell_quality(data, measures=None)

How well shaped each cell is: a dict of cell fields, one per measure named in measures (all of MEASURES by default). Each is 1 on the ideal cell of its shape (equilateral; a cube, a regular tetrahedron, an equilateral wedge and pyramid), as Verdict (VTK's vtkMeshQuality) defines it where it does:

  • scaled_jacobian: the least, over the cell's corners (and a hexahedron's center), of the determinant of its edges there, each scaled to unit length, relative to the ideal cell's: 1 ideal, 0 degenerate, negative inverted.
  • jacobian_ratio: the least corner determinant over the greatest; 1 for a cell whose Jacobian is constant (every simplex), negative where one corner is inverted.
  • condition: the greatest, over the corners, condition number of the corner's edges against the ideal cell's (|T| |T^-1| / d, T = A W^-1, Frobenius norms): 1 ideal, growing as the cell distorts, infinite (the largest float) where a corner is degenerate or inverted.
  • edge_ratio: the longest edge over the shortest.
  • min_edge, max_edge: the shortest and longest edge, in the spatial unit.

The corners' edges are the geometry's Jacobian there applied to the reference cell's edge directions, so an order-2 (curved) cell is measured along its curved edges at its corners. A pyramid, which has no Jacobian at its apex, is measured as Verdict does, by its four tetrahedra over both base diagonals. Polyhedra and polygons, which have no reference cell, are measured by the fan tetrahedra and triangles Tack evaluates them in, each against the regular tetrahedron or equilateral triangle: a valid star-shaped polyhedron has a positive scaled Jacobian, though an ideal-looking one (a cube as a polyhedron) scores below 1. Lines and vertices have no Jacobian measures (NaN); a line's edges are its segments. Hidden cells are measured like the rest.

Matched to VTK 9.7.1's vtkMeshQuality to rounding for the scaled Jacobian of every linear shape, the condition number of triangles, quads, tetrahedra and hexahedra, and the edge ratio of triangles, quads, tetrahedra and wedges. VTK's wedge condition number and hexahedron edge ratio differ (see the tests).

Source code in packages/tack-vis/src/tack/dataset/algorithms/cell_quality.py 
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
@on_inputs_device
def cell_quality(data, measures=None):
    """How well shaped each cell is: a dict of cell fields, one per measure named in
    ``measures`` (all of ``MEASURES`` by default). Each is 1 on the ideal cell of its
    shape (equilateral; a cube, a regular tetrahedron, an equilateral wedge and
    pyramid), as Verdict (VTK's ``vtkMeshQuality``) defines it where it does:

    - ``scaled_jacobian``: the least, over the cell's corners (and a hexahedron's
      center), of the determinant of its edges there, each scaled to unit length,
      relative to the ideal cell's: 1 ideal, 0 degenerate, negative inverted.
    - ``jacobian_ratio``: the least corner determinant over the greatest; 1 for a
      cell whose Jacobian is constant (every simplex), negative where one corner
      is inverted.
    - ``condition``: the greatest, over the corners, condition number of the
      corner's edges against the ideal cell's (``|T| |T^-1| / d``, ``T = A W^-1``,
      Frobenius norms): 1 ideal, growing as the cell distorts, infinite (the
      largest float) where a corner is degenerate or inverted.
    - ``edge_ratio``: the longest edge over the shortest.
    - ``min_edge``, ``max_edge``: the shortest and longest edge, in the spatial
      unit.

    The corners' edges are the geometry's Jacobian there applied to the reference
    cell's edge directions, so an order-2 (curved) cell is measured along its
    curved edges at its corners. A pyramid, which has no Jacobian at its apex, is
    measured as Verdict does, by its four tetrahedra over both base diagonals.
    Polyhedra and polygons, which have no reference cell, are measured by the fan
    tetrahedra and triangles Tack evaluates them in, each against the regular
    tetrahedron or equilateral triangle: a valid star-shaped polyhedron has a
    positive scaled Jacobian, though an ideal-looking one (a cube as a polyhedron)
    scores below 1. Lines and vertices have no Jacobian measures (NaN); a line's
    edges are its segments. Hidden cells are measured like the rest.

    Matched to VTK 9.7.1's ``vtkMeshQuality`` to rounding for the scaled Jacobian
    of every linear shape, the condition number of triangles, quads, tetrahedra
    and hexahedra, and the edge ratio of triangles, quads, tetrahedra and wedges.
    VTK's wedge condition number and hexahedron edge ratio differ (see the
    tests)."""
    names = list(MEASURES if measures is None else
                 [measures] if isinstance(measures, str) else measures)
    unknown = [m for m in names if m not in MEASURES]
    if unknown:
        raise ValueError(f"no measure {', '.join(map(repr, unknown))}: the measures are "
                         f"{', '.join(MEASURES)}")
    n = data.num_cells
    dtype = data.dtype
    big = float(np.finfo(dtype.numpy_dtype).max)
    out = {m: tack.empty((n,), dtype) for m in ("scaled_jacobian", "jacobian_ratio",
                                                "condition", "min_edge", "max_edge")}
    if data.topology.dimension > 1 and not getattr(data.topology, "reference_cells", True):
        data.topology.faces()                      # polyhedra and polygons walk theirs
    jacobian = {"scaled_jacobian", "jacobian_ratio", "condition"} & set(names)
    for group in data.launch_groups("cells", []):
        if not group.count:
            continue
        view = data.domain_view("cells", group)
        shape = group.shape
        sj, ratio, cond = out["scaled_jacobian"], out["jacobian_ratio"], out["condition"]
        lo, hi = out["min_edge"], out["max_edge"]
        if shape in _REFERENCE:
            if jacobian:
                pcs, directions, norms, inverses, corner, _ = _tables(shape, dtype)
                _reference_quality(view, shape.DIMENSION, pcs, directions, norms, inverses,
                                   corner, big, sj, ratio, cond)
            _edges_by_table(view, _tables(shape, dtype)[5], lo, hi)
        elif shape is shapes.Pyramid:
            if jacobian:
                norms, inverses = _pyramid_tables(dtype)
                _pyramid_quality(view, norms, inverses, big, sj, ratio, cond)
            edges = _table("pyramid edges", _PYRAMID_EDGES)
            _edges_by_table(view, edges, lo, hi)
        elif shape is shapes.Polyhedron:
            if jacobian:
                w = tack.array(np.linalg.inv(_REGULAR).astype(dtype.numpy_dtype).reshape(-1))
                _polyhedron_quality(view, w, big, sj, ratio, cond)
            _edges_by_sides(view, lo, hi)
        elif shape is shapes.Polygon:
            if jacobian:
                _polygon_quality(view, big, sj, ratio, cond)
            _edges_by_sides(view, lo, hi)
        else:
            # Lines, polylines, vertices: no Jacobian measures.
            _no_jacobian(view, float("nan"), sj, ratio, cond)
            if shape is shapes.PolyLine:
                _edges_of_lines(view, lo, hi)
            else:
                _edges_by_table(view, _table(("line edges", shape is shapes.Line),
                                             [0, 1] if shape is shapes.Line else [0, 0]), lo, hi)
    if "edge_ratio" in names:
        out["edge_ratio"] = tack.empty((n,), dtype)
        if n:
            _ratio(out["min_edge"], out["max_edge"], big, float("nan"), out["edge_ratio"])
    space = Values(data, "cells")
    units = {"min_edge": data.unit, "max_edge": data.unit}
    return {m: Field(space, out[m], units.get(m, "dimensionless")) for m in names}

point_elevation 

point_elevation(
    data,
    low_point=(0, 0, 0),
    high_point=(0, 0, 1),
    range=(0.0, 1.0),
)

Each geometry value's position along the line from low_point to high_point, clamped to it and mapped to range (VTK's vtkElevationFilter, Viskores' PointElevation): a field in the geometry's space. The two points are world positions, taken relative to the dataset's origin in double precision.

Source code in packages/tack-vis/src/tack/dataset/algorithms/point_elevation.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
@on_inputs_device
def point_elevation(data, low_point=(0, 0, 0), high_point=(0, 0, 1), range=(0.0, 1.0)):
    """Each geometry value's position along the line from ``low_point`` to
    ``high_point``, clamped to it and mapped to ``range`` (VTK's
    vtkElevationFilter, Viskores' PointElevation): a field in the geometry's
    space. The two points are world positions, taken relative to the dataset's
    origin in double precision."""
    low = np.asarray(low_point, float) - np.asarray(data.origin)
    direction = np.asarray(high_point, float) - np.asarray(data.origin) - low
    length2 = float(direction @ direction)
    if length2 == 0:
        raise ValueError("the low and high points must differ")
    n = size_of(data.geometry.values)
    out = _scalars(n, data.dtype)
    if n:
        _elevations(materialize(data.geometry.values), *map(float, low), *map(float, direction),
                    length2, float(range[0]), float(range[1] - range[0]), out)
    return Field(data.geometry.space, out)

implicit_values 

implicit_values(data, function)

function, an implicit function (tack.dataset.implicit: Plane, Sphere, Box, ...), at the geometry's values: a field in the geometry's space and type, negative inside the function's region. On an H1 geometry that is one value per point. The function is in world positions, taken relative to the dataset's origin in double precision.

Source code in packages/tack-vis/src/tack/dataset/algorithms/implicit_values.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
@on_inputs_device
def implicit_values(data, function):
    """``function``, an implicit function (``tack.dataset.implicit``: ``Plane``,
    ``Sphere``, ``Box``, ...), at the geometry's values: a field in the
    geometry's space and type, negative inside the function's region. On an
    ``H1`` geometry that is one value per point. The function is in world
    positions, taken relative to the dataset's origin in double precision."""
    positions = materialize(data.geometry.values)
    n = size_of(data.geometry.values)
    dtype = arrays.dtype_of(data.geometry.values)
    out = tack.empty((n,), dtype)
    # A function holding arrays (Planes) gives a copy at the data's precision and
    # on its device; the others hold Python numbers, which take it anyway.
    from tack.dataset.implicit import relative_to

    function = relative_to(function, data.origin)    # the geometry's values are relative
    matching = getattr(type(function), "_matching", None)
    if matching is not None:
        function = matching(function, dtype, data.device)
    if n:
        _implicit_values(positions, function, out)
    return Field(data.geometry.space, out)

cell_ids 

cell_ids(data)

Each cell's id, as values on cells: an implicit array, nothing stored.

Stored as a field (data.fields["source"] = cell_ids(data)), it reaches every filter's output as the id of the input cell each output cell came from -- kept, a piece of, or a face of -- with no option on the filter.

Source code in packages/tack-vis/src/tack/dataset/algorithms/cell_ids.py 
11
12
13
14
15
16
17
18
@on_inputs_device
def cell_ids(data):
    """Each cell's id, as values on cells: an implicit array, nothing stored.

    Stored as a field (``data.fields["source"] = cell_ids(data)``), it reaches
    every filter's output as the id of the input cell each output cell came
    from -- kept, a piece of, or a face of -- with no option on the filter."""
    return Field(Values(data, "cells"), CountingArray(data.num_cells, dtype=data.id_dtype))

point_ids 

point_ids(data)

Each point's id, as values on points: an implicit array, nothing stored.

Stored as a field, it reaches every filter's output as the id of the input point each kept point is (cell_ids does the same for cells).

Source code in packages/tack-vis/src/tack/dataset/algorithms/point_ids.py 
11
12
13
14
15
16
17
@on_inputs_device
def point_ids(data):
    """Each point's id, as values on points: an implicit array, nothing stored.

    Stored as a field, it reaches every filter's output as the id of the input
    point each kept point is (``cell_ids`` does the same for cells)."""
    return Field(Values(data, "points"), CountingArray(data.num_points, dtype=data.id_dtype))

Fields from fields

New fields computed from a dataset's fields.

Algorithm Gives Takes VTK Viskores
gradients the gradient of a scalar or vector field, at cells or points dataset vtkGradientFilter Gradient
flow_quantities divergence, vorticity and Q-criterion from a velocity gradient dataset vtkGradientFilter Gradient
project a field moved into another space (cells to points, points to cells, ...) dataset vtkCellDataToPointData, vtkPointDataToCellData PointAverage, CellAverage
values_at_centers a field at each cell's parametric center dataset — —
interpolate a @tack.func of position at a space's nodes dataset — —
magnitude each vector's length dataset — VectorMagnitude
dot the dot product of two vector fields dataset — DotProduct
cross the cross product of two 3-vector fields dataset — CrossProduct
composite scalar fields as the components of a vector field dataset — CompositeVectors
log_values the logarithm of a field, in any base dataset — LogValues
convert a field, or the geometry, in another unit dataset — —

gradients 

gradients(data, field, at='cells')

The gradient of field -- H1, L2 or Constant, scalar or vector -- at each cell's parametric center (at="cells") or at the points (at="points": each cell's gradient at the point, averaged over the cells around it, as VTK's vtkGradientFilter and Viskores' point gradients do).

A scalar gives a 3-vector per entity; a vector of k components 3k values, row by row -- d(u_i)/d(x_j) at 3 * i + j, the layout of VTK's and Viskores' gradient tensors. The field's basis gives its derivative in parametric coordinates and the geometry's Jacobian turns that into a world gradient, so the two need not share a space. A 2D cell's gradient lies in its plane (VTK's and Viskores' in-plane gradient), lines and vertices get zero, and a cell of no volume or area NaN (its Jacobian is singular). flow_quantities derives divergence, vorticity and the Q-criterion from a vector field's gradient.

On a polyhedral topology there is no reference element. A cell's gradient of point data is then Green-Gauss over its boundary: the faces walked outward and fanned from their points' means (the triangles cell_geometry measures the volume by), each triangle's area vector times the mean of its values, over the volume. That is exact for a linear field and needs no interpolant. Cell data, constant in each cell, has a zero gradient there, as on shapes. A polygon's is the same in its plane, over its edges; polylines get zero. At the points, the cells' gradients are averaged over the cells around each point.

Source code in packages/tack-vis/src/tack/dataset/algorithms/gradients.py 
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
@on_inputs_device
def gradients(data, field, at="cells"):
    """The gradient of ``field`` -- ``H1``, ``L2`` or ``Constant``, scalar or vector --
    at each cell's parametric center (``at="cells"``) or at the points
    (``at="points"``: each cell's gradient at the point, averaged over the cells
    around it, as VTK's vtkGradientFilter and Viskores' point gradients do).

    A scalar gives a 3-vector per entity; a vector of ``k`` components ``3k``
    values, row by row -- ``d(u_i)/d(x_j)`` at ``3 * i + j``, the layout of
    VTK's and Viskores' gradient tensors. The field's basis gives its
    derivative in parametric coordinates and the geometry's Jacobian turns that
    into a world gradient, so the two need not share a space. A 2D cell's
    gradient lies in its plane (VTK's and Viskores' in-plane gradient), lines
    and vertices get zero, and a cell of no volume or area NaN (its Jacobian is
    singular). ``flow_quantities`` derives divergence, vorticity and
    the Q-criterion from a vector field's gradient.

    On a polyhedral topology there is no reference element. A cell's gradient
    of point data is then Green-Gauss over its boundary: the faces walked
    outward and fanned from their points' means (the triangles ``cell_geometry``
    measures the volume by), each triangle's area vector times the mean of its
    values, over the volume. That is exact for a linear field and needs no
    interpolant. Cell data, constant in each cell, has a zero gradient there,
    as on shapes. A polygon's is the same in its plane, over its edges; polylines
    get zero. At the points, the cells' gradients are averaged over the cells
    around each point."""
    field = _with_basis(data, _field(data, field))
    if (not getattr(data.topology, "reference_cells", True)
            and not isinstance(field.space, (H1, Constant))):
        raise NotImplementedError(f"gradients of a field on {field.space!r} on a mesh "
                                  "without reference elements: point or cell data only")
    if at not in ("cells", "points"):
        raise ValueError(f"gradients are at 'cells' or 'points', not {at!r}")
    space = Values(data, at)
    # Per unit length; a difference of degrees Celsius is in kelvin.
    unit = units.derived(units.delta(field.unit), (data.unit, -1))
    width = width_of(field.values)
    if not width:
        return Field(space, _scalar_gradients(data, field, at), unit)
    n = data.num_cells if at == "cells" else data.num_points
    out = tack.empty((n,), tack.vec(3 * width, data.dtype))
    values = materialize(field.values)
    for k in range(width):
        part = tack.empty((size_of(field.values),), arrays.dtype_of(field.values))
        if part.shape[0]:
            _component(values, k, part)
        rows = _scalar_gradients(data, Field(field.space, part), at)
        if n:
            _place_rows(rows, k, out)
    return Field(space, out, unit)

flow_quantities 

flow_quantities(data, gradient)

Divergence, vorticity and the Q-criterion from the gradient of a 3-vector field (gradients of it: 9 values per entity), in the gradient's space, as VTK's vtkGradientFilter and Viskores' Gradient compute them, as FlowQuantities: divergence, vorticity and q_criterion.

Source code in packages/tack-vis/src/tack/dataset/algorithms/flow_quantities.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
@on_inputs_device
def flow_quantities(data, gradient):
    """Divergence, vorticity and the Q-criterion from the gradient of a 3-vector field
    (``gradients`` of it: 9 values per entity), in the gradient's space, as VTK's
    vtkGradientFilter and Viskores' Gradient compute them, as ``FlowQuantities``:
    ``divergence``, ``vorticity`` and ``q_criterion``."""
    gradient = _field(data, gradient)
    if width_of(gradient.values) != 9:
        raise TypeError("flow_quantities takes the gradient of a 3-vector field (9 values)")
    n = size_of(gradient.values)
    dtype = arrays.dtype_of(gradient.values)
    divergence = tack.empty((n,), dtype)
    vorticity = tack.empty((n,), tack.vec(3, dtype))
    q = tack.empty((n,), dtype)
    if n:
        _flow_quantities(materialize(gradient.values), divergence, vorticity, q)
    unit = gradient.unit
    return FlowQuantities(Field(gradient.space, divergence, unit),
                          Field(gradient.space, vorticity, unit),
                          Field(gradient.space, q, units.derived((unit, 2))))

FlowQuantities

Bases: NamedTuple

flow_quantities' fields, each in the gradient's space.

project 

project(data, field, space)

field (a Field or a name) in space, a space of the same topology: a new field. Values on points or cells count as order-1 H1 and Constant, both as source and as target.

Into H1, each node takes the mean of the source evaluated there by the cells around it: the source's own value wherever it is continuous (so a continuous field of lower order comes across exactly), and an average where cells disagree -- a cell field's, as vtkCellDataToPointData takes it, or a DG field's corner values. Into Constant (or order-0 L2), each cell takes the mean of its corner values, as vtkPointDataToCellData does. Into order-1 L2, each cell takes the source's values at its corners, exactly whatever the source's order is up to 1. Into SideTraces, each face side the source's values from its cell (traces). Sums run in a fixed order, so results do not vary from run to run.

Values on faces or edges have no basis to project, and are refused as a source or target. On polyhedra, polygons and polylines -- no reference element -- only cell values and point values project, one into the other, each cell averaging over its distinct points. The values must be floating point, scalars or vectors.

Source code in packages/tack-vis/src/tack/dataset/algorithms/project.py 
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
@on_inputs_device
def project(data, field, space):
    """``field`` (a ``Field`` or a name) in ``space``, a space of the same topology:
    a new field. Values on points or cells count as order-1 ``H1`` and
    ``Constant``, both as source and as target.

    Into ``H1``, each node takes the mean of the source evaluated there by the
    cells around it: the source's own value wherever it is continuous (so a
    continuous field of lower order comes across exactly), and an average where
    cells disagree -- a cell field's, as vtkCellDataToPointData takes it, or a DG
    field's corner values. Into ``Constant`` (or order-0 ``L2``), each cell takes
    the mean of its corner values, as vtkPointDataToCellData does. Into order-1
    ``L2``, each cell takes the source's values at its corners, exactly whatever
    the source's order is up to 1. Into ``SideTraces``, each face side the
    source's values from its cell (``traces``). Sums run in a fixed order, so
    results do not vary from run to run.

    Values on faces or edges have no basis to project, and are refused as a
    source or target. On polyhedra, polygons and polylines -- no reference
    element -- only cell values and point values project, one into the other,
    each cell averaging over its distinct points. The values must be floating
    point, scalars or vectors.
    """
    field = _with_basis(data, _field(data, field))
    source = field.space
    if space.topology is not data.topology:
        raise ValueError("project moves a field between spaces of one topology")
    if isinstance(space, Values) and space.on in ("points", "cells"):
        projected = project(data, field, H1(data) if space.on == "points" else Constant(data))
        return Field(space, projected.values, field.unit)
    for which, s in (("from", source), ("into", space)):
        if not s.interpolated:
            raise TypeError(f"project works {which} a space with a basis, or values on the "
                            f"points or cells; {s!r} has no basis to project")
    if isinstance(source, SideTraces):
        raise TypeError("a field of traces has two values on each face; project the field "
                        "they were taken from")
    values = field.values
    if space is source:
        return Field(space, _take(values, tack.arange(space.size, dtype=data.id_dtype)),
                     field.unit)
    if dtype_of(values) not in (tack.f32, tack.f64):
        raise TypeError(f"project averages and interpolates floating-point values, not "
                        f"{dtype_of(values).name}")
    if isinstance(space, SideTraces):
        return traces(data, field)
    if not getattr(data.topology, "reference_cells", True):
        if isinstance(space, H1) and space.order == 1 and isinstance(source, Constant):
            return Field(space, _onto_points(data, field, values), field.unit)
        if isinstance(space, Constant) and isinstance(source, H1) and source.order == 1:
            out = _like(values, data.num_cells)
            for_each(_distinct_point_means, data, "cells", materialize(values), out)
            return Field(space, out, field.unit)
        reference_cells(data.topology, f"project from {source!r} to {space!r}")
    if isinstance(space, H1):
        if space.order == 1:
            return Field(space, _onto_points(data, field, values), field.unit)
        return Field(space, _onto_nodes(data, field, space, values), field.unit)
    out = _like(values, space.size)
    if isinstance(space, Constant) or (isinstance(space, L2) and not space.varies
                                       and space.order == 0):
        for group in data.launch_groups("cells", [field]):
            if group.count:
                _cell_means(data.domain_view("cells", group), field.view(group), out)
        return Field(space, out, field.unit)
    target = Field(space, out)
    for group in data.launch_groups("cells", [field, target]):
        if group.count:
            _cell_nodes(data.domain_view("cells", group), field.view(group), target.view(group),
                        space.offsets, out)
    return Field(space, out, field.unit)

values_at_centers 

values_at_centers(data, field)

field evaluated at each cell's parametric center, through its basis.

On a polyhedron, polygon or polyline, which have no reference element, the center is the cell point (see cell_centers) and point data (H1) takes its value there, the mean of the cell's distinct points' values: the value of the cell's interpolant, linear on tetrahedra fanned from that point. Cell data (Constant) is its own value.

Source code in packages/tack-vis/src/tack/dataset/algorithms/values_at_centers.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
@on_inputs_device
def values_at_centers(data, field):
    """``field`` evaluated at each cell's parametric center, through its basis.

    On a polyhedron, polygon or polyline, which have no reference element, the
    center is the cell point (see ``cell_centers``) and point data (``H1``)
    takes its value there, the mean of the cell's distinct points' values: the
    value of the cell's interpolant, linear on tetrahedra fanned from that
    point. Cell data (``Constant``) is its own value."""
    field = _with_basis(data, _field(data, field))
    out = _like(field.values, data.num_cells)
    if getattr(data.topology, "reference_cells", True):
        for_each(_at_centers, data, "cells", field, out)
    elif isinstance(field.space, H1) and field.space.order == 1:
        from tack.dataset.algorithms.project import _distinct_point_means

        for_each(_distinct_point_means, data, "cells", materialize(field.values), out)
    elif isinstance(field.space, Constant):
        tack.copyto(out, materialize(field.values))
    else:
        raise NotImplementedError(f"values_at_centers of {field.space!r} on polyhedra, "
                                  "polygons or polylines: only point data (H1, order 1) and "
                                  "cell data (Constant) have values there")
    return Field(Values(data, "cells"), out, field.unit)

interpolate 

interpolate(data, function, space, dtype=None, unit=None)

A Field in space whose values are function at each of the space's nodes: function is a @tack.func of a world position (a 3-vector of the geometry's type) returning a scalar or a vector, and dtype its result's type -- data.dtype by default, tack.vec3() for 3-vectors.

The nodes are where the geometry puts them, evaluated as needed and never stored: H1's points (and at order 2 its edge middles, quad face centers and cell centers), L2's cell corners (order 1) or parametric centers (order 0), Constant's and Values(data, "cells")' parametric centers, Values(data, "points")' points. A node shared by several cells is evaluated once, by the first cell that has it, so a discontinuous (L2) geometry gives a shared node the first cell's position. A point no cell uses takes function at its own position when the geometry is H1, otherwise zero.

A function of position is what the space can represent exactly when it is a polynomial of the space's order on straight cells; otherwise this is the space's nodal interpolant. On polyhedra, polygons and polylines -- no reference element -- only point values (H1 of order 1 and Values(data, "points")) of an H1 geometry interpolate.

unit is the values' unit (Field.unit): a function's code does not say it, so by default it is unknown.

Source code in packages/tack-vis/src/tack/dataset/algorithms/interpolate.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
@on_inputs_device
def interpolate(data, function, space, dtype=None, unit=None):
    """A ``Field`` in ``space`` whose values are ``function`` at each of the space's
    nodes: ``function`` is a ``@tack.func`` of a world position (a 3-vector of
    the geometry's type) returning a scalar or a vector, and ``dtype`` its
    result's type -- ``data.dtype`` by default, ``tack.vec3()`` for 3-vectors.

    The nodes are where the geometry puts them, evaluated as needed and never
    stored: ``H1``'s points (and at order 2 its edge middles, quad face centers
    and cell centers), ``L2``'s cell corners (order 1) or parametric centers
    (order 0), ``Constant``'s and ``Values(data, "cells")``' parametric
    centers, ``Values(data, "points")``' points. A node shared by several
    cells is evaluated once, by the first cell that has it, so a discontinuous
    (``L2``) geometry gives a shared node the first cell's position. A point no
    cell uses takes ``function`` at its own position when the geometry is
    ``H1``, otherwise zero.

    A function of position is what the space can represent exactly when it is
    a polynomial of the space's order on straight cells; otherwise this is the
    space's nodal interpolant. On polyhedra, polygons and polylines -- no
    reference element -- only point values (``H1`` of order 1 and
    ``Values(data, "points")``) of an ``H1`` geometry interpolate.

    ``unit`` is the values' unit (``Field.unit``): a function's code does not say
    it, so by default it is unknown."""
    if not isinstance(function, Func):
        raise TypeError(f"interpolate takes a @tack.func of position, not {function!r}")
    if space.topology is not data.topology:
        raise ValueError("interpolate takes a space of the dataset's own topology")
    if isinstance(space, Values) and space.on not in ("points", "cells"):
        raise TypeError(f"interpolate gives values at points or cells, not on {space.on}")
    if not isinstance(space, (H1, L2, Constant, Values)):
        raise TypeError(f"interpolate takes H1, L2, Constant or Values, not {space!r}")
    fn = _Of(function, data.origin)
    values = tack.zeros((space.size,), dtype or data.dtype)
    out = Field(space, values, unit)
    shared = isinstance(space, H1) or (isinstance(space, Values) and space.on == "points")
    if not (shared and not isinstance(data.geometry.space, H1)):
        # Every node once, from the space's node views.
        for_each(_at_space_nodes, data, space, fn, values)
        return out
    # An L2 geometry gives a shared node no position of its own: the first cell
    # that has it decides, through its corners.
    reference_cells(data.topology, f"interpolate into {space!r}")
    if isinstance(space, Values):
        raise ValueError("an L2 geometry has positions per cell corner, not per point")
    starts, writers = _writers(data, space, out)
    for group in data.launch_groups("cells", [out]):
        if not group.count:
            continue
        cells = data.domain_view("cells", group)
        if space.order == 1:
            _at_corners(cells, fn, writers, starts[group.root], values)
        else:
            _at_nodes(cells, out.view(group), fn, writers, starts[group.root], 0, values)
    return out

magnitude 

magnitude(data, field)

Each vector's length (Viskores' VectorMagnitude): a scalar field in the space and unit of field, a vector field or its name. The values convert to f32 first, or stay f64, so integer vectors give floating-point lengths.

Source code in packages/tack-vis/src/tack/dataset/algorithms/magnitude.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
@on_inputs_device
def magnitude(data, field):
    """Each vector's length (Viskores' VectorMagnitude): a scalar field in the
    space and unit of ``field``, a vector field or its name. The values convert
    to f32 first, or stay f64, so integer vectors give floating-point lengths."""
    field = _field(data, field)
    if not width_of(field.values):
        raise TypeError("magnitude takes a vector field")
    n = size_of(field.values)
    dtype = _floating_dtype(field.values)
    out = _scalars(n, dtype)
    if n:
        _magnitudes(_converted(field.values, dtype), out)
    return Field(field.space, out, field.unit)

dot 

dot(data, a, b)

a . b at every value of two vector fields a and b (fields or names) of one width, dtype and space: a scalar field in that space and dtype (tack.arrays.Cast converts one first). Its unit is the product of theirs, by its SI name where it has one ("N" times "m" is "J").

Source code in packages/tack-vis/src/tack/dataset/algorithms/dot.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
@on_inputs_device
def dot(data, a, b):
    """``a . b`` at every value of two vector fields ``a`` and ``b`` (fields or
    names) of one width, dtype and space: a scalar field in that space and
    dtype (``tack.arrays.Cast`` converts one first). Its unit is the product of
    theirs, by its SI name where it has one (``"N"`` times ``"m"`` is ``"J"``)."""
    a, b = _field(data, a), _field(data, b)
    space, n = _same_space(a, b)
    if not width_of(a.values) or width_of(a.values) != width_of(b.values):
        raise TypeError("dot takes two vector fields of one width")
    _same_dtype("dot", a, b)
    out = _scalars(n, arrays.dtype_of(a.values))
    if n:
        _dots(materialize(a.values), materialize(b.values), out)
    return Field(space, out, units.named(units.derived(a.unit, b.unit)))

cross 

cross(data, a, b)

a x b at every value of two 3-vector fields a and b (fields or names) of one dtype and space: a 3-vector field in that space and dtype (tack.arrays.Cast converts one first). Its unit is the product of theirs, by its SI name where it has one.

Source code in packages/tack-vis/src/tack/dataset/algorithms/cross.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
@on_inputs_device
def cross(data, a, b):
    """``a x b`` at every value of two 3-vector fields ``a`` and ``b`` (fields
    or names) of one dtype and space: a 3-vector field in that space and dtype
    (``tack.arrays.Cast`` converts one first). Its unit is the product of
    theirs, by its SI name where it has one."""
    a, b = _field(data, a), _field(data, b)
    space, n = _same_space(a, b)
    if width_of(a.values) != 3 or width_of(b.values) != 3:
        raise TypeError("cross takes two 3-vector fields")
    _same_dtype("cross", a, b)
    out = _vectors(n, arrays.dtype_of(a.values))
    if n:
        _crosses(materialize(a.values), materialize(b.values), out)
    return Field(space, out, units.named(units.derived(a.unit, b.unit)))

composite 

composite(data, *fields)

Scalar fields of one space as the components of one vector field (Viskores' CompositeVectors), in the order given. The fields must be in one unit, the vector's (an unknown one is not checked).

Source code in packages/tack-vis/src/tack/dataset/algorithms/composite.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
@on_inputs_device
def composite(data, *fields):
    """Scalar fields of one space as the components of one vector field (Viskores'
    CompositeVectors), in the order given. The fields must be in one unit, the
    vector's (an unknown one is not checked)."""
    fields = [_field(data, f) for f in fields]
    if len(fields) < 2:
        raise ValueError("composite takes at least two scalar fields")
    space, n = _same_space(*fields)
    if any(width_of(f.values) for f in fields):
        raise TypeError("composite takes scalar fields")
    _same_dtype("composite", *fields)
    unit = units.common("composite's fields", [f.unit for f in fields])
    out = _vectors(n, arrays.dtype_of(fields[0].values), len(fields))
    for k, f in enumerate(fields):
        if n:
            _set_component(materialize(f.values), k, out)
    return Field(space, out, unit)

log_values 

log_values(data, field, base=math.e, min_value=None)

The logarithm of a scalar field in base (e, 2, 10, or any), values below min_value -- by default the smallest positive normal float, as in Viskores' LogValues -- taken as min_value. The field must be a pure number ("dimensionless") or of unknown unit: a logarithm of metres is refused.

Source code in packages/tack-vis/src/tack/dataset/algorithms/log_values.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
@on_inputs_device
def log_values(data, field, base=math.e, min_value=None):
    """The logarithm of a scalar field in ``base`` (e, 2, 10, or any), values below
    ``min_value`` -- by default the smallest positive normal float, as in
    Viskores' LogValues -- taken as ``min_value``. The field must be a pure number
    (``"dimensionless"``) or of unknown unit: a logarithm of metres is refused."""
    field = _field(data, field)
    if width_of(field.values):
        raise TypeError("log_values takes a scalar field")
    unit = units.dimensionless_or_none("log_values", field.unit)
    if not base > 0 or base == 1:
        raise ValueError(f"no logarithm in base {base}")
    dtype = _floating_dtype(field.values)
    if min_value is None:
        min_value = float(np.finfo(dtype.numpy_dtype).tiny)
    n = size_of(field.values)
    out = _scalars(n, dtype)
    if n:
        _logs(_converted(field.values, dtype), 1.0 / math.log(base), float(min_value), out)
    return Field(field.space, out, unit)

convert 

convert(data, unit, field=None)

data with field (a name) in unit, or with the geometry in unit when field is None: a new dataset, the values scaled (and offset, as from degrees Celsius to kelvin) on the device, the topology, sets, marks and other fields kept.

The field's unit, or the spatial unit, must be known (Field.unit, DataSet.unit) and of unit's dimensions. Converting the geometry scales the origin too, in double precision on the host, so survey coordinates in US survey feet become metres without losing digits; an image stays an image and a rectilinear grid keeps its axes. Integer fields are refused unless the conversion is a relabelling ("1000*m" to "km").

Source code in packages/tack-vis/src/tack/dataset/algorithms/convert.py 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
@on_inputs_device
def convert(data, unit, field=None):
    """``data`` with ``field`` (a name) in ``unit``, or with the geometry in ``unit``
    when ``field`` is ``None``: a new dataset, the values scaled (and offset,
    as from degrees Celsius to kelvin) on the device, the topology, sets, marks
    and other fields kept.

    The field's unit, or the spatial unit, must be known (``Field.unit``,
    ``DataSet.unit``) and of ``unit``'s dimensions. Converting the geometry
    scales the origin too, in double precision on the host, so survey
    coordinates in US survey feet become metres without losing digits; an
    image stays an image and a rectilinear grid keeps its axes. Integer fields
    are refused unless the conversion is a relabelling (``"1000*m"`` to
    ``"km"``)."""
    from tack.dataset.algorithms._moving import _moved

    target = units.canonical(unit)
    if field is None:
        factor, offset = units.conversion(data.unit, target)
        if offset:
            raise ValueError(f"positions in {data.unit!r} are not lengths with a zero")
        values = data.geometry.values
        if isinstance(values, UniformCoordinates):
            moved = UniformCoordinates(values.dims, np.asarray(values._origin) * factor,
                                       np.asarray(values._spacing) * factor,
                                       values._matrix if values.ORIENTED else None,
                                       dtype=values.dtype, device=data.device)
        elif isinstance(values, CartesianProduct):
            axes = [arrays.to_host(axis).astype(np.float64).reshape(-1) * factor
                    for axis in (values.xs, values.ys, values.zs)]
            moved = CartesianProduct(*axes, dtype=values.dtype, device=data.device)
        else:
            moved = _times(values, factor, 0.0)
        out = _moved(data, moved)
        out._geometry = Field(out.geometry.space, out.geometry.values, target)
        out.origin = np.asarray(data.origin) * factor
        return out
    if not isinstance(field, str):
        raise TypeError(f"convert takes a field's name, not {field!r}")
    if field not in data.fields:
        raise KeyError(f"no field {field!r}: the dataset has "
                       f"{', '.join(map(repr, data.fields)) or 'none'}")
    given = data.fields[field]
    factor, offset = units.conversion(given.unit, target)
    if factor == 1.0 and offset == 0.0:
        values = given.values
    elif arrays.dtype_of(given.values).numpy_dtype.kind != "f":
        raise TypeError(f"field {field!r} holds integers, which cannot hold its values in "
                        f"{target!r}: give it floating-point values first")
    else:
        values = _times(given.values, factor, offset)
    fields = dict(data.fields)
    fields[field] = Field(given.space, values, target)
    from tack.dataset.ghosts import keep_marks

    return keep_marks(data, DataSet(data.topology, data.geometry, fields=fields,
                                    sets=dict(data.sets)))

Finite volumes

Values on faces and the cells around them.

Algorithm Gives Takes VTK Viskores
jump a field's difference across each face dataset — —
upwind_flux a field's upwind flux through each face dataset — —
net_flux each cell's outward sum of face fluxes dataset — —
divergence each cell's finite-volume divergence of face fluxes dataset — —
perot cell vectors from face-normal components (C-grid models) dataset — —

jump 

jump(data, field)

A field's difference across each face at the face's center, side 1 minus side 0 (zero on the boundary): a field on faces. Each side is the field as its own cell has it -- traces -- so a DG field's jumps are its own, not its cells' averages; a continuous field's are zero. Cell data (Constant, values on cells) is read from each side's cell directly, on any topology, polyhedral included.

Source code in packages/tack-vis/src/tack/dataset/algorithms/jump.py 
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
@on_inputs_device
def jump(data, field):
    """A field's difference across each face at the face's center, side 1 minus side
    0 (zero on the boundary): a field on faces. Each side is the field as its own
    cell has it -- ``traces`` -- so a DG field's jumps are its own, not its cells'
    averages; a continuous field's are zero. Cell data (``Constant``, values on
    cells) is read from each side's cell directly, on any topology, polyhedral
    included."""
    field = _field(data, field)
    out = _like(field.values, data.topology.faces().num_faces)
    if _on_cells(data, field):
        for_each(_cell_jump, data, "faces", materialize(field.values), out)
    else:
        for_each(_jump, data, "faces", traces(data, field), out)
    return Field(Values(data, "faces"), out, units.delta(field.unit))

upwind_flux 

upwind_flux(data, field, velocity)

The flux of field carried by a constant velocity through each face, out of side 0: (v . n) * area * u, u taken from the upwind side at the face's center -- the side the flow leaves -- as a DG advection scheme does. A boundary face takes its one side's value whichever way the flow goes. divergence of the result is each cell's net outflow.

Source code in packages/tack-vis/src/tack/dataset/algorithms/upwind_flux.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
@on_inputs_device
def upwind_flux(data, field, velocity):
    """The flux of ``field`` carried by a constant ``velocity`` through each face, out
    of side 0: ``(v . n) * area * u``, ``u`` taken from the upwind side at the face's
    center -- the side the flow leaves -- as a DG advection scheme does. A boundary
    face takes its one side's value whichever way the flow goes. ``divergence`` of
    the result is each cell's net outflow."""
    field = _field(data, field)
    normals, areas = face_geometry(data)
    out = _like(field.values, data.topology.faces().num_faces)
    velocity, speed = units.number_and_unit(velocity)
    vx, vy, vz = (float(v) for v in velocity)
    if _on_cells(data, field):
        for_each(_cell_upwind, data, "faces", materialize(field.values), normals.values,
                 areas.values, vx, vy, vz, out)
    else:
        for_each(_upwind, data, "faces", traces(data, field), normals.values, areas.values,
                 vx, vy, vz, out)
    # (v . n) area u: a velocity given as a bare number leaves it unknown.
    area = units.derived((data.unit, data.topology.dimension - 1))
    return Field(Values(data, "faces", oriented=True), out,
                 units.named(units.derived(field.unit, speed, area)))

net_flux 

net_flux(data, flux)

Each cell's outward sum of flux, values on faces measured out of side 0 (Values(data, "faces", oriented=True), or plain face values taken so): a face's value counts plus for its side-0 cell and minus for its side-1 cell.

Given each face's integrated flux (upwind_flux, or a normal component times the area), it is the net outflow of each cell, V times the divergence; divergence divides by the cells' measure. Nothing is weighted: the values are summed as they are, so face area vectors sum to zero on a closed cell. Cells of a lower dimension than the mesh's have no faces: zero.

Source code in packages/tack-vis/src/tack/dataset/algorithms/net_flux.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
@on_inputs_device
def net_flux(data, flux):
    """Each cell's outward sum of ``flux``, values on faces measured out of side 0
    (``Values(data, "faces", oriented=True)``, or plain face values taken so): a
    face's value counts plus for its side-0 cell and minus for its side-1 cell.

    Given each face's integrated flux (``upwind_flux``, or a normal component times
    the area), it is the net outflow of each cell, ``V`` times the divergence;
    ``divergence`` divides by the cells' measure. Nothing is weighted: the
    values are summed as they are, so face area vectors sum to zero on a closed
    cell. Cells of a lower dimension than the mesh's have no faces: zero."""
    flux = _field(data, flux)
    if not _on_faces(data, flux):
        raise TypeError("net_flux sums a field on faces")
    data.topology.faces()                          # the incidence the cell views need
    out = _like(flux.values, data.num_cells)
    out.fill(0)
    # Cells of a lower dimension (a surface's triangles among tetrahedra) have
    # no faces of the mesh's, so no flux: zero.
    values = materialize(flux.values)
    for group in data.launch_groups("cells", []):
        if group.count and group.shape.DIMENSION == data.topology.dimension:
            _outward_sums(data.domain_view("cells", group), values, out)
    return Field(Values(data, "cells"), out, flux.unit)

divergence 

divergence(data, flux)

Each cell's finite-volume divergence of flux, the integrated flux through each face out of its side 0 (upwind_flux, or a normal component times the area): net_flux over the cell's measure (cell_geometry), so for a linear field it is the field's divergence. A cell of no measure gives no finite value, as perot and gradients give. flow_quantities gives the finite-element divergence instead, the trace of a field's gradient.

Source code in packages/tack-vis/src/tack/dataset/algorithms/divergence.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
@on_inputs_device
def divergence(data, flux):
    """Each cell's finite-volume divergence of ``flux``, the integrated flux through
    each face out of its side 0 (``upwind_flux``, or a normal component times the
    area): ``net_flux`` over the cell's measure (``cell_geometry``), so for a
    linear field it is the field's divergence. A cell of no measure gives no
    finite value, as ``perot`` and ``gradients`` give. ``flow_quantities`` gives
    the finite-element divergence instead, the trace of a field's gradient."""
    net = net_flux(data, flux)
    measure = cell_geometry(data).measure
    out = _like(net.values, data.num_cells)
    if data.num_cells:
        _per_measure(net.values, measure.values, out)
    return Field(Values(data, "cells"), out, units.derived(net.unit, (measure.unit, -1)))

perot 

perot(data, flux)

A vector per cell from each face's normal component -- Perot's reconstruction, u_c = (1/V) sum_f A_f u_f (r_f - r_c), as C-grid models (MPAS, ICON) need to show their face-normal velocity at cells.

flux is a scalar field on faces, the component along each face's normal out of its side 0; the side each cell is on gives the sign. Face centers, areas and cell centroids and volumes are this module's, so a uniform field's normal components give it back exactly on cells with planar faces. First order, as Perot's method is. A cell of no volume gives NaN, and one of a lower dimension than the mesh's (no faces of it) zero.

Source code in packages/tack-vis/src/tack/dataset/algorithms/perot.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
@on_inputs_device
def perot(data, flux):
    """A vector per cell from each face's normal component -- Perot's reconstruction,
    ``u_c = (1/V) sum_f A_f u_f (r_f - r_c)``, as C-grid models (MPAS, ICON) need to
    show their face-normal velocity at cells.

    ``flux`` is a scalar field on faces, the component along each face's normal out
    of its side 0; the side each cell is on gives the sign. Face centers, areas
    and cell centroids and volumes are this module's, so a uniform field's normal
    components give it back exactly on cells with planar faces. First order, as
    Perot's method is. A cell of no volume gives NaN, and one of a lower dimension
    than the mesh's (no faces of it) zero."""
    flux = _field(data, flux)
    if not _on_faces(data, flux):
        raise TypeError("perot reconstructs from a field on faces")
    _, areas = face_geometry(data)
    centers = face_centers(data)
    volumes, centroids = cell_geometry(data)
    out = _vectors(data.num_cells, data.dtype)
    out.fill(0)                                    # cells of a lower dimension: no faces
    for group in data.launch_groups("cells", []):
        if group.count and group.shape.DIMENSION == data.topology.dimension:
            _perot(data.domain_view("cells", group), materialize(flux.values), areas.values,
                   centers.values, centroids.values, volumes.values, out)
    # Areas times offsets over volumes cancel: the components' unit.
    return Field(Values(data, "cells"), out, flux.unit)

Summaries

Python numbers about a field, over every piece and rank.

Algorithm Gives Takes VTK Viskores
statistics count, extrema, mean, moments, variance, skewness, kurtosis; weighted dataset, pieces, AMR, octree vtkDescriptiveStatistics Statistics
entropy Shannon entropy of a field's histogram dataset, pieces, AMR, octree — Entropy

statistics 

statistics(data, field, weights=None, level=None)

Descriptive statistics of a scalar field's values, as Viskores' Statistics filter gives them: a Statistics of n, min, max, sum, mean, the central moment sums m2, m3, m4, sample_variance (m2 / (n - 1)), population_variance (m2 / n), their standard deviations, skewness (sqrt(n) m3 / m2^1.5) and kurtosis (n m4 / m2^2, not the excess), the last two 0 for a constant field. Moments are taken about the mean, in a second pass. Ghosts are left out: cells and points data.ghosts marks a duplicate (another piece owns them) or hidden.

Over a PartitionedDataSet, field names a field every piece has: each piece's moments, over what it owns, are taken on its own device and combined pairwise, in the order of the pieces over every rank of data.world -- the same result however the pieces are shared out. Where every shared cell and point has one owner, the result is the whole's up to rounding. Over an AMRDataSet, the same over every box of every level: refined cells and points are left out, so each region of space counts once, at its finest level. Over an OctreeDataSet, its leaves (or the nodes of one level), every rank's.

weights weighs each value: "volume", each cell's measure (a cell field's), or a scalar field on the same entities, by name (a density, for mass-weighted moments). Then mean is sum(w x) / W, sum is sum(w x), the moments are sum(w (x - mean)^k), weight is W, the sum of the weights, the population variance m2 / W, the sample variance m2 / (W - sum(w^2) / W) (reliability weights, NumPy's aweights), and skewness and kurtosis take W for n. Cells of different sizes -- an AMR hierarchy's levels, an octree's leaves -- then count by their size.

Source code in packages/tack-vis/src/tack/dataset/algorithms/statistics.py 
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
@on_inputs_device
def statistics(data, field, weights=None, level=None):
    """Descriptive statistics of a scalar field's values, as Viskores' Statistics
    filter gives them: a ``Statistics`` of ``n``, ``min``, ``max``, ``sum``, ``mean``, the
    central moment sums ``m2``, ``m3``, ``m4``, ``sample_variance``
    (``m2 / (n - 1)``), ``population_variance`` (``m2 / n``), their standard
    deviations, ``skewness`` (``sqrt(n) m3 / m2^1.5``) and ``kurtosis``
    (``n m4 / m2^2``, not the excess), the last two 0 for a constant field.
    Moments are taken about the mean, in a second pass. Ghosts are left out:
    cells and points ``data.ghosts`` marks a duplicate (another piece owns
    them) or hidden.

    Over a ``PartitionedDataSet``, ``field`` names a field every piece has: each
    piece's moments, over what it owns, are taken on its own device and
    combined pairwise, in the order of the pieces over every rank of
    ``data.world`` -- the same result however the pieces are shared out. Where
    every shared cell and point has one owner, the result is the whole's up to
    rounding. Over an ``AMRDataSet``, the same over every box of every level:
    refined cells and points are left out, so each region of space counts
    once, at its finest level. Over an ``OctreeDataSet``, its leaves (or the
    nodes of one ``level``), every rank's.

    ``weights`` weighs each value: ``"volume"``, each cell's measure (a cell
    field's), or a scalar field on the same entities, by name (a density, for
    mass-weighted moments). Then ``mean`` is ``sum(w x) / W``, ``sum`` is ``sum(w
    x)``, the moments are ``sum(w (x - mean)^k)``, ``weight`` is ``W``, the sum of
    the weights, the population variance ``m2 / W``, the sample variance
    ``m2 / (W - sum(w^2) / W)`` (reliability weights, NumPy's ``aweights``), and
    skewness and kurtosis take ``W`` for ``n``. Cells of different sizes -- an
    AMR hierarchy's levels, an octree's leaves -- then count by their size."""
    from tack.dataset.amr import AMRDataSet
    from tack.dataset.octree import OctreeDataSet

    if isinstance(data, OctreeDataSet):
        return data._statistics(field, weights, level)
    if level is not None:
        raise TypeError("level= picks an octree's level; this is not an OctreeDataSet")
    if isinstance(data, AMRDataSet):
        data = data._as_pieces()
    if isinstance(data, PartitionedDataSet):
        fields_of(data, field, "statistics")
        from tack.dataset import units

        # Every rank's pieces' units, so a rank holding none agrees with the rest.
        every = [pair for part in data.world._allgather(
            [_units(piece, field, weights) for piece in data]) for pair in part]
        unit = units.common("statistics' pieces' field", [u for u, _ in every])
        weight_unit = units.common("statistics' pieces' weights", [w for _, w in every])
        moments = None
        for part in _over_pieces(data, lambda piece: _moments(piece, field, weights)):
            if part is not None:
                moments = part if moments is None else _merged(moments, part)
    else:
        unit, weight_unit = _units(data, field, weights)
        moments = _moments(data, field, weights)
    if moments is None:
        raise ValueError("statistics of no values")
    return _summary(*moments)._replace(unit=unit, weight_unit=weight_unit)

Statistics

Bases: NamedTuple

What statistics returns: Python numbers, by name (s.mean), and the field's unit and the weights' (weight_unit: "dimensionless" without weights, a measure's unit for "volume"); unit_of(name) gives each number's.

unit_of 

unit_of(name)

The unit of the statistic name ("mean", "sum", ...), or None if unknown: min, max and mean are in the values' unit, sum in it times the weights', the moment sums m2 ... m4 in its powers times the weights', the variances in its square and the standard deviations in it (a difference of degrees Celsius is in kelvin), skewness, kurtosis and n pure numbers, weight in the weights' unit.

Source code in packages/tack-vis/src/tack/dataset/algorithms/statistics.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def unit_of(self, name):
    """The unit of the statistic ``name`` (``"mean"``, ``"sum"``, ...), or
    ``None`` if unknown: ``min``, ``max`` and ``mean`` are in the values'
    unit, ``sum`` in it times the weights', the moment sums ``m2`` ... ``m4``
    in its powers times the weights', the variances in its square and the
    standard deviations in it (a difference of degrees Celsius is in
    kelvin), ``skewness``, ``kurtosis`` and ``n`` pure numbers, ``weight``
    in the weights' unit."""
    from tack.dataset import units

    u, w = self.unit, self.weight_unit
    if name not in self._fields or name in ("unit", "weight_unit"):
        raise KeyError(f"no statistic {name!r}")
    if name in ("min", "max", "mean"):
        return u
    if name == "sum":
        return u if w == "dimensionless" else units.derived(u, w)
    if name in ("m2", "m3", "m4"):
        return units.derived((units.delta(u), int(name[1])), w)
    if name.endswith("variance"):
        return units.derived((units.delta(u), 2))
    if name.endswith("stddev"):
        return units.delta(u)
    if name == "weight":
        return w
    return "dimensionless"                   # n, skewness, kurtosis

entropy 

entropy(data, field, bins=10, weights=None, level=None)

The Shannon entropy, in bits, of a scalar field's values in bins equal bins over their range: -sum(p log2 p) over the bins' shares, as Viskores' Entropy filter. Ghosts are left out: cells and points data.ghosts marks a duplicate or hidden.

Over a PartitionedDataSet, field names a field every piece has: the range is the whole's, and each piece's histogram over what it owns, made on its own device, is summed, over every rank of data.world. Over an AMRDataSet, the same over every box of every level, refined cells and points left out; over an OctreeDataSet, its leaves (or one level's nodes). With weights ("volume" or a field, as statistics takes them) each bin's share is the weights' sum in it over theirs.

Source code in packages/tack-vis/src/tack/dataset/algorithms/entropy.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
@on_inputs_device
def entropy(data, field, bins=10, weights=None, level=None):
    """The Shannon entropy, in bits, of a scalar field's values in ``bins`` equal bins
    over their range: ``-sum(p log2 p)`` over the bins' shares, as Viskores'
    Entropy filter. Ghosts are left out: cells and points ``data.ghosts``
    marks a duplicate or hidden.

    Over a ``PartitionedDataSet``, ``field`` names a field every piece has: the
    range is the whole's, and each piece's histogram over what it owns, made on
    its own device, is summed, over every rank of ``data.world``. Over an
    ``AMRDataSet``, the same over every box of every level, refined cells and
    points left out; over an ``OctreeDataSet``, its leaves (or one ``level``'s
    nodes). With ``weights`` (``"volume"`` or a field, as ``statistics`` takes
    them) each bin's share is the weights' sum in it over theirs."""
    from tack.dataset.amr import AMRDataSet
    from tack.dataset.octree import OctreeDataSet

    if isinstance(data, OctreeDataSet):
        return data._entropy(field, bins, weights, level)
    if level is not None:
        raise TypeError("level= picks an octree's level; this is not an OctreeDataSet")
    if isinstance(data, AMRDataSet):
        data = data._as_pieces()
    if isinstance(data, PartitionedDataSet):
        fields_of(data, field, "entropy")
        extrema = [e for e in _over_pieces(data, lambda piece: _extrema(piece, field))
                   if e is not None]
        if not extrema:
            raise ValueError("entropy of no values")
        value_range = (float(np.min([e[0] for e in extrema])),
                       float(np.max([e[1] for e in extrema])))

        def counts(piece):
            if _extrema(piece, field) is None:
                return np.zeros(bins)
            return _counts(piece, field, bins, value_range, weights)
        p = sum(_over_pieces(data, counts))
    else:
        if _extrema(data, field) is None:
            raise ValueError("entropy of no values")
        # The histogram's own range, as before: it refuses one that is not finite.
        p = _counts(data, field, bins, None, weights)
    return _bits(p)

Probing, sampling and particles

Values at other points, and paths through a velocity field.

Algorithm Gives Takes VTK Viskores
probe fields at given points or another dataset's points dataset, pieces, AMR, octree vtkProbeFilter, vtkAMRResampleFilter, vtkHyperTreeGridProbeFilter Probe
voxelize the share of each voxel the data covers, and fields' means over it dataset, pieces, AMR, octree — —
advect particles moved through a velocity field by fixed steps dataset, pieces — ParticleAdvection
streamlines the particles' paths, as polylines dataset, pieces vtkStreamTracer Streamline

probe 

probe(data, where, fields=None, locator=None)

data's fields at the points of where -- a dataset, or an (n, 3) array of points, which becomes vertex cells -- as Viskores' Probe and VTK's vtkProbeFilter: the result is where with each field as point data, a field with a basis (H1, L2, Constant, values on points) evaluated through it in the cell holding the point, values on cells taken from that cell, and valid (1, or 0 where no cell holds the point and the values are 0). locator reuses a CellLocator of data.

data may be an AMRDataSet: each point goes to the finest level whose cell holds it (refined cells hold nothing), so a probe onto an image resamples the hierarchy, as VTK's vtkAMRResampleFilter does.

data may be an OctreeDataSet: each point takes the values of the leaf holding it (cell data, constant over the leaf), found by one binary search of the keys, as VTK's vtkHyperTreeGridProbeFilter does; a hidden leaf holds nothing. Over ranks, each point goes to the rank whose trees hold it and the values come back, as for a PartitionedDataSet.

data may be a PartitionedDataSet: each point goes to the piece one of whose own cells (not a ghost) holds it, the lowest such piece, which evaluates it there, and the values come back. where is then this rank's points -- a dataset, an array, or a PartitionedDataSet, which gives one back, each piece probed -- and every rank's points are probed over every rank's pieces of data.world. A point on a face two pieces share takes the lower piece's values, as one dataset takes its lower cell's.

Points given as an array are world positions; a dataset's are its own, relative to its origin. Either way they are taken relative to data's origin in double precision before they are located (the pieces of a PartitionedDataSet share one).

Source code in packages/tack-vis/src/tack/dataset/algorithms/probe.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
@on_inputs_device
def probe(data, where, fields=None, locator=None):
    """``data``'s fields at the points of ``where`` -- a dataset, or an ``(n, 3)``
    array of points, which becomes vertex cells -- as Viskores' Probe and VTK's
    vtkProbeFilter: the result is ``where`` with each field as point data, a
    field with a basis (``H1``, ``L2``, ``Constant``, values on points)
    evaluated through it in the cell holding the point, values on cells taken
    from that cell, and ``valid`` (1, or 0 where no cell holds the point and the
    values are 0). ``locator`` reuses a ``CellLocator`` of ``data``.

    ``data`` may be an ``AMRDataSet``: each point goes to the finest level whose
    cell holds it (refined cells hold nothing), so a probe onto an image
    resamples the hierarchy, as VTK's ``vtkAMRResampleFilter`` does.

    ``data`` may be an ``OctreeDataSet``: each point takes the values of the leaf
    holding it (cell data, constant over the leaf), found by one binary search
    of the keys, as VTK's ``vtkHyperTreeGridProbeFilter`` does; a hidden leaf
    holds nothing. Over ranks, each point goes to the rank whose trees hold it
    and the values come back, as for a ``PartitionedDataSet``.

    ``data`` may be a ``PartitionedDataSet``: each point goes to the piece one
    of whose own cells (not a ghost) holds it, the lowest such piece, which
    evaluates it there, and the values come back. ``where`` is then this
    rank's points -- a dataset, an array, or a ``PartitionedDataSet``, which
    gives one back, each piece probed -- and every rank's points are probed
    over every rank's pieces of ``data.world``. A point on a face two pieces
    share takes the lower piece's values, as one dataset takes its lower
    cell's.

    Points given as an array are world positions; a dataset's are its own,
    relative to its origin. Either way they are taken relative to ``data``'s
    origin in double precision before they are located (the pieces of a
    ``PartitionedDataSet`` share one)."""
    from tack.dataset.amr import AMRDataSet
    from tack.dataset.octree import OctreeDataSet
    from tack.dataset.partitioned import PartitionedDataSet

    if isinstance(data, OctreeDataSet):
        if locator is not None:
            raise ValueError("probe of an OctreeDataSet searches its keys: it takes no "
                             "locator")
        return data._probe(where, fields)
    if isinstance(data, AMRDataSet):
        if locator is not None:
            raise ValueError("probe of an AMRDataSet uses each box's own locator")
        return _probe_pieces(data._as_pieces(), where, fields,
                             levels=[level for level, _ in data])
    if isinstance(data, PartitionedDataSet):
        if locator is not None:
            raise ValueError("probe of a PartitionedDataSet uses each piece's own locator")
        return _probe_pieces(data, where, fields)
    from tack.dataset.ghosts import refuse_hidden, visible

    if locator is None:
        data = visible(data)
    else:
        refuse_hidden(data, "probe with a locator")
    from tack.dataset.carry import _point_values, selected

    if not isinstance(where, DataSet):
        where = _vertices(relative_positions(where, data.origin, data.unit), data)
    if not isinstance(where.geometry.space, H1):
        raise TypeError("probe needs the probed points: an H1 geometry")
    units.check("probe's points and the dataset's spatial unit", where.unit, data.unit)
    locator = locator or cached_locator(data)
    if not locator.serves(data):
        raise ValueError("the locator was built for another topology or geometry, or "
                         "Tack has written the positions since: make a new one")
    points = arrays.materialize(_point_values(where.geometry.values, where.num_points))
    if where.origin != data.origin:
        from tack.dataset.algorithms._moving import _shifted_values

        points = _shifted_values(points, where.dtype,
                                 np.asarray(where.origin) - np.asarray(data.origin))
    found, pcs, tets = locator.find(points)
    m = found.shape[0]
    # Fields read by cell id where find searched from the points; built only
    # then, since it needs the cell bins.
    search = locator.__dict__.get("_mesh_search")
    grouped = []

    def by_cell():
        # The points grouped by the cell found, for the fields evaluated from
        # the cells (through each group's view), made once if any is.
        if not grouped:
            keys = tack.empty((m,), found.dtype)
            if m:
                _found_keys(found, data.num_cells, keys)
            # The extra cell: not found, never evaluated, so not filled.
            grouped.append(bucket_groups(keys, data.num_cells + 1,
                                         ids.at_least(data.id_dtype, m), skip=data.num_cells))
        return grouped[0]

    out = {}
    for name in selected(data, fields):
        f = _field(data, name)
        space = f.space
        if isinstance(space, Values) and space.on in ("points", "cells"):
            f = Field((H1 if space.on == "points" else Constant)(data), f.values, f.unit)
            space = f.space
        if not getattr(space, "interpolated", False) or space.on not in ("points", "cells"):
            continue
        width = arrays.width_of(f.values)
        dtype = arrays.dtype_of(f.values)
        values = (tack.empty((m,), tack.vec(width, dtype)) if width
                  else tack.empty((m,), dtype))
        by_id = field_by_id(data, f) if search is not None else None
        if m and by_id is not None:
            values.fill(0)
            # tets is read only in polyhedra; found stands in for it elsewhere.
            _values_from_points(search, by_id, found, pcs,
                                tets if tets is not None else found, values)
        elif m:
            values.fill(0)
            order, links = by_cell()
            for group in data.launch_groups("cells", [f]):
                if group.count and locator.polyhedral:
                    _evaluate_in_tets(data.domain_view("cells", group), f.view(group), links,
                                      order, tets, pcs, values)
                elif group.count:
                    _evaluate(data.domain_view("cells", group), f.view(group), links, order,
                              pcs, values)
        out[name] = Field(H1(where), values, f.unit)
    valid = tack.empty((m,), tack.i32)
    if m:
        _valid(found, valid)
    out["valid"] = Field(Values(where, "points"), valid)
    kept = dict(where.fields)
    kept.update(out)
    from tack.dataset.ghosts import keep_marks

    return keep_marks(where, DataSet(where.topology, where.geometry, fields=kept,
                                     sets=dict(where.sets)))

voxelize 

voxelize(data, grid, fields=None, samples=4)

data on the voxels (cells) of grid, an image (uniform_grid): the image with fraction, the share of each voxel's volume the data covers, and each of fields' mean over that share (point and cell fields with a basis, all by default), as cell data -- what a simulation's fields are on a uniform grid for analysis or learning. A voxel the data does not reach is hidden (HIDDENCELL), its means 0.

Each voxel is sampled at samples points along each axis it spans (samples ** 3 in 3D, ** 2 on a 2D image), evenly, at the centres of its sub-voxels; a sample is inside where probe finds a cell holding it, and takes probe's value there. So fraction is right to about one sample's share, and a mean is the samples' mean, not an integral; more samples cost more probes. data is anything probe takes: a dataset, pieces over MPI ranks, an AMR hierarchy (its finest level at each sample) or an octree. The image keeps its fields; the result is on its device, at its origin and time. data's and grid's spatial units must agree.

Source code in packages/tack-vis/src/tack/dataset/algorithms/voxelize.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
def voxelize(data, grid, fields=None, samples=4):
    """``data`` on the voxels (cells) of ``grid``, an image (``uniform_grid``): the
    image with ``fraction``, the share of each voxel's volume the data covers, and
    each of ``fields``' mean over that share (point and cell fields with a basis,
    all by default), as cell data -- what a simulation's fields are on a uniform
    grid for analysis or learning. A voxel the data does not reach is hidden
    (``HIDDENCELL``), its means 0.

    Each voxel is sampled at ``samples`` points along each axis it spans
    (``samples ** 3`` in 3D, ``** 2`` on a 2D image), evenly, at the centres of
    its sub-voxels; a sample is inside where ``probe`` finds a cell holding it,
    and takes ``probe``'s value there. So ``fraction`` is right to about one
    sample's share, and a mean is the samples' mean, not an integral; more
    samples cost more probes. ``data`` is anything ``probe`` takes: a dataset,
    pieces over MPI ranks, an AMR hierarchy (its finest level at each sample)
    or an octree. The image keeps its fields; the result is on its device, at
    its origin and time. ``data``'s and ``grid``'s spatial units must agree."""
    from tack.dataset.algorithms._extract import _vertices
    from tack.dataset.algorithms.probe import probe
    from tack.dataset.ghosts import Ghosts, keep_marks
    from tack.dataset.topology import StructuredTopology

    coords = grid.geometry.values
    if not isinstance(coords, UniformCoordinates) or \
            not isinstance(grid.topology, StructuredTopology):
        raise TypeError("voxelize goes onto an image's voxels (uniform_grid)")
    k = int(samples)
    if k < 1:
        raise ValueError(f"samples is at least 1 a side, not {samples!r}")
    like = _like(data)
    units.check("voxelize's data and image spatial units", like.unit, grid.unit)
    cells = [max(d - 1, 1) for d in coords.dims]
    spans = [d > 1 for d in coords.dims]
    ks = [k if s else 1 for s in spans]
    halves = [0.5 if s else 0.0 for s in spans]
    per = int(np.prod(ks))
    n = grid.num_cells
    shift = np.asarray(grid.origin, np.float64) - np.asarray(like.origin, np.float64)
    device = grid.device
    with device:
        fraction = tack.empty((n,), like.dtype)
        hidden = tack.empty((n,), tack.u8)
        means = {}
        batch = max(1, _BATCH // per)
        for v0 in range(0, n, batch):
            nb = min(batch, n - v0)
            points = tack.empty((nb * per,), tack.vec3(like.dtype))
            _samples(coords, v0, cells[0], cells[1], ks[0], ks[1], ks[2], *halves, *shift,
                     per, points)
            probed = probe(data, _vertices(points, like), fields=fields)
            valid = probed.fields["valid"].values
            _fractions(valid, v0, per, fraction, hidden, nb)
            for name, field in probed.fields.items():
                if name == "valid" or name in grid.fields:
                    continue
                values = field.values
                if name not in means:
                    width = arrays.width_of(values)
                    dtype = arrays.dtype_of(values)
                    dtype = dtype if dtype.numpy_dtype.kind == "f" else like.dtype
                    means[name] = (tack.zeros((n,), tack.vec(width, dtype) if width else dtype),
                                   field.unit)
                _means(values, valid, v0, per, means[name][0], nb)
    out = dict(grid.fields)
    space = Values(grid, "cells")
    out["fraction"] = Field(space, fraction, "dimensionless")
    for name, (values, unit) in means.items():
        out[name] = Field(space, values, unit)
    result = keep_marks(grid, DataSet(grid.topology, grid.geometry, fields=out,
                                      sets=dict(grid.sets)))
    # Empty voxels hidden, beside whatever marks the image had; on its device.
    marks = hidden.numpy() * np.uint8(Ghosts.HIDDENCELL)
    if grid.ghosts is not None and grid.ghosts.cells is not None:
        marks |= grid.ghosts.cells.numpy()
    result.ghosts = Ghosts(cells=tack.array(marks.astype(np.uint8), device=device),
                           points=None if grid.ghosts is None else grid.ghosts.points)
    return result

advect 

advect(
    data,
    field,
    seeds,
    step_size,
    steps,
    integrator="rk4",
    locator=None,
)

Seeds moved through the velocity field (a 3-vector field or its name: point data, or cell data) by up to steps fixed steps of step_size, as Viskores' ParticleAdvection: a dataset of the particles where they stopped, one vertex cell each, with steps (taken), time and status (the bits SUCCESS, TERMINATE, SPATIAL_BOUNDS, TOOK_ANY_STEPS, ZERO_VELOCITY). seeds is an (n, 3) array or a field of 3-vectors; integrator "rk4" or "euler"; locator reuses a CellLocator of data.

seeds may also be what advect returned: the particles go on from where they stopped, with their times and the steps they took (steps counts those too), and a particle that stopped -- its steps done, at zero velocity, or outside -- stays as it is. A seed point field on it comes along.

Over a PartitionedDataSet (field by name), particles cross from piece to piece, in rounds: each goes to the lowest piece one of whose own cells (not a ghost) holds it, steps until it leaves that piece, and is handed on, until none is in flight on any rank.

Each particle is one thread that takes all its steps, so on a GPU a few seeds use a few threads of a device built for thousands, each slower than a CPU core: below a few thousand seeds a CPU copy of the mesh is faster (the Visualization Guide's "Integrate a steady velocity field").

Source code in packages/tack-vis/src/tack/dataset/algorithms/advect.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
@on_inputs_device
def advect(data, field, seeds, step_size, steps, integrator="rk4", locator=None):
    """Seeds moved through the velocity ``field`` (a 3-vector field or its name: point
    data, or cell data) by up to ``steps`` fixed steps of ``step_size``, as
    Viskores' ParticleAdvection: a dataset of the particles where they stopped,
    one vertex cell each, with ``steps`` (taken), ``time`` and ``status`` (the
    bits ``SUCCESS``, ``TERMINATE``, ``SPATIAL_BOUNDS``, ``TOOK_ANY_STEPS``,
    ``ZERO_VELOCITY``). ``seeds`` is an ``(n, 3)`` array or a field of
    3-vectors; ``integrator`` ``"rk4"`` or ``"euler"``; ``locator`` reuses a
    ``CellLocator`` of ``data``.

    ``seeds`` may also be what ``advect`` returned: the particles go on from
    where they stopped, with their times and the steps they took (``steps``
    counts those too), and a particle that stopped -- its steps done, at zero
    velocity, or outside -- stays as it is. A ``seed`` point field on it comes
    along.

    Over a ``PartitionedDataSet`` (``field`` by name), particles cross from
    piece to piece, in rounds: each goes to the lowest piece one of whose own
    cells (not a ghost) holds it, steps until it leaves that piece, and is
    handed on, until none is in flight on any rank.

    Each particle is one thread that takes all its steps, so on a GPU a few
    seeds use a few threads of a device built for thousands, each slower than
    a CPU core: below a few thousand seeds a CPU copy of the mesh is faster
    (the Visualization Guide's "Integrate a steady velocity field")."""
    # A step as a time ("0.5 s") is in the spatial unit over the velocity's.
    if _has_quantity(step_size):
        step_size = units.magnitude(step_size, _time_unit(data, field), "advect's step size")
    if isinstance(data, PartitionedDataSet):
        from tack.dataset.algorithms._particles import advect_pieces
        return advect_pieces(data, field, seeds, step_size, steps, integrator)
    from tack.dataset.ghosts import refuse_hidden, visible

    if locator is None:
        data, field = visible(data, field)
    else:
        refuse_hidden(data, "advect with a locator")
    seeds, resume, seed_ids = _resumed(seeds, data)
    xs, times, taken, status, _, _, _ = _run(data, field, seeds, step_size, steps,
                                             integrator, locator, record=False,
                                             resume=resume)
    out = _vertices(xs, data)                  # the particles stay on the device
    points = Values(out, "points")
    for name, values, unit in (("steps", taken, None), ("time", times, _time_unit(data, field)),
                               ("status", status, None)):
        out.fields[name] = Field(points, values, unit)
    if seed_ids is not None:
        out.fields["seed"] = Field(points, seed_ids.to_device(out.device))
    return out

ParticleStatus

Bases: IntFlag

The bits of advect's status field, Viskores' ParticleStatus: SUCCESS (cleared when the seed was outside the mesh), TERMINATE (it stopped after its steps or at zero velocity), SPATIAL_BOUNDS (it left the mesh), TOOK_ANY_STEPS, ZERO_VELOCITY.

streamlines 

streamlines(
    data,
    field,
    seeds,
    step_size,
    steps,
    integrator="rk4",
    locator=None,
    step_unit="time",
)

The paths of advect's particles, as Viskores' Streamline: each seed's positions after every step, the seed first, as one polyline per seed (a PolylineTopology), with seed (which seed's path) on the cells and time on the points. A seed that takes no step has a point and no line. The arguments are advect's.

step_unit is what step_size measures: "time" (the default, as advect and Viskores step), or "length", an arc length along the path: the particle follows the velocity's direction, as VTK's vtkStreamTracer does by default, so slow flow takes as long steps as fast, and steps * step_size is the longest path. time is the time the particle takes either way.

Over a PartitionedDataSet (field by name), paths cross from piece to piece, each piece holding the stretches it traced.

Source code in packages/tack-vis/src/tack/dataset/algorithms/streamlines.py 
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
@on_inputs_device
def streamlines(data, field, seeds, step_size, steps, integrator="rk4", locator=None,
                step_unit="time"):
    """The paths of ``advect``'s particles, as Viskores' Streamline: each seed's
    positions after every step, the seed first, as one polyline per seed
    (a ``PolylineTopology``), with ``seed`` (which seed's path) on the cells and
    ``time`` on the points. A seed that takes no step has a point and no line.
    The arguments are ``advect``'s.

    ``step_unit`` is what ``step_size`` measures: ``"time"`` (the default, as
    ``advect`` and Viskores step), or ``"length"``, an arc length along the
    path: the particle follows the velocity's direction, as VTK's
    ``vtkStreamTracer`` does by default, so slow flow takes as long steps as
    fast, and ``steps * step_size`` is the longest path. ``time`` is the time
    the particle takes either way.

    Over a ``PartitionedDataSet`` (``field`` by name), paths cross from piece to
    piece, each piece holding the stretches it traced."""
    # A step as a time ("0.5 s") is in the spatial unit over the velocity's; as
    # a length ("1 mm"), in the spatial unit.
    if step_unit == "length":
        step_size = _length(data, step_size, "streamlines's step size")
    elif _has_quantity(step_size):
        step_size = units.magnitude(step_size, _time_unit(data, field), "streamlines's step size")
    from tack.dataset.partitioned import PartitionedDataSet

    if isinstance(data, PartitionedDataSet):
        from tack.dataset.algorithms._particles import streamlines_pieces
        return streamlines_pieces(data, field, seeds, step_size, steps, integrator, step_unit)
    from tack.dataset.ghosts import refuse_hidden, visible

    if locator is None:
        data, field = visible(data, field)
    else:
        refuse_hidden(data, "streamlines with a locator")
    from tack.dataset.polyhedra import PolylineTopology

    _, _, _, _, (trace, trace_times), lengths, capacity = _run(
        data, field, seeds, step_size, steps, integrator, locator, record=True, unit=step_unit)
    # Everything is chosen and numbered on the device; only counts come back.
    m = int(lengths.shape[0])
    idt = ids.for_output(data, m * capacity)
    starts = tack.empty((m,), idt)
    total = exclusive_scan(lengths, starts, m) if m else 0
    points = tack.empty((total,), tack.vec(3, data.dtype))
    point_times = tack.empty((total,), data.dtype)
    if m:
        _paths(trace, trace_times, lengths, starts, capacity, points, point_times)
    # The paths that moved, their points numbered as they are, make the chains.
    flags = tack.empty((m,), tack.i32)
    slots = tack.empty((m,), idt)
    if m:
        _moved(lengths, flags)
    lines = exclusive_scan(flags, slots, m) if m else 0
    kept = tack.empty((lines,), idt)
    seed_of = tack.empty((lines,), tack.i32)
    line_lengths = tack.empty((lines,), idt)
    if lines:
        _moved_paths(flags, slots, lengths, kept, seed_of, line_lengths)
    line_offsets = tack.empty((lines + 1,), idt)
    length = exclusive_scan(line_lengths, line_offsets, lines, close=True)
    line_points = tack.empty((length,), idt)
    if lines:
        _chains(kept, starts, lengths, line_offsets, line_points)
    topology = PolylineTopology(line_offsets, line_points, num_points=total, id_dtype=idt)
    out = DataSet(topology, points, **data._frame())
    out.fields["seed"] = Field(Values(out, "cells"), seed_of)
    out.fields["time"] = Field(Values(out, "points"), point_times, _time_unit(data, field))
    return out

Point clouds

Points as vertex cells: neighbours found through a PointLocator.

Algorithm Gives Takes VTK Viskores
extract_points points inside an implicit function dataset vtkExtractPoints ExtractPoints
threshold_points points whose field lies in a range dataset vtkThresholdPoints ThresholdPoints
mask_points every n-th point dataset vtkMaskPoints MaskPoints
point_normals normals from each point's neighbourhood (PCA) dataset vtkPCANormalEstimation —
outliers points marked by a radius or a statistical test dataset, pieces vtkRadiusOutlierRemoval, vtkStatisticalOutlierRemoval —
voxel_downsample one point per occupied voxel dataset, pieces vtkVoxelGrid —
deposit points counted, summed or averaged onto an image (NGP, CIC) dataset, pieces — ParticleDensityNearestGridPoint, ParticleDensityCloudInCell
sph_density SPH densities of particles dataset — —
sph_interpolate particles' fields at other points by SPH kernel sums dataset vtkSPHInterpolator —
halo_catalog a halo catalogue from friends-of-friends labels dataset, pieces — cosmotools

extract_points 

extract_points(data, function, inside=True, fields=None)

The points inside function's region (at or below zero), or outside it (above), as vertex cells. Needs an H1 geometry: one position per point.

The point twin of extract_geometry: the points as vertex cells, with point fields carried and cell fields left behind (a vertex cell has no cell of the input).

Source code in packages/tack-vis/src/tack/dataset/algorithms/extract_points.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
@on_inputs_device
def extract_points(data, function, inside=True, fields=None):
    """The points inside ``function``'s region (at or below zero), or outside it
    (above), as vertex cells. Needs an ``H1`` geometry: one position per point.

    The point twin of ``extract_geometry``: the points as vertex cells, with point fields
    carried and cell fields left behind (a vertex cell has no cell of the input)."""

    if not isinstance(data.geometry.space, H1):
        raise TypeError("extract_points needs an H1 geometry: one position per point")
    values = _point_values(implicit_values(data, function).values, data.num_points)
    # Outside is what is not inside: "above zero" as a range would need a bound
    # between zero and the smallest positive value, which f32 kernels (and GPUs
    # that flush denormals) do not keep.
    flags = _point_flags(data, values, -np.inf, 0.0)
    if not inside:
        _complement(flags)
    return _keep_points(data, flags, fields)

threshold_points 

threshold_points(data, field, lower, upper, fields=None)

The points whose field (point data: H1, values on points) lies in [lower, upper], as vertex cells.

The point twin of threshold: the points as vertex cells, with point fields carried and cell fields left behind (a vertex cell has no cell of the input).

Source code in packages/tack-vis/src/tack/dataset/algorithms/threshold_points.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
@on_inputs_device
def threshold_points(data, field, lower, upper, fields=None):
    """The points whose ``field`` (point data: ``H1``, values on points) lies in
    ``[lower, upper]``, as vertex cells.

    The point twin of ``threshold``: the points as vertex cells, with point fields
    carried and cell fields left behind (a vertex cell has no cell of the input)."""
    field = _field(data, field)
    lower, upper = (_in_unit_of(field, v, "threshold_points' bounds") for v in (lower, upper))
    space = field.space
    if not (isinstance(space, H1) or (isinstance(space, Values) and space.on == "points")):
        raise TypeError(f"threshold_points needs point data, not {space!r}")
    if arrays.width_of(field.values):
        raise TypeError("threshold_points needs a scalar field")
    values = _point_values(field.values, data.num_points)
    return _keep_points(data, _point_flags(data, values, lower, upper), fields)

mask_points 

mask_points(data, stride, fields=None)

Every stride-th point (0, stride, 2 * stride, ...) as vertex cells, as Viskores' MaskPoints.

The point twin of mask: the points as vertex cells, with point fields carried and cell fields left behind (a vertex cell has no cell of the input).

Source code in packages/tack-vis/src/tack/dataset/algorithms/mask_points.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
@on_inputs_device
def mask_points(data, stride, fields=None):
    """Every ``stride``-th point (0, stride, 2 * stride, ...) as vertex cells, as
    Viskores' ``MaskPoints``.

    The point twin of ``mask``: the points as vertex cells, with point fields
    carried and cell fields left behind (a vertex cell has no cell of the input)."""
    if int(stride) < 1:
        raise ValueError(f"stride must be at least 1, not {stride!r}")
    flags = tack.empty((data.num_points,), tack.i32)
    if data.num_points:
        _every(flags, int(stride))
    return _keep_points(data, flags, fields)

point_normals 

point_normals(
    data,
    k=None,
    radius=None,
    viewpoint=None,
    direction=(0.0, 0.0, 1.0),
    periodic=None,
    locator=None,
)

data with normals, each point's unit normal: the direction in which its neighbourhood -- its k nearest points (at most 32), itself among them, or the points within radius -- spreads least, from the eigenvectors of their covariance, as VTK's vtkPCANormalEstimation and PDAL's filters.normal. A point with fewer than three neighbours, or whose neighbours lie on a line, gets a zero normal: so does one whose neighbourhood's two least spreads (variances) differ by less than 4 sqrt(eps) of the greatest, which the closed-form eigenvalues cannot tell from a line (in f32, a strip about 4% as wide as it is long).

A normal's sign is chosen to point toward viewpoint (a world position) when one is given, else along direction (up, by default: the outside of terrain). periodic (a box) takes neighbours across its faces. locator reuses a PointLocator of data.

Source code in packages/tack-vis/src/tack/dataset/algorithms/point_normals.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
@on_inputs_device
def point_normals(data, k=None, radius=None, viewpoint=None, direction=(0.0, 0.0, 1.0),
                  periodic=None, locator=None):
    """``data`` with ``normals``, each point's unit normal: the direction in which
    its neighbourhood -- its ``k`` nearest points (at most 32), itself among
    them, or the points within ``radius`` -- spreads least, from the
    eigenvectors of their covariance, as VTK's vtkPCANormalEstimation and
    PDAL's ``filters.normal``. A point with fewer than three neighbours, or
    whose neighbours lie on a line, gets a zero normal: so does one whose
    neighbourhood's two least spreads (variances) differ by less than
    4 sqrt(eps) of the greatest, which the closed-form eigenvalues cannot tell
    from a line (in f32, a strip about 4% as wide as it is long).

    A normal's sign is chosen to point toward ``viewpoint`` (a world position)
    when one is given, else along ``direction`` (up, by default: the outside
    of terrain). ``periodic`` (a box) takes neighbours across its faces.
    ``locator`` reuses a ``PointLocator`` of ``data``."""
    radius = _length(data, radius, "point_normals's radius")
    periodic = _length(data, periodic, "point_normals's periodic box")
    if (k is None) == (radius is None):
        raise ValueError("point_normals takes k or radius")
    if not isinstance(data.geometry.space, H1):
        raise TypeError("point_normals needs an H1 geometry: one position per point")
    locator = _locator(data, locator, periodic)
    n = data.num_points
    positions = locator.positions
    if k is not None:
        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}")
        nearest, _ = locator.nearest(int(k), positions)      # itself included
        # Rows of k already: CSR offsets 0, k, 2k, ..., stored nowhere.
        offsets = CountingArray(n + 1, start=0, step=int(k), dtype=data.id_dtype)
        ids = nearest.reshape((n * k,))
    else:
        if not radius >= 0:
            raise ValueError(f"radius must be at least zero, not {radius!r}")
        offsets, ids = locator.neighbors(float(radius), positions)
    if viewpoint is not None:
        toward = np.asarray(viewpoint, np.float64) - np.asarray(data.origin)
        toward_point = 1
    else:
        toward = np.asarray(direction, np.float64)
        toward_point = 0
    normals = tack.empty((n,), tack.vec(3, data.dtype))
    if n:
        # Lines' gaps reach 1.22 sqrt(eps); see _smallest_axis.
        tolerance = 4.0 * float(np.sqrt(np.finfo(data.dtype.numpy_dtype).eps))
        _normals(locator.grid(), offsets, ids, *map(float, toward), toward_point, tolerance,
                 normals)
    out = DataSet(data.topology, data.geometry, fields=dict(data.fields),
                  sets=dict(data.sets), dtype=data.dtype, **data._frame())
    out.fields["normals"] = Field(H1(out), normals, "dimensionless")
    from tack.dataset.ghosts import keep_marks

    return keep_marks(data, out)

outliers 

outliers(
    data,
    radius=None,
    min_neighbors=None,
    k=None,
    deviations=1.0,
    periodic=None,
    locator=None,
)

data with each point marked an outlier or not: outlier, a u8 point field, 1 for an outlier. Two tests, as VTK's two filters:

  • radius (radius and min_neighbors): fewer than min_neighbors other points within radius, as vtkRadiusOutlierRemoval; neighbors gives each point's count.
  • statistical (k, at most 32, and deviations): a mean distance to its k nearest other points more than deviations standard deviations above the mean of those means, as vtkStatisticalOutlierRemoval and PCL; mean_distance gives each point's. The mean and deviation are taken in double precision on the host.

The points are marked, not removed, as PDAL marks them (class 7): threshold_points(out, "outlier", 0, 0) keeps the rest. periodic (a box) measures across its faces. locator reuses a PointLocator of data.

Over a PartitionedDataSet of point clouds, each piece is searched with the copies of the other pieces' points it needs (ghost_points: within radius, or its points' k nearest), and the statistical test's mean and deviation are over every piece's points: each point gets the whole's answer. The pieces come back without the copies.

Source code in packages/tack-vis/src/tack/dataset/algorithms/outliers.py 
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
@on_inputs_device
def outliers(data, radius=None, min_neighbors=None, k=None, deviations=1.0, periodic=None,
             locator=None):
    """``data`` with each point marked an outlier or not: ``outlier``, a u8 point
    field, 1 for an outlier. Two tests, as VTK's two filters:

    - **radius** (``radius`` and ``min_neighbors``): fewer than ``min_neighbors``
      other points within ``radius``, as vtkRadiusOutlierRemoval; ``neighbors``
      gives each point's count.
    - **statistical** (``k``, at most 32, and ``deviations``): a mean distance
      to its ``k`` nearest other points more than ``deviations`` standard
      deviations above the mean of those means, as vtkStatisticalOutlierRemoval
      and PCL; ``mean_distance`` gives each point's. The mean and deviation are
      taken in double precision on the host.

    The points are marked, not removed, as PDAL marks them (class 7):
    ``threshold_points(out, "outlier", 0, 0)`` keeps the rest. ``periodic`` (a
    box) measures across its faces. ``locator`` reuses a ``PointLocator`` of
    ``data``.

    Over a ``PartitionedDataSet`` of point clouds, each piece is searched with
    the copies of the other pieces' points it needs (``ghost_points``: within
    ``radius``, or its points' ``k`` nearest), and the statistical test's mean
    and deviation are over every piece's points: each point gets the whole's
    answer. The pieces come back without the copies."""
    radius = _length(data, radius, "outliers's radius")
    periodic = _length(data, periodic, "outliers's periodic box")
    radial = radius is not None or min_neighbors is not None
    if radial == (k is not None):
        raise ValueError("outliers takes radius and min_neighbors, or k (and deviations)")
    if radial and (radius is None or min_neighbors is None or not radius >= 0):
        raise ValueError("the radius test takes a radius of at least zero and "
                         "min_neighbors")
    if not radial and 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}")
    if isinstance(data, PartitionedDataSet):
        if locator is not None:
            raise ValueError("a PointLocator finds one piece's points; pieces take none")
        return _outliers_pieces(data, radius, min_neighbors, k, deviations, periodic)
    if not isinstance(data.geometry.space, H1):
        raise TypeError("outliers needs an H1 geometry: one position per point")
    locator = _locator(data, locator, periodic)
    n = data.num_points
    marks = tack.zeros((n,), tack.u8)
    out = DataSet(data.topology, data.geometry, fields=dict(data.fields),
                  sets=dict(data.sets), dtype=data.dtype, **data._frame())
    points = Values(out, "points")
    if radial:
        locator.check_radius(radius)
        counts = tack.empty((n,), data.id_dtype)
        if n:
            _count_within(locator.grid(), locator.positions, float(radius), 1, counts)
            _below(counts, int(min_neighbors), marks)
        out.fields["neighbors"] = Field(points, counts)
    else:
        means = _means(locator, n, int(k), data.dtype)
        if n:
            host = arrays.to_host(means).astype(np.float64)
            limit = host.mean() + float(deviations) * host.std()
            _above(means, float(limit), marks)
        out.fields["mean_distance"] = Field(points, means, data.unit)
    out.fields["outlier"] = Field(points, marks)
    from tack.dataset.ghosts import keep_marks

    return keep_marks(data, out)

voxel_downsample 

voxel_downsample(data, size, keep='centroid', fields=None)

One point per occupied voxel of a grid of size (a voxel's side, or three sides, in the geometry's units) anchored at the points' lower corner, as VTK's vtkVoxelGrid and PDAL's filters.voxeldownsize: vertex cells, with the input's origin.

keep chooses the point: "centroid", the mean of the voxel's points (new points: floating-point point fields averaged, integer ones -- a classification, intensities -- the voxel's first point's); "first", the voxel's lowest-numbered point; "nearest", its point nearest the centroid (the lowest-numbered on a tie). The last two keep input points, in input order, with their fields and ghost marks. Point fields come along (fields choosing, as filters do); cell fields do not. Each voxel is reduced in index order, so the result has the same bits from run to run, and the voxels and kept points are the same on every backend: the nearest point's distances are rounded alike everywhere (squared_length), so a tie falls the same way.

Over a PartitionedDataSet, the grid is anchored at the lower corner of every piece's points, and each voxel's owned points are brought to one block by a sort over the blocks, in order of piece and then index, and reduced there. A voxel belongs to the piece of its first point, which keeps its centroid; "first" and "nearest" keep input points on their own pieces. Each piece of the result holds its share of the whole's points.

Source code in packages/tack-vis/src/tack/dataset/algorithms/voxel_downsample.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
@on_inputs_device
def voxel_downsample(data, size, keep="centroid", fields=None):
    """One point per occupied voxel of a grid of ``size`` (a voxel's side, or three
    sides, in the geometry's units) anchored at the points' lower corner, as
    VTK's vtkVoxelGrid and PDAL's ``filters.voxeldownsize``: vertex cells, with
    the input's origin.

    ``keep`` chooses the point: ``"centroid"``, the mean of the voxel's points
    (new points: floating-point point fields averaged, integer ones -- a
    classification, intensities -- the voxel's first point's); ``"first"``,
    the voxel's lowest-numbered point; ``"nearest"``, its point nearest the
    centroid (the lowest-numbered on a tie). The last two keep input points,
    in input order, with their fields and ghost marks. Point fields come along
    (``fields`` choosing, as filters do); cell fields do not. Each voxel is
    reduced in index order, so the result has the same bits from run to run,
    and the voxels and kept points are the same on every backend: the
    nearest point's distances are rounded alike everywhere
    (``squared_length``), so a tie falls the same way.

    Over a ``PartitionedDataSet``, the grid is anchored at the lower corner of
    every piece's points, and each voxel's owned points are brought to one
    block by a sort over the blocks, in order of piece and then index, and
    reduced there. A voxel belongs to the piece of its first point, which keeps
    its centroid; ``"first"`` and ``"nearest"`` keep input points on their own
    pieces. Each piece of the result holds its share of the whole's points."""
    size = _length(data, size, "voxel_downsample's size")
    if keep not in _KEEP:
        raise ValueError(f"keep is one of {', '.join(map(repr, _KEEP))}, not {keep!r}")
    sizes = np.broadcast_to(np.asarray(size, np.float64), (3,)).astype(float)
    if not (np.all(np.isfinite(sizes)) and np.all(sizes > 0)):
        raise ValueError(f"size is a positive voxel side, or three, not {size!r}")
    if isinstance(data, PartitionedDataSet):
        return _downsample_pieces(data, size, sizes, keep, fields)
    if not isinstance(data.geometry.space, H1):
        raise TypeError("voxel_downsample needs an H1 geometry: one position per point")
    n = data.num_points
    positions = arrays.materialize(_point_values(data.geometry.values, n))
    if n:
        x = arrays.to_host(positions).reshape(-1, 3).astype(np.float64)
        low, high = x.min(axis=0), x.max(axis=0)
    else:
        low = high = np.zeros(3)
    counts = np.floor((high - low) / sizes).astype(np.int64) + 1
    if float(np.prod(counts.astype(float))) >= 2.0 ** 62:
        raise ValueError(f"{size!r} makes too many voxels over the points' box")
    keys = tack.empty((n,), tack.i64)
    if n:
        _voxel_keys(positions, *map(float, low), *map(float, sizes), int(counts[0]),
                    int(counts[1]), keys)
    order = argsort(keys)
    offsets, nruns = _run_offsets(gather(keys, order), n)
    if keep == "centroid":
        centroids = tack.empty((nruns,), tack.vec(3, data.dtype))
        if nruns:
            _run_means(positions, order, offsets, centroids)
        return _centroids(data, centroids, order, offsets, nruns, fields)
    flags = tack.zeros((n,), tack.i32)
    if nruns and keep == "first":
        _run_firsts(order, offsets, flags)
    elif nruns:
        centroids = tack.empty((nruns,), tack.vec(3, data.dtype))
        _run_means(positions, order, offsets, centroids)
        _run_nearest(positions, centroids, order, offsets, flags)
    return _keep_points(data, flags, fields)

deposit 

deposit(
    data,
    grid,
    field=None,
    reduce=None,
    method="nearest",
    per_volume=False,
)

data's points deposited on grid, an image (uniform_grid, axis- aligned): a field on the image.

method="nearest" puts each point in the cell holding it (nearest grid point, Viskores' ParticleDensityNearestGridPoint): the field is on the cells. method="cic" spreads it over that cell's corners by trilinear weights (cloud in cell, Viskores' ParticleDensityCloudInCell, yt): the field is on the points. Points outside the image are left out.

Without field, each cell's (or corner's) count of points, or sum of weights. With one -- a point field of data or its name -- reduce takes "sum" (the default), "mean", "count", or, for "nearest", "min" and "max": a LiDAR surface model is the "max" of the heights on a 2D image, a terrain model the "min" of the ground's. A cell nothing reaches holds 0 for counts and sums, NaN for the rest. per_volume divides by a cell's volume (area, length on images of fewer axes): a number or mass density. Each target sums its entries in point order, so the result has the same bits from run to run. Counts, minima and maxima are the same on every backend; sums and means agree to rounding (a GPU may fuse a weight's multiply-add).

data and grid may have different origins: the points are moved to the image's by the difference, taken in double precision.

Over a PartitionedDataSet, each piece deposits its owned points (not the copies ghost_points adds) and the pieces' grids combine over every rank, in order of piece: sums and counts add, means divide the sums by the counts, minima and maxima take the least and the greatest. Every rank gets the whole grid's field, on grid's device; field is a name.

Source code in packages/tack-vis/src/tack/dataset/algorithms/deposit.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
@on_inputs_device
def deposit(data, grid, field=None, reduce=None, method="nearest", per_volume=False):
    """``data``'s points deposited on ``grid``, an image (``uniform_grid``, axis-
    aligned): a field on the image.

    ``method="nearest"`` puts each point in the cell holding it (nearest grid
    point, Viskores' ParticleDensityNearestGridPoint): the field is on the
    cells. ``method="cic"`` spreads it over that cell's corners by trilinear
    weights (cloud in cell, Viskores' ParticleDensityCloudInCell, yt): the
    field is on the points. Points outside the image are left out.

    Without ``field``, each cell's (or corner's) count of points, or sum of
    weights. With one -- a point field of ``data`` or its name -- ``reduce``
    takes ``"sum"`` (the default), ``"mean"``, ``"count"``, or, for
    ``"nearest"``, ``"min"`` and ``"max"``: a LiDAR surface model is the
    ``"max"`` of the heights on a 2D image, a terrain model the ``"min"`` of the
    ground's. A cell nothing reaches holds 0 for counts and sums, NaN for the
    rest. ``per_volume`` divides by a cell's volume (area, length on images of
    fewer axes): a number or mass density. Each target sums its entries in
    point order, so the result has the same bits from run to run. Counts,
    minima and maxima are the same on every backend; sums and means agree to
    rounding (a GPU may fuse a weight's multiply-add).

    ``data`` and ``grid`` may have different origins: the points are moved to
    the image's by the difference, taken in double precision.

    Over a ``PartitionedDataSet``, each piece deposits its owned points (not
    the copies ``ghost_points`` adds) and the pieces' grids combine over every
    rank, in order of piece: sums and counts add, means divide the sums by the
    counts, minima and maxima take the least and the greatest. Every rank gets
    the whole grid's field, on ``grid``'s device; ``field`` is a name."""
    if method not in _METHOD:
        raise ValueError(f"method is one of {', '.join(map(repr, _METHOD))}, not {method!r}")
    reduce = reduce or ("sum" if field is not None else "count")
    if reduce not in _REDUCE:
        raise ValueError(f"reduce is one of {', '.join(map(repr, _REDUCE))}, not {reduce!r}")
    if field is None and reduce != "count":
        raise ValueError(f"reduce={reduce!r} needs a field to reduce")
    if method == "cic" and reduce in ("min", "max"):
        raise ValueError("cloud in cell weighs each point: it sums, counts or averages")
    coordinates = grid.geometry.values
    if not isinstance(coordinates, UniformCoordinates) or coordinates.ORIENTED == 1:
        raise TypeError("deposit goes onto an axis-aligned image (uniform_grid)")
    if isinstance(data, PartitionedDataSet):
        return _deposit_pieces(data, grid, field, reduce, method, per_volume)
    unit = _unit(data, grid, field, reduce, per_volume)
    out, on = _deposited(data, grid, field, reduce, method, per_volume)
    return Field(Values(grid, on), out, unit)

sph_density 

sph_density(
    particles,
    kernel="cubic",
    smoothing_length="h",
    mass="mass",
    dimension=3,
    periodic=None,
    locator=None,
)

particles with density, each particle's SPH density rho_i = sum_j m_j W(|x_i - x_j|, h_i) -- itself included, from its own smoothing length (the gather form, as SPLASH computes it) -- for particles that come without one. The arguments are sph_interpolate's.

Source code in packages/tack-vis/src/tack/dataset/algorithms/sph.py 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
@on_inputs_device
def sph_density(particles, kernel="cubic", smoothing_length="h", mass="mass", dimension=3,
                periodic=None, locator=None):
    """``particles`` with ``density``, each particle's SPH density ``rho_i = sum_j m_j
    W(|x_i - x_j|, h_i)`` -- itself included, from its own smoothing length (the
    gather form, as SPLASH computes it) -- for particles that come without one.
    The arguments are ``sph_interpolate``'s."""
    smoothing_length = _smoothing(particles, smoothing_length)
    mass, mass_unit = _given(particles, mass, "mass")
    periodic = _length(particles, periodic, "the periodic box")
    from tack.dataset.ghosts import keep_marks

    w, support = _checked(kernel, dimension)
    if not isinstance(particles.geometry.space, H1):
        raise TypeError("sph_density needs the particles' positions: an H1 geometry")
    n = particles.num_points
    density = tack.zeros((n,), particles.dtype)
    if n:
        found = _locator(particles, locator, periodic)
        smoothing = _per_particle(particles, smoothing_length, "smoothing_length")
        found.check_radius(support * float(arrays.to_host(smoothing).max()))
        _densities(found.grid(), w, support, smoothing,
                   _per_particle(particles, mass, "mass"), density)
    out = DataSet(particles.topology, particles.geometry, fields=dict(particles.fields),
                  sets=dict(particles.sets), dtype=particles.dtype, **particles._frame())
    out.fields["density"] = Field(H1(out), density, units.derived(
        mass_unit, (particles.unit, -dimension)))
    return keep_marks(particles, out)

sph_interpolate 

sph_interpolate(
    particles,
    target,
    fields=None,
    kernel="cubic",
    smoothing_length="h",
    mass="mass",
    density="density",
    normalize=False,
    dimension=3,
    periodic=None,
    locator=None,
)

particles' point fields at target's points by SPH kernel sums: A(x) = sum_j (m_j / rho_j) A_j W(|x - x_j|, h_j), each particle reaching to its kernel's support (2, 2.5 or 3 smoothing lengths), as VTK's vtkSPHInterpolator and SPLASH.

target is a dataset -- an image to resample onto, a slice, a point cloud -- or an (n, 3) array of world positions, which becomes vertex cells; the result is target with each field (fields chooses; all point fields but the smoothing length, mass and density by default) on its points, and shepard, the sum of the particles' (m / rho) W there: where it is near 1 the particles cover the point, near 0 they do not. normalize divides by it (VTK's ShepardProcessing), which makes a constant exact near the particles' edge.

smoothing_length, mass and density are the particles' point fields (names or Fields) or numbers for all of them. kernel is "cubic", "quartic", "quintic" or "wendland", normalized in dimension dimensions. periodic (a box) reaches across its faces. locator reuses a PointLocator of the particles. The target's points are moved to the particles' origin by the difference of the two, in double precision.

Source code in packages/tack-vis/src/tack/dataset/algorithms/sph.py 
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
@on_inputs_device
def sph_interpolate(particles, target, fields=None, kernel="cubic", smoothing_length="h",
                    mass="mass", density="density", normalize=False, dimension=3,
                    periodic=None, locator=None):
    """``particles``' point fields at ``target``'s points by SPH kernel sums: ``A(x) =
    sum_j (m_j / rho_j) A_j W(|x - x_j|, h_j)``, each particle reaching to its
    kernel's support (2, 2.5 or 3 smoothing lengths), as VTK's
    vtkSPHInterpolator and SPLASH.

    ``target`` is a dataset -- an image to resample onto, a slice, a point cloud
    -- or an ``(n, 3)`` array of world positions, which becomes vertex cells;
    the result is ``target`` with each field (``fields`` chooses; all point
    fields but the smoothing length, mass and density by default) on its points,
    and ``shepard``, the sum of the particles' ``(m / rho) W`` there: where it
    is near 1 the particles cover the point, near 0 they do not. ``normalize``
    divides by it (VTK's ShepardProcessing), which makes a constant exact near
    the particles' edge.

    ``smoothing_length``, ``mass`` and ``density`` are the particles' point
    fields (names or ``Field``s) or numbers for all of them. ``kernel`` is
    ``"cubic"``, ``"quartic"``, ``"quintic"`` or ``"wendland"``, normalized in
    ``dimension`` dimensions. ``periodic`` (a box) reaches across its faces.
    ``locator`` reuses a ``PointLocator`` of the particles. The target's points are moved to the particles' origin by the
    difference of the two, in double precision."""
    smoothing_length = _smoothing(particles, smoothing_length)
    mass, mass_unit = _given(particles, mass, "mass")
    density, density_unit = _given(particles, density, "density")
    periodic = _length(particles, periodic, "the periodic box")
    from tack.dataset.algorithms._extract import _vertices
    from tack.dataset.dataset import relative_positions
    from tack.dataset.ghosts import keep_marks

    w, support = _checked(kernel, dimension)
    if not isinstance(particles.geometry.space, H1):
        raise TypeError("sph_interpolate needs the particles' positions: an H1 geometry")
    if not isinstance(target, DataSet):
        target = _vertices(relative_positions(target, particles.origin, particles.unit), particles)
    if not isinstance(target.geometry.space, H1):
        raise TypeError("sph_interpolate needs the target's points: an H1 geometry")
    smoothing = _per_particle(particles, smoothing_length, "smoothing_length")
    reserved = {v for v in (smoothing_length, mass, density) if isinstance(v, str)}
    names = [name for name in (fields if fields is not None else particles.fields)
             if name not in reserved or fields is not None]
    n = particles.num_points
    volumes = tack.empty((n,), particles.dtype)
    if n:
        _volumes(_per_particle(particles, mass, "mass"),
                 _per_particle(particles, density, "density"), volumes)
    units.check("sph_interpolate's target and particles' spatial unit", target.unit,
                particles.unit)
    # Each term is f times m/rho, a volume, times W, per volume: m / (rho L^d) is
    # a pure number for consistent units, unknown where the mass or density is a
    # number.
    weight_unit = units.derived(mass_unit, (density_unit, -1), (particles.unit, -dimension))
    locator = _locator(particles, locator, periodic)
    m = target.num_points
    queries = arrays.materialize(_point_values(target.geometry.values, m))
    if target.origin != particles.origin:
        from tack.dataset.algorithms._moving import _shifted_values

        queries = _shifted_values(queries, target.dtype,
                                  np.asarray(target.origin) - np.asarray(particles.origin))
    if queries.dtype != particles.geometry.values.dtype:
        queries = queries.astype(tack.vec(3, particles.dtype))
    reach = support * (float(arrays.to_host(smoothing).max()) if n else 0.0)
    locator.check_radius(reach)
    out = DataSet(target.topology, target.geometry, fields=dict(target.fields),
                  sets=dict(target.sets), dtype=target.dtype, **target._frame())
    shepard = tack.zeros((m,), particles.dtype)
    grid = locator.grid()
    for name in names:
        field = particles.fields[name]
        on_points = (isinstance(field.space, H1) and field.space.order == 1) or (
            isinstance(field.space, Values) and field.space.on == "points")
        if not on_points:
            continue
        values = arrays.materialize(_point_values(field.values, n))
        scalar = getattr(values.dtype, "scalar_dtype", values.dtype)
        if scalar.numpy_dtype.kind != "f":
            continue                            # an integer field has no kernel average
        result = tack.zeros((m,), values.dtype)
        if n and m:
            _interpolated(grid, w, support, reach, queries, smoothing, volumes, values,
                          result, shepard)
            if normalize:
                _normalized(result, shepard)
        same = normalize or weight_unit == "dimensionless"
        out.fields[name] = Field(H1(out), result, field.unit if same else
                                 units.derived(field.unit, weight_unit))
    if not names and n and m:                      # no fields: the Shepard sum alone
        _interpolated(grid, w, support, reach, queries, smoothing, volumes, volumes,
                      tack.zeros((m,), particles.dtype), shepard)
    out.fields["shepard"] = Field(H1(out), shepard, weight_unit)
    return keep_marks(target, out)

halo_catalog 

halo_catalog(
    particles,
    region="point_region",
    mass=1.0,
    velocity=None,
    center="mbp",
    softening=0.0,
    G=1.0,
    periodic=None,
)

The halos of particles -- the groups of region, a point field of halo labels (connected_components' point_region: friends-of- friends with radius the linking length; -1 for no halo) -- as a point cloud of one point per halo, ported from Viskores' cosmotools and ParaView's LANL halo finder. Each halo's fields:

  • halo, its label; count; mass, the sum of its particles' (mass a point field or one number for all, as Viskores takes it);
  • com, its centre of mass;
  • mbp, its most bound particle (the index of the particle of least potential -G sum m_j / sqrt(r^2 + softening^2) over the halo), and potential, that least potential;
  • with velocity (a 3-vector point field): velocity, the mass- weighted mean, and dispersion, the one-dimensional velocity dispersion sqrt(sum m |v - mean|^2 / 3M);
  • radius, the farthest particle from the halo's point.

The point is the MBP's position (center="mbp") or the centre of mass ("com"). In a periodic box (a side, three, or corners, as PointLocator takes it) each halo is unwrapped about its first particle before anything is summed, so a halo across a face is whole, and its points are put back in the box. Positions are relative to the particles' origin, which the catalogue shares. Every sum runs in particle order, so the result has the same bits from run to run; counts and masses are the same on every backend, the other sums agree to rounding (a GPU may fuse a multiply-add), and so the most bound particle can differ between backends only where two potentials tie to rounding.

Over a PartitionedDataSet (labels numbered across the pieces, as connected_components numbers them), each halo's owned particles are brought to one block by a sort over the blocks, in order of piece and then index, and catalogued there: a PartitionedDataSet of catalogues, the halos in order of label across them, each mbp an index in the piece mbp_piece (a global piece id).

Source code in packages/tack-vis/src/tack/dataset/algorithms/halo_catalog.py 
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
@on_inputs_device
def halo_catalog(particles, region="point_region", mass=1.0, velocity=None, center="mbp",
                 softening=0.0, G=1.0, periodic=None):
    """The halos of ``particles`` -- the groups of ``region``, a point field of
    halo labels (``connected_components``' ``point_region``: friends-of-
    friends with ``radius`` the linking length; -1 for no halo) -- as a point
    cloud of one point per halo, ported from Viskores' cosmotools and ParaView's
    LANL halo finder. Each halo's fields:

    - ``halo``, its label; ``count``; ``mass``, the sum of its particles'
      (``mass`` a point field or one number for all, as Viskores takes it);
    - ``com``, its centre of mass;
    - ``mbp``, its most bound particle (the index of the particle of least
      potential ``-G sum m_j / sqrt(r^2 + softening^2)`` over the halo), and
      ``potential``, that least potential;
    - with ``velocity`` (a 3-vector point field): ``velocity``, the mass-
      weighted mean, and ``dispersion``, the one-dimensional velocity
      dispersion ``sqrt(sum m |v - mean|^2 / 3M)``;
    - ``radius``, the farthest particle from the halo's point.

    The point is the MBP's position (``center="mbp"``) or the centre of mass
    (``"com"``). In a ``periodic`` box (a side, three, or corners, as
    ``PointLocator`` takes it) each halo is unwrapped about its first particle
    before anything is summed, so a halo across a face is whole, and its points
    are put back in the box. Positions are relative to the particles' origin,
    which the catalogue shares. Every sum runs in particle order, so the
    result has the same bits from run to run; counts and masses are the same
    on every backend, the other sums agree to rounding (a GPU may fuse a
    multiply-add), and so the most bound particle can differ between
    backends only where two potentials tie to rounding.

    Over a ``PartitionedDataSet`` (labels numbered across the pieces, as
    ``connected_components`` numbers them), each halo's owned particles are
    brought to one block by a sort over the blocks, in order of piece and then
    index, and catalogued there: a ``PartitionedDataSet`` of catalogues, the
    halos in order of label across them, each ``mbp`` an index in the piece
    ``mbp_piece`` (a global piece id)."""
    softening = _length(particles, softening, "the softening")
    periodic = _length(particles, periodic, "the periodic box")
    # One mass for every particle may bring its unit ("1e10 Msun", (1e10, "Msun")).
    given_mass_unit = None
    if units.quantity(mass) is not None:
        mass, given_mass_unit = units.number_and_unit(mass)
    if center not in _CENTER:
        raise ValueError(f"center is one of {', '.join(map(repr, _CENTER))}, not {center!r}")
    if isinstance(particles, PartitionedDataSet):
        return _catalog_pieces(particles, region, mass, velocity, center, softening, G,
                               periodic, given_mass_unit)
        raise ValueError(f"center is one of {', '.join(map(repr, _CENTER))}, not {center!r}")
    if not isinstance(particles.geometry.space, H1):
        raise TypeError("halo_catalog needs the particles' positions: an H1 geometry")
    n = particles.num_points
    labels_field = particles.fields[region] if isinstance(region, str) else region
    labels = arrays.materialize(_point_values(labels_field.values, n))
    dtype = particles.dtype
    positions = arrays.materialize(_point_values(particles.geometry.values, n))
    box = periodic_box(periodic)
    if box is not None:
        lower = np.asarray(box[0]) - np.asarray(particles.origin)
        sides = np.asarray(box[1]) - np.asarray(box[0])
    else:
        lower, sides = np.zeros(3), np.zeros(3)
    order = argsort(labels)
    sorted_labels = gather(labels, order)
    offsets, nruns = _run_offsets(sorted_labels, n)
    start = 1 if nruns and int(sorted_labels[0]) < 0 else 0
    halos = max(nruns - start, 0)
    masses = _taken(_per_particle(particles, mass, "mass"), order, n)
    u = tack.empty((n,), tack.vec(3, dtype))
    run_of = tack.full((n,), -1, offsets.dtype)
    if halos:
        _unwrapped(positions, order, offsets, start, int(box is not None), *map(float, sides),
                   u, run_of)
    counts = tack.empty((halos,), particles.id_dtype)
    total = tack.empty((halos,), dtype)
    com = tack.empty((halos,), tack.vec(3, dtype))
    if halos:
        _masses(u, masses, offsets, start, counts, total, com)
    # Potentials: every pair on small halos, the pruned search on large ones.
    sizes = np.diff(arrays.to_host(offsets).astype(np.int64)) if nruns else np.zeros(0)
    large = [r for r in range(start, nruns) if sizes[r] > _EXACT]
    exact = tack.array((sizes <= _EXACT).astype(np.int32)) if nruns else tack.zeros((1,),
                                                                                    tack.i32)
    phi = tack.zeros((n,), dtype)
    g, eps2 = float(G), float(softening) ** 2
    if halos:
        _exact_potentials(u, masses, offsets, run_of, exact, g, eps2, phi)
        bounds = arrays.to_host(offsets).astype(np.int64)
        for r in large:
            _pruned(u, masses, int(bounds[r]), int(sizes[r]), g, eps2, phi, dtype)
    entry = tack.empty((halos,), offsets.dtype)
    mbp = tack.empty((halos,), tack.i64)
    potential = tack.empty((halos,), dtype)
    centers = tack.empty((halos,), tack.vec(3, dtype))
    radius = tack.empty((halos,), dtype)
    if halos:
        _least(phi, order, offsets, start, entry, mbp, potential)
        _centers(u, com, entry, int(center == "mbp"), offsets, start, centers, radius)
        _wrapped(centers, int(box is not None), *map(float, lower), *map(float, sides))
        _wrapped(com, int(box is not None), *map(float, lower), *map(float, sides))
    out = _vertices(centers, particles)
    points = Values(out, "points")
    halo_labels = tack.empty((halos,), labels.dtype)
    if halos:
        _halo_labels(sorted_labels, offsets, start, halo_labels)
    # The mass is the particles' mass field's; the potential's G is a number, so
    # its unit is unknown; labels, counts and indices are not quantities.
    mass_unit = _field(particles, mass).unit if isinstance(mass, (str, Field)) else \
        given_mass_unit
    for name, values, unit in (("halo", halo_labels, None), ("count", counts, None),
                               ("mass", total, mass_unit), ("mbp", mbp, None),
                               ("potential", potential, None),
                               ("radius", radius, particles.unit)):
        out.fields[name] = Field(points, values, unit)
    out.fields["com"] = Field(H1(out), com, particles.unit)
    if velocity is not None:
        field = particles.fields[velocity] if isinstance(velocity, str) else velocity
        v = _taken(arrays.materialize(_point_values(field.values, n)), order, n)
        mean = tack.empty((halos,), v.dtype)
        dispersion = tack.empty((halos,), dtype)
        if halos:
            _velocities(v, masses, offsets, start, total, mean, dispersion)
        out.fields["velocity"] = Field(H1(out), mean, field.unit)
        out.fields["dispersion"] = Field(points, dispersion, field.unit)
    return out

Pieces and ghosts

A PartitionedDataSet's pieces, over MPI ranks or not.

Algorithm Gives Takes VTK Viskores
split a structured grid cut into pieces, with ghost layers structured dataset — —
ghost_cells pieces with layers of their neighbours' cells pieces, AMR vtkGhostCellsGenerator —
ghost_points point-cloud pieces with copies of their neighbours' points pieces — —
update_ghosts ghost values refreshed from their owners, in place pieces, AMR — —
remove_ghosts a piece without its ghost cells dataset (a piece) vtkRemoveGhosts GhostCellRemove
generate_global_ids points and cells numbered across every piece pieces vtkGenerateGlobalIds —
merge_pieces the pieces as one dataset, on one rank pieces vtkMergeBlocks MergeDataSets

split 

split(data, pieces, ghost_layers=0)

data, a structured grid, cut by index ranges into a PartitionedDataSet of pieces structured grids, x fastest.

pieces is a count, factored over the axes (each prime factor, largest first, cuts the axis whose pieces have most cells), or a tuple of counts, one per axis of the topology. Each axis's cells are shared out as evenly as they go, every piece getting at least one; neighbouring pieces share the points between them. A piece keeps the geometry's kind -- uniform coordinates with their origin moved, rectilinear axes cut, explicit points taken -- and carries its block of every point, cell and L2 field and every set, as filters carry them; fields on faces or edges are left behind.

ghost_layers adds that many layers of each neighbour's cells around a piece, as far as the grid goes, marked duplicates (DUPLICATECELL in ghosts), as VTK's vtkGhostCellsGenerator adds them. Filters need one to give each piece's share of what they make on the faces and edges pieces share. Every point has one owner: a point of several pieces' own cells is the first such piece's, x fastest, and a duplicate (DUPLICATEPOINT) in the others, as is every point of the ghost layers. The input's own ghost marks come along. An unstructured mesh is partitioned by VTK (vtkRedistributeDataSetFilter) or by its host application.

Source code in packages/tack-vis/src/tack/dataset/algorithms/split.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
@on_inputs_device
def split(data, pieces, ghost_layers=0):
    """``data``, a structured grid, cut by index ranges into a ``PartitionedDataSet``
    of ``pieces`` structured grids, x fastest.

    ``pieces`` is a count, factored over the axes (each prime factor, largest
    first, cuts the axis whose pieces have most cells), or a tuple of counts,
    one per axis of the topology. Each axis's cells are shared out as evenly
    as they go, every piece getting at least one; neighbouring pieces share
    the points between them. A piece keeps the geometry's kind -- uniform
    coordinates with their origin moved, rectilinear axes cut, explicit
    points taken -- and carries its block of every point, cell and L2 field
    and every set, as filters carry them; fields on faces or edges are left
    behind.

    ``ghost_layers`` adds that many layers of each neighbour's cells around a
    piece, as far as the grid goes, marked duplicates (``DUPLICATECELL`` in
    ``ghosts``), as VTK's ``vtkGhostCellsGenerator`` adds them. Filters need
    one to give each piece's share of what they make on the faces and edges
    pieces share. Every point has one owner: a point of several pieces' own
    cells is the first such piece's, x fastest, and a duplicate
    (``DUPLICATEPOINT``) in the others, as is every point of the ghost
    layers. The input's own ghost marks come along. An unstructured mesh is
    partitioned by VTK (``vtkRedistributeDataSetFilter``) or by its host
    application."""
    topology = data.topology
    if not isinstance(topology, StructuredTopology):
        raise TypeError("split cuts structured grids; partition an unstructured mesh "
                        "with VTK (vtkRedistributeDataSetFilter) or the host application")
    layers = int(ghost_layers)
    if layers < 0 or layers != ghost_layers:
        raise ValueError(f"ghost_layers is a count of layers, not {ghost_layers!r}")
    dims = topology.point_dims
    cells = tuple(d - 1 for d in dims)
    counts = _counts(pieces, cells)
    coordinates = data.geometry.values
    implicit = (isinstance(data.geometry.space, H1)
                and isinstance(coordinates, (UniformCoordinates, CartesianProduct)))
    if implicit:
        # The topology's axes are the coordinates' axes of more than one point.
        axes = [a for a, d in enumerate(coordinates.dims) if d != 1]
    device = data.device
    bounds = [[k * n // c for k in range(c + 1)] for n, c in zip(cells, counts)]
    out = []
    for index in np.ndindex(*reversed(counts)):
        index = tuple(reversed(index))                # x fastest
        own = [(b[i], b[i + 1]) for b, i in zip(bounds, index)]        # its cells
        # With its ghost layers, as far as the grid goes.
        start = tuple(max(lo - layers, 0) for lo, _ in own)
        stop = tuple(min(hi + layers, n) for (_, hi), n in zip(own, cells))
        block_cells = tuple(b - a for a, b in zip(start, stop))
        block_points = tuple(n + 1 for n in block_cells)
        piece = StructuredTopology(block_points, id_dtype=topology.id_dtype, device=device,
                                   offset=[o + s for o, s in zip(topology.offset, start)])
        geometry = None
        if implicit:
            geometry = Field(H1(piece), _block_geometry(coordinates, axes, start,
                                                        block_points, device))
        point_ids = _ids(dims, start, block_points, topology.id_dtype)
        cell_ids = _ids(cells, start, block_cells, topology.id_dtype)
        part = carry(data, piece, points=Take(point_ids), cells=Take(cell_ids),
                     geometry=geometry)
        part.ghosts = _block_ghosts(data.ghosts, point_ids, cell_ids, block_points,
                                    block_cells, start, own)
        out.append(part)
    return PartitionedDataSet(out)

ghost_cells 

ghost_cells(
    data, layers=1, world=None, fields=None, tolerance=0.0
)

data's pieces, each with a layer of its neighbours' cells added as ghosts: a new PartitionedDataSet, as VTK's vtkGhostCellsGenerator makes.

data holds this rank's pieces, all unstructured (UnstructuredTopology) or all structured (below), without ghost cells, on any devices; world is the tack.communication.World they live in (data.world, this process alone, by default). Pieces neighbour where they share points at exactly the same position, or with tolerance at most that far apart (in the spatial unit, or a quantity; chains of them too), for pieces whose copies of a shared point differ in their last bits. Each piece gains the cells of other pieces that use a point it shares with them (DUPLICATECELL), and their points it lacks (DUPLICATEPOINT); a point several pieces share belongs to the lowest of them, and is DUPLICATEPOINT in the others. Input marks come along: HIDDENCELL and the other bits are copied with each cell.

Point and cell fields (H1 of order 1, Constant, values on points or cells) are carried, fields choosing them as filters do (None: all); point and cell sets keep their members, and global ids come along where every piece has them, points then matched by their ids rather than their positions. update_ghosts refreshes a field's ghost values from their owners afterwards, after a filter that is right only where it sees every neighbour.

layers layers: the first, every filter's need, is the cells that use a point a piece shares; each further layer the cells, of any piece, that use a point of the layer before (a corner's diagonal neighbours included), as split's ghost layers and vtkGhostCellsGenerator's grow. Pieces thinner than layers cells are fine: a layer reaches past them to the pieces beyond.

Structured pieces (StructuredTopology, each placed by its offset, together tiling a box of cells, as split's and VTK's image pieces do) stay structured: each becomes its block grown by layers cells, as far as the box goes, with its values from the pieces that own them, and its implicit coordinates grown with it.

An AMRDataSet's boxes each grow by layers cells within their level's domain, the ghosts marked DUPLICATECELL and DUPLICATEPOINT: a ghost a sibling box holds is its value, copied (AMReX's FillBoundary); one no sibling holds takes the coarser level's -- the coarse cell's value for a cell, the coarse cell's corners' by linear weights for a point (FillPatchTwoLevels) -- and one neither holds is hidden. A grid's ghost positions are filled the same way, so rectilinear and curvilinear boxes grow too. Ghosts a finer box covers are marked refined. update_ghosts refreshes them, average_down puts fine data into coarse cells.

Source code in packages/tack-vis/src/tack/dataset/algorithms/ghost_cells.py 
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
def ghost_cells(data, layers=1, world=None, fields=None, tolerance=0.0):
    """``data``'s pieces, each with a layer of its neighbours' cells added as ghosts:
    a new ``PartitionedDataSet``, as VTK's ``vtkGhostCellsGenerator`` makes.

    ``data`` holds this rank's pieces, all unstructured
    (``UnstructuredTopology``) or all structured (below), without ghost cells,
    on any devices; ``world`` is the
    ``tack.communication.World`` they live in (``data.world``, this process
    alone, by default). Pieces neighbour where they share points at exactly
    the same position, or with ``tolerance`` at most that far apart (in the
    spatial unit, or a quantity; chains of them too), for pieces whose copies
    of a shared point differ in their last bits. Each piece gains the cells of other pieces that use a
    point it shares with them (``DUPLICATECELL``), and their points it lacks
    (``DUPLICATEPOINT``); a point several pieces share belongs to the lowest
    of them, and is ``DUPLICATEPOINT`` in the others. Input marks come along:
    ``HIDDENCELL`` and the other bits are copied with each cell.

    Point and cell fields (``H1`` of order 1, ``Constant``, values on points
    or cells) are carried, ``fields`` choosing them as filters do (``None``:
    all); point and cell sets keep their members, and global ids come along
    where every piece has them, points then matched by their ids rather than
    their positions. ``update_ghosts`` refreshes
    a field's ghost values from their owners afterwards, after a filter that
    is right only where it sees every neighbour.

    ``layers`` layers: the first, every filter's need, is the cells that use
    a point a piece shares; each further layer the cells, of any piece, that use
    a point of the layer before (a corner's diagonal neighbours included), as
    ``split``'s ghost layers and ``vtkGhostCellsGenerator``'s grow. Pieces
    thinner than ``layers`` cells are fine: a layer reaches past them to the
    pieces beyond.

    Structured pieces (``StructuredTopology``, each placed by its ``offset``,
    together tiling a box of cells, as ``split``'s and VTK's image pieces do)
    stay structured: each becomes its block grown by ``layers`` cells, as far
    as the box goes, with its values from the pieces that own them, and its
    implicit coordinates grown with it.

    An ``AMRDataSet``'s boxes each grow by ``layers`` cells within their level's
    domain, the ghosts marked ``DUPLICATECELL`` and ``DUPLICATEPOINT``: a ghost a
    sibling box holds is its value, copied (AMReX's ``FillBoundary``); one no
    sibling holds takes the coarser level's -- the coarse cell's value for a
    cell, the coarse cell's corners' by linear weights for a point
    (``FillPatchTwoLevels``) -- and one neither holds is hidden. A grid's ghost
    positions are filled the same way, so rectilinear and curvilinear boxes
    grow too. Ghosts a finer box covers are marked refined. ``update_ghosts``
    refreshes them, ``average_down`` puts fine data into coarse cells.
    """
    if isinstance(layers, bool) or not isinstance(layers, (int, np.integer)) or layers < 1:
        raise ValueError(f"ghost_cells makes one or more layers, not {layers!r}")
    from tack.dataset.amr import AMRDataSet

    if isinstance(data, AMRDataSet):
        from tack.dataset.algorithms.amr_ghosts import amr_ghost_cells

        return amr_ghost_cells(data, int(layers), fields)
    if not isinstance(data, PartitionedDataSet):
        raise TypeError(f"ghost_cells takes a PartitionedDataSet, not {type(data).__name__}")
    world = data.world if world is None else world
    common_origin(data, "ghost_cells", world)       # positions are matched as they are
    tolerance = _tolerance(data, tolerance, "ghost_cells")
    pieces = list(data)
    names = _carried(pieces, fields)
    for i, piece in enumerate(pieces):
        _check_piece(i, piece)
    blocks = tc.Blocks([piece.device for piece in pieces], world)
    if blocks.count > 1 << (63 - INDEX_BITS):
        raise ValueError(f"ghost_cells handles up to {1 << (63 - INDEX_BITS)} pieces")
    structured = {g: int(isinstance(p.topology, StructuredTopology))
                  for g, p in zip(blocks.local, pieces)}
    every = blocks.all_gather(structured)
    if blocks.count and all(every):
        out, plans = structured_ghost_cells(blocks, pieces, int(layers), names, _kind)
        result = PartitionedDataSet(out, world=world)
        result._ghost_plans = plans
        return result
    if any(every):
        raise TypeError("ghost_cells takes pieces all structured or all unstructured, "
                        "not a mix")

    hosts = {g: _host(piece) for g, piece in zip(blocks.local, pieces)}
    # Global ids come along where every piece has them, and points are matched
    # by them then, by position otherwise.
    numbered = {kind: bool(blocks.all_reduce({
        g: int(p.global_ids is not None and getattr(p.global_ids, kind) is not None)
        for g, p in zip(blocks.local, pieces)}, "min")) if blocks.count else False
        for kind in ("points", "cells")}

    def entry(g, piece):
        chosen = hosts[g]["candidates"]
        unmarked = np.zeros(len(chosen), np.uint8)
        if numbered["points"]:
            return (chosen, [piece.global_ids.points.numpy()[chosen].astype(np.uint64)],
                    unmarked)
        return chosen, hosts[g]["positions"][chosen], unmarked

    entries = {g: entry(g, piece) for g, piece in zip(blocks.local, pieces)}
    near = tolerance and not numbered["points"]
    partners = (_matching.match_positions(blocks, entries, tolerance) if near
                else _matching.match(blocks, {g: e if numbered["points"] else
                                              (e[0], _matching.position_bits(e[1]), e[2])
                                              for g, e in entries.items()}))
    owners = {g: _matching.owners(g, hosts[g]["num_points"], partners[g],
                                  np.zeros(hosts[g]["num_points"], bool))
              for g in blocks.local}
    if near:
        _matching.settled(blocks, owners, partners)
    needed = _needed(blocks, hosts, partners, owners, int(layers))
    sent = _ghost_parts(blocks, pieces, hosts, needed, owners, names, numbered)
    incoming = blocks.exchange(sent)

    out, receive_points, receive_cells = [], {}, {}
    for g, piece in zip(blocks.local, pieces):
        result, points_from, cells_from = _appended(g, piece, hosts[g], owners[g],
                                                    incoming[g], names, numbered, blocks)
        out.append(result)
        receive_points[g], receive_cells[g] = points_from, cells_from
    point_plan = _matching.plan(blocks, receive_points)
    cell_plan = _matching.plan(blocks, receive_cells)
    result = PartitionedDataSet(out, world=world)
    result._ghost_plans = {"points": point_plan, "cells": cell_plan}
    return result

ghost_points 

ghost_points(
    data,
    radius=None,
    k=None,
    periodic=None,
    around=None,
    fields=None,
    world=None,
)

data's pieces, point clouds (point_cloud: one vertex cell per point), each with copies of the other pieces' points it needs appended as ghosts: a new PartitionedDataSet.

With radius, every point within radius of the piece's box: a filter over the points within radius (outliers' radius test, point_normals, sph_density, sph_interpolate with its kernel's reach) then gives each owned point the whole's answer. around (a PartitionedDataSet with a piece for each of data's) widens each piece's box to its own piece's points too: the targets sph_interpolate resamples onto there. With k (at most 32), every point among an owned point's k nearest, wherever it is: point_normals and outliers by the k nearest then see the whole's neighbours.

Each copy is a vertex cell marked DUPLICATECELL at a point marked DUPLICATEPOINT, owned by the piece it came from, after the piece's own points, which keep their order. Point and cell fields come along (fields choosing, as filters do), sets keep their members, and global ids come along where every piece has them. update_ghosts refreshes a field's copies from their owners, after a per-piece filter that is right only on what a piece owns (sph_density before sph_interpolate).

In a periodic box (as PointLocator takes it) a copy across a face is an image, moved by the box's side: run filters on the result without periodic. data holds this rank's pieces, on any devices, in world (data.world by default), on one origin, without ghosts.

Source code in packages/tack-vis/src/tack/dataset/algorithms/ghost_points.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def ghost_points(data, radius=None, k=None, periodic=None, around=None, fields=None,
                 world=None):
    """``data``'s pieces, point clouds (``point_cloud``: one vertex cell per point),
    each with copies of the other pieces' points it needs appended as ghosts: a
    new ``PartitionedDataSet``.

    With ``radius``, every point within ``radius`` of the piece's box: a filter
    over the points within ``radius`` (``outliers``' radius test,
    ``point_normals``, ``sph_density``, ``sph_interpolate`` with its kernel's
    reach) then gives each owned point the whole's answer. ``around`` (a
    ``PartitionedDataSet`` with a piece for each of ``data``'s) widens each
    piece's box to its own piece's points too: the targets ``sph_interpolate``
    resamples onto there. With ``k`` (at most 32), every point among an owned
    point's ``k`` nearest, wherever it is: ``point_normals`` and ``outliers``
    by the ``k`` nearest then see the whole's neighbours.

    Each copy is a vertex cell marked ``DUPLICATECELL`` at a point marked
    ``DUPLICATEPOINT``, owned by the piece it came from, after the piece's own
    points, which keep their order. Point and cell fields come along
    (``fields`` choosing, as filters do), sets keep their members, and global
    ids come along where every piece has them. ``update_ghosts`` refreshes a
    field's copies from their owners, after a per-piece filter that is right
    only on what a piece owns (``sph_density`` before ``sph_interpolate``).

    In a ``periodic`` box (as ``PointLocator`` takes it) a copy across a face is
    an image, moved by the box's side: run filters on the result without
    ``periodic``. ``data`` holds this rank's pieces, on any devices, in
    ``world`` (``data.world`` by default), on one origin, without ghosts."""
    radius = _length(data, radius, "ghost_points' radius")
    periodic = _length(data, periodic, "ghost_points' periodic box")
    ghosted = _ghosted(data, radius, k, periodic, around, fields, world)
    result = PartitionedDataSet(ghosted.pieces, world=ghosted.blocks.world)
    result._ghost_plans = {"points": ghosted.plan, "cells": ghosted.plan}
    return result

update_ghosts 

update_ghosts(data, *fields)

Copy the named fields' values on ghost points and cells from the pieces that own them, in place, for a PartitionedDataSet that ghost_cells or ghost_points made.

A filter that needs each entity's neighbours -- project to points, point gradients -- is right on what a piece owns and wrong on the ghost layer's rim; refreshing the ghosts makes a second such filter right too. Each named field is a point or cell field of explicit values on every piece. One exchange moves all the point fields, one all the cell fields.

Source code in packages/tack-vis/src/tack/dataset/algorithms/ghost_cells.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
def update_ghosts(data, *fields):
    """Copy the named fields' values on ghost points and cells from the pieces that
    own them, in place, for a ``PartitionedDataSet`` that ``ghost_cells`` or
    ``ghost_points`` made.

    A filter that needs each entity's neighbours -- ``project`` to points, point
    gradients -- is right on what a piece owns and wrong on the ghost layer's
    rim; refreshing the ghosts makes a second such filter right too. Each named
    field is a point or cell field of explicit values on every piece. One
    exchange moves all the point fields, one all the cell fields."""
    from tack.dataset.amr import AMRDataSet

    if isinstance(data, AMRDataSet):
        from tack.dataset.algorithms.amr_ghosts import amr_update_ghosts

        amr_update_ghosts(data, fields)
        return
    plans = getattr(data, "_ghost_plans", None)
    if plans is None:
        raise ValueError("update_ghosts needs the PartitionedDataSet ghost_cells or "
                         "ghost_points made, which keeps the exchanges")
    pieces = list(data)
    blocks = plans["points"].blocks
    by_kind = {"points": {g: {} for g in blocks.local}, "cells": {g: {} for g in blocks.local}}
    for name in fields:
        for g, piece in zip(blocks.local, pieces):
            if name not in piece.fields:
                raise KeyError(f"partition {g - blocks.local.start} has no field {name!r}")
            field = piece.fields[name]
            kind = _kind(field)
            if kind is None:
                raise TypeError(f"update_ghosts refreshes point and cell fields; {name!r} "
                                f"is in {field.space!r}")
            if not isinstance(field.values, tack.Array):
                raise TypeError(f"update_ghosts writes the values in place; {name!r} holds "
                                f"a {type(field.values).__name__}, not a Tack array")
            by_kind[kind][g][name] = field.values
    for kind, values in by_kind.items():
        if any(values.values()):
            plans[kind].run(values)

remove_ghosts 

remove_ghosts(data, fields=None)

data without its ghost cells -- the cells another piece owns (DUPLICATECELL) and the hidden ones (HIDDENCELL) -- and the points no remaining cell uses, as VTK's vtkRemoveGhosts: the last step before rendering, writing one file, or handing a piece to code that does not know ghost marks. Filters keep the ghost layers' output, which the next filter near a piece's edge needs, until this.

The points kept keep their marks: a point a piece shares with another is still that one's, so reductions still count it once. A point this piece owns stays even when only ghost cells used it (a point a filter made on an edge goes to the piece owning the edge's smallest point, whose cells there may all be ghosts), so every point keeps its one owner. A structured grid whose remaining cells are a block stays a structured grid; otherwise the cells are kept as extract_cells keeps them. data itself comes back when it has no ghost cells.

Source code in packages/tack-vis/src/tack/dataset/algorithms/remove_ghosts.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
@on_inputs_device
def remove_ghosts(data, fields=None):
    """``data`` without its ghost cells -- the cells another piece owns
    (``DUPLICATECELL``) and the hidden ones (``HIDDENCELL``) -- and the points no
    remaining cell uses, as VTK's vtkRemoveGhosts: the last step before
    rendering, writing one file, or handing a piece to code that does not know
    ghost marks. Filters keep the ghost layers' output, which the next filter
    near a piece's edge needs, until this.

    The points kept keep their marks: a point a piece shares with another is
    still that one's, so reductions still count it once. A point this piece
    owns stays even when only ghost cells used it (a point a filter made on an
    edge goes to the piece owning the edge's smallest point, whose cells there
    may all be ghosts), so every point keeps its one owner. A structured grid
    whose remaining cells are a block stays a structured grid; otherwise the
    cells are kept as ``extract_cells`` keeps them. ``data`` itself comes back
    when it has no ghost cells."""
    if data._ghosts_lost:
        raise ValueError("this dataset is a filter's output that lost its input's ghost "
                         "marks, so which of its cells are ghosts is not known")
    ghosts = data.ghosts
    if ghosts is None or ghosts.cells is None:
        return data
    keep = tack.empty((data.num_cells,), tack.i32)
    if data.num_cells:
        _unmarked(ghosts.cells, Ghosts.SKIPPED_CELLS, keep)
    owned = None
    if ghosts.points is not None:
        owned = tack.empty((data.num_points,), tack.i32)
        if data.num_points:
            _unmarked(ghosts.points, Ghosts.DUPLICATEPOINT, owned)
    if isinstance(data.topology, StructuredTopology):
        block = _kept_block(data, keep, owned)
        if block is not None:
            return block
    return _keep_cells(data, keep, fields, keep_points=owned)

generate_global_ids 

generate_global_ids(data, world=None, tolerance=0.0)

data's pieces numbered across every piece, as VTK's vtkGenerateGlobalIds: a new PartitionedDataSet of the same pieces with global_ids (tack.dataset.GlobalIds), every point and cell one id, the same on every piece holding it.

Each piece numbers what it owns in its own order, after the pieces before it -- piece 0's owned points from 0, then piece 1's -- and copies take their owner's ids. Owners come from the ghost marks where pieces have them (a copy marked DUPLICATEPOINT or DUPLICATECELL is another piece's); a point several pieces hold unmarked, as pieces without ghost cells share their boundary points, belongs to the lowest of them, and is marked DUPLICATEPOINT in the others, as VTK marks it. Points are copies of each other when their positions are exactly equal -- or, with tolerance, at most that far apart (in the spatial unit, or a quantity), chains of them too, for pieces whose copies of a shared point differ in their last bits -- and cells when their points' global ids are; pieces from ghost_cells reuse the exchanges it kept instead of matching again. Points a piece repeats within itself stay apart: merge_points merges them first.

data holds this rank's pieces, on any devices; world is the tack.communication.World they live in (data.world by default). The matching runs on the host and the exchanges on the pieces' devices.

Source code in packages/tack-vis/src/tack/dataset/algorithms/generate_global_ids.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def generate_global_ids(data, world=None, tolerance=0.0):
    """``data``'s pieces numbered across every piece, as VTK's
    ``vtkGenerateGlobalIds``: a new ``PartitionedDataSet`` of the same pieces with
    ``global_ids`` (``tack.dataset.GlobalIds``), every point and cell one id,
    the same on every piece holding it.

    Each piece numbers what it owns in its own order, after the pieces before
    it -- piece 0's owned points from 0, then piece 1's -- and copies take their
    owner's ids. Owners come from the ghost marks where pieces have them (a
    copy marked ``DUPLICATEPOINT`` or ``DUPLICATECELL`` is another piece's);
    a point several pieces hold unmarked, as pieces without ghost cells share
    their boundary points, belongs to the lowest of them, and is marked
    ``DUPLICATEPOINT`` in the others, as VTK marks it. Points are copies of
    each other when their positions are exactly equal -- or, with
    ``tolerance``, at most that far apart (in the spatial unit, or a
    quantity), chains of them too, for pieces whose copies of a shared point
    differ in their last bits -- and cells when their points' global ids are;
    pieces from ``ghost_cells`` reuse the exchanges it kept instead of
    matching again. Points a piece repeats within itself stay apart:
    ``merge_points`` merges them first.

    ``data`` holds this rank's pieces, on any devices; ``world`` is the
    ``tack.communication.World`` they live in (``data.world`` by default).
    The matching runs on the host and the exchanges on the pieces' devices.
    """
    if not isinstance(data, PartitionedDataSet):
        raise TypeError(f"generate_global_ids takes a PartitionedDataSet, not "
                        f"{type(data).__name__}")
    world = data.world if world is None else world
    common_origin(data, "generate_global_ids", world)
    tolerance = _tolerance(data, tolerance, "generate_global_ids")
    pieces = list(data)
    blocks = tc.Blocks([piece.device for piece in pieces], world)
    plans = getattr(data, "_ghost_plans", None)
    marks = {g: _marks(piece) for g, piece in zip(blocks.local, pieces)}

    # Points: who owns each, then the owners' numbers, then the copies'.
    if plans is not None:
        point_plan = plans["points"]
        point_owned = {g: (marks[g][1] & Ghosts.DUPLICATEPOINT) == 0 for g in blocks.local}
        point_shared = None
    else:
        point_owned, point_plan, point_shared = _point_owners(blocks, pieces, marks,
                                                                 tolerance)
    point_ids = _numbered(blocks, pieces, point_owned, point_plan, "points")

    # Cells: owned cells numbered, ghost cells from their owners.
    cell_owned = {g: (marks[g][0] & Ghosts.DUPLICATECELL) == 0 for g in blocks.local}
    copies = blocks.all_reduce({g: int((~cell_owned[g]).sum()) for g in blocks.local})
    if plans is not None:
        cell_plan = plans["cells"]
    elif copies:
        cell_plan = _cell_owners(blocks, pieces, marks, point_ids, point_shared)
    else:
        cell_plan = None
    cell_ids = _numbered(blocks, pieces, cell_owned, cell_plan, "cells")

    out = []
    for g, piece in zip(blocks.local, pieces):
        points = marks[g][1] | np.where(point_owned[g], 0,
                                        Ghosts.DUPLICATEPOINT).astype(np.uint8)
        with blocks.device(g):
            ghosts = piece.ghosts
            if not np.array_equal(points, marks[g][1]):     # copies newly marked
                ghosts = Ghosts(cells=None if ghosts is None else ghosts.cells, points=points)
            copy = DataSet(piece.topology, piece.geometry, fields=dict(piece.fields),
                           sets=dict(piece.sets), dtype=piece.dtype, ghosts=ghosts,
                           global_ids=GlobalIds(cells=cell_ids[g], points=point_ids[g]),
                           **piece._frame())
        copy._ghosts_lost = piece._ghosts_lost
        out.append(copy)
    result = PartitionedDataSet(out, world=world)
    if plans is not None:
        result._ghost_plans = plans
    return result

merge_pieces 

merge_pieces(
    data, root=0, fields=None, device=None, tolerance=0.0
)

data's pieces as one dataset: every cell a piece owns, once, and every point once, with the point and cell fields (fields choosing them as filters do), as VTK's vtkMergeBlocks after vtkRemoveGhosts.

Copies of a point several pieces hold become one point, matched as generate_global_ids matches them (the ghost marks' owners, else exact positions, or positions at most tolerance apart); ghost cells are left out, and hidden cells kept, marked. The points and cells are numbered as generate_global_ids numbers them: each piece's owned ones after the pieces before it. The pieces' own global ids, if they have them, come along as the result's.

Structured pieces that tile a box of cells (each placed by its topology.offset) give one structured grid of the box -- split's input back -- with its offset and its coordinates (uniform, rectilinear or explicit); other pieces give an unstructured dataset.

Across MPI ranks every rank passes its pieces and the result is on the root rank (None elsewhere), on device (by default the device of the root rank's first piece).

Source code in packages/tack-vis/src/tack/dataset/algorithms/merge_pieces.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def merge_pieces(data, root=0, fields=None, device=None, tolerance=0.0):
    """``data``'s pieces as one dataset: every cell a piece owns, once, and every point
    once, with the point and cell fields (``fields`` choosing them as filters
    do), as VTK's vtkMergeBlocks after vtkRemoveGhosts.

    Copies of a point several pieces hold become one point, matched as
    ``generate_global_ids`` matches them (the ghost marks' owners, else exact
    positions, or positions at most ``tolerance`` apart); ghost cells are left out, and hidden cells kept, marked. The
    points and cells are numbered as ``generate_global_ids`` numbers them:
    each piece's owned ones after the pieces before it. The pieces' own global
    ids, if they have them, come along as the result's.

    Structured pieces that tile a box of cells (each placed by its
    ``topology.offset``) give one structured grid of the box -- ``split``'s
    input back -- with its offset and its coordinates (uniform, rectilinear
    or explicit); other pieces give an unstructured dataset.

    Across MPI ranks every rank passes its pieces and the result is on the
    ``root`` rank (``None`` elsewhere), on ``device`` (by default the device
    of the root rank's first piece)."""
    if not isinstance(data, PartitionedDataSet):
        raise TypeError(f"merge_pieces takes a PartitionedDataSet, not {type(data).__name__}")
    from tack.dataset.algorithms.ghost_cells import _carried, _kind

    pieces = list(data)
    names = _carried(pieces, fields)
    common_origin(data, "merge_pieces")
    tolerance = _tolerance(data, tolerance, "merge_pieces")
    blocks = tc.Blocks([piece.device for piece in pieces], data.world)
    world = blocks.world
    if not 0 <= root < world.size:
        raise ValueError(f"root is a rank, 0 to {world.size - 1}, not {root}")
    firsts = blocks.all_gather({g: g == blocks.local.start for g in blocks.local})
    starts = [g for g, first in enumerate(firsts) if first]
    target = next((g for g in starts if blocks.rank_of(g) == root), None)
    if target is None:
        raise ValueError(f"merge_pieces gathers onto a piece of rank {root}, which has none")
    on_root = target in blocks.local
    out_device = (tack.device(device) if device is not None else
                  blocks.device(target)) if on_root else None
    structured = blocks.all_gather({g: isinstance(p.topology, StructuredTopology)
                                    for g, p in zip(blocks.local, pieces)})
    if all(structured):
        tiled = _box(blocks, pieces)
        if tiled is not None:
            return _merged_structured(blocks, pieces, names, _kind, target, out_device,
                                      tiled)
    return _merged_unstructured(blocks, pieces, names, _kind, target, out_device,
                                tolerance)

Adaptive mesh refinement

Algorithms only an AMRDataSet takes; ghost_cells, probe, statistics and the rest above take one too.

Algorithm Gives Takes VTK Viskores
average_down finer levels' data averaged into the coarse cells they cover, in place AMR — —

average_down 

average_down(amr, field)

Put the finer levels' data into the coarse cells and points they cover, in place, as AMReX's average_down: each refined coarse cell (REFINEDCELL) takes the volume-weighted average of the fine cells over it, each refined coarse point (REFINEDPOINT) the value of the fine point at its location. field names a floating-point point or cell field every box has; finest level first, so a value reaches the coarsest level in one call.

Source code in packages/tack-vis/src/tack/dataset/algorithms/amr_ghosts.py 
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
def average_down(amr, field):
    """Put the finer levels' data into the coarse cells and points they cover, in
    place, as AMReX's ``average_down``: each refined coarse cell (``REFINEDCELL``)
    takes the volume-weighted average of the fine cells over it, each refined
    coarse point (``REFINEDPOINT``) the value of the fine point at its location.
    ``field`` names a floating-point point or cell field every box has; finest
    level first, so a value reaches the coarsest level in one call."""
    from tack.dataset.algorithms.cell_geometry import cell_geometry
    from tack.dataset.algorithms.ghost_cells import _kind

    if not isinstance(amr, AMRDataSet):
        raise TypeError(f"average_down takes an AMRDataSet, not {type(amr).__name__}")
    boxes = [box for _, box in amr]
    kinds = {_kind(box.fields[field]) for box in boxes}
    if len(kinds) != 1 or None in kinds:
        raise TypeError(f"average_down averages a point or cell field; {field!r} is not one")
    kind = kinds.pop()
    layout = _Layout(amr)
    d = amr.dimension
    for level in range(amr.num_levels - 2, -1, -1):
        rows = {}
        for g, (lv, k) in zip(layout.blocks.local, layout.local):
            rows[g] = _down_rows(amr, layout, lv, k, kind, d) if lv == level else \
                (_extent_size(amr._extents[lv][k], kind, d), np.zeros(0, np.int64), [])
        fill = _Fill(layout.blocks, rows)
        values = {g: box.fields[field].values for g, box in zip(layout.blocks.local, boxes)}
        for g, array in values.items():
            scalar = getattr(array.dtype, "scalar_dtype", array.dtype)
            if not isinstance(array, tack.Array) or scalar.numpy_dtype.kind != "f":
                raise TypeError(f"average_down averages floating-point values of explicit "
                                f"arrays; {field!r} is not one")
        if kind == "points":
            fill.run(values)
            continue
        # Cells: sums of measure times value over sums of measure.
        weighted, measures = {}, {}
        for g, box in zip(layout.blocks.local, boxes):
            with box.device:
                measures[g] = arrays.materialize(cell_geometry(box).measure.values)
                weighted[g] = _times(values[g], measures[g])
        fill.run(weighted)
        fill.run(measures)
        for g in layout.blocks.local:
            targets = fill.boxes[g]["targets"]
            if len(targets):
                with layout.blocks.device(g):
                    _divide(weighted[g], measures[g], tack.array(targets, tack.i64), values[g])