Mixed cells plan

Mixed cells: plan

Issue #307, 0.8. Status: in review, #389 (the six phases, clip's docstring, and the changes from Jeff Lee's review), with the test fixes it depends on in #388. The point gradients at bad corners moved to a pull request of their own (vis/point-gradients-bad-corners), since they change every mesh's point gradients. The performance work that followed is on two branches of berkgeveci/tack-fork, stacked on #389 and opened after it merges: perf/vtk-reading (reading VTK grids, on the CPU and a discrete GPU) and perf/analysis (the analysis filters, and the AhmedML sweep benchmark). The AhmedML times on this page were measured with those two branches applied; #389 alone reads and analyses more slowly.

One UnstructuredTopology holds shaped cells and polyhedra together, and, in 2D, shaped cells and polygons. It has one point numbering and one cell numbering. Each face is stored once, with its two sides, whatever kinds of cells it lies between. Filters return one dataset.

Why

OpenFOAM's meshes are mostly hexahedra and wedges, with polyhedra where the mesher refines or snaps. AhmedML's volume_1_appended.vtu has 20,897,282 hexahedra, 66,606 wedges and 508,298 polyhedra. Today Tack reads such a mesh one of two ways:

• Every cell as a polyhedron (vtk_to_dataset's default when a polyhedron is present). This takes 27–36 s and about 19.7 GB to read, and every algorithm then runs its polyhedral path on 21M cells, of which 97% are hexahedra.
• Split by type (polyhedral="split"), into a PartitionedDataSet of a shaped piece and a polyhedral piece. Interface points then exist twice, the pieces' surfaces differ where a hexahedron and a polyhedron share a non-planar face (#337: a sliver no piece holds), and every filter, probe and advection works over pieces.

Mixed cells replace both readings. The shaped cells keep their fast paths, and only the polyhedra pay for being polyhedra.

The data model

td.UnstructuredTopology(
    types, offsets, connectivity,     # as now; types may hold 42 (polyhedron), 7 (polygon)
    num_points=None, id_dtype=None,
    polyhedra=None,                   # the polyhedra's faces, when there are polyhedra:
                                      #   (face_offsets, face_points,
                                      #    cell_face_offsets, cell_faces, cell_face_sides)
)

• types, offsets and connectivity are VTK's. A polyhedron's connectivity is its distinct points and a polygon's is its loop, so the three arrays stay shareable with VTK (copy=False) and with warp-simdata.

• The polyhedra's faces are stored the way PolyhedralTopology stores them, with each face once:

  • face_offsets/face_points give each face's points, in the order that winds the face out of its side 0;
  • cell_faces/cell_face_sides give each polyhedron's faces and which side of each it is on;
  • cell_face_offsets has one entry per cell plus one, and is empty for a shaped cell, as VTK's face locations are.

  A polyhedral view then needs no change: its cell(c) is the global cell id. On AhmedML the per-cell offsets cost 86 MB in i32, against the 19.7 GB of today's all-polyhedra reading.

• Polygons need no extra arrays. A polygon's connectivity is its loop, so its edges come from consecutive points, as a PolygonalTopology's do.

• Every cell is either a shape or a polyhedral cell of its own group. groups() gives one group per shape present, as now, plus a Polyhedron group and a Polygon group of the cells of type 42 and 7. The latter are selected by id (_SelectedPolyhedra; a polygon view of the same kind).

• Checks. validate() and debug mode also check the face arrays (CSR ranges, sides 0/1, at most two sides a face), that only polyhedra have faces, and that a polyhedron's connectivity is its faces' points.

• PolyhedralTopology, PolygonalTopology and PolylineTopology stay as they are. Folding the first two in as the all-polyhedra case is decided once the mixed path works, and goes through the deprecation mechanism below if it is done.

Faces, edges and point links

faces() gives one record for every face of the cells of the highest dimension:

• Triangles and quadrilaterals stay in fixed rows, the record shapes have today (rows, kinds, sides, the per-(cell, local face) arrays). A polyhedron's triangular or quadrilateral face goes into the same match as the shaped cells' faces, by its point set. So a face between a hexahedron and a polyhedron is one face with two sides, and local on the polyhedron's side is its index in cell_faces.
• Faces of five or more points can only be polyhedral. They go in a CSR group of their own (_PolygonFaces), matched as PolygonFaces matches them now.
• A mesh without polyhedra derives its faces exactly as now.
• The per-(cell, local face) arrays (side_face, side_slot, side_orientation) keep their fixed stride for shape groups. For the polyhedral group they are CSR over cell_faces, as PolygonFaces has them; group_starts says which layout a group has.
• T-junctions (#319), where one polyhedral face covers several shaped faces, are designed with this record. A shaped face with no partner keeps one side, and a later step matches it to the larger face. Building that is not part of this plan.
• edges() and point_links() join the shaped groups' corner layout with the polyhedral group's cell-point entries. The polyhedral group's cell points are its connectivity, not derived.

Every caller that tests which record it has (hasattr(faces, "rows"), extract_surface, generate_global_ids, connected_components, sets) reads the groups instead.

Polyhedral or not, per group

About 40 places decide "polyhedral or not" for the whole topology: reference_cells, by_id.NO_REFERENCE and CellLocator.polyhedral. Each becomes a property of a group's shape (shape.REFERENCE, which Polyhedron, Polygon and PolyLine set false). A topology-wide reference_cells stays as "every group has one", for the spaces that need a reference element (H1 above order 1, L2), which refuse a topology with polyhedra as now.

The algorithms

Phase	What	Modules
1	The topology, faces(), edges(), point_links(), groups, checks; vtk_to_dataset and dataset_to_vtk; keeping a subset of cells	topology.py, polyhedra.py, views.py, dataset.py, _extract._keep_cells (so extract_cells, threshold's output and remove_ghosts), interop/vtk
2	Per-cell work per group	gradients, integrate, cell_geometry, cell_centers, values_at_centers, project, cell_quality, surface_normals, threshold's flags, feature_edges, jump, traces, spaces
3	Filters that build a new topology, their shaped and polyhedral parts joined into one	contour and slice (triangles and polygons), clip (shapes and polyhedra), extract_surface, merge_points, triangulate, tetrahedralize, as_polyhedra, as_polygons, check_winding, orient
4	Cells by id: the choice per cell type inside the kernel	by_id, locator, probe, voxelize, advection and streamlines
5	Pieces	ghost_cells, merge_pieces, generate_global_ids, connected_components; polyhedral="split" deprecated (#337)
6	The whole	every call of test_dataset_devices.CALLS on a mixed mesh; AhmedML measured

Done so far. Phase 1: the topology, faces(), edges(), groups, checks, VTK both ways, keeping a subset of cells; AhmedML reads in 10 s (40-54 s as polyhedra). Phase 2: every algorithm in its row, except surface_normals and feature_edges, which take surfaces (phase 3's 2D mixed meshes), and point_links, moved here from phase 1. traces refuses polyhedra with faces of more than four points, which SideTraces cannot hold; point data on a mixed mesh is continuous, so jump and upwind_flux take it without traces. On AhmedML, warm: cell gradients 0.1 s, point gradients 1.1 s, cell_quality 1.1 s, integrate 0.9 s. Phase 3: polygons beside triangles and quads (type 7, edges derived from the loops), and every filter in its row except orient, which the VTK reader already applies. clip keeps whole shapes and makes every other kept cell the polyhedron of its kept part, since a cut shape's table pieces would not match a neighbouring polyhedron's clipped faces; tetrahedralize adds a point per polyhedron and per face of five or more points. Not yet: shrink, subdivide and clip of polygons refuse. (vtk_to_dataset read polygons as a PolygonalTopology until the review of #389; see below.) On AhmedML, warm, mixed against all polyhedra: surface 0.1 s against 0.0, slice 0.1 against 0.1, contour 0.2 against 0.3, clip 0.6 against 0.4, with the same points and measures; peak memory 22 GB against 28. Phase 4: the locator, probe, voxelize, advect and streamlines. Each cell is located by its type, a shape by its coordinates and a polyhedron or polygon in its fan; the per-cell test is compiled in only for mixed meshes, and the optimized code of the search and advection kernels on shaped and polyhedral meshes is unchanged. The walk to face neighbours (advect(..., neighbors=True)) stays off with polyhedra. On AhmedML (240K probes, 200 streamlines), mixed against all polyhedra: locator 0.2 s against 1.3, first probe 0.5 against 0.7, first streamlines 0.6 against 2.8, both warm under 0.05 s; probes in polyhedra equal to rounding, in hexahedra trilinear (up to 2% of the velocity apart); peak memory 17 GB against 25. Phase 5: ghost_cells, update_ghosts, merge_pieces, generate_global_ids and connected_components over pieces, in one process and across MPI ranks. A polyhedron travels with its faces wound out of it and is matched again where it lands, so a piece without polyhedra of its own takes them as ghosts, and a merge gives the mixed mesh back. The deprecation mechanism is in as proposed below, except that deprecated(...) declares a deprecation (in a declaration module per package, so the registry needs no VTK) and the deprecated code calls its warn(); polyhedral="split" is deprecated in 0.8, to be removed in 0.9. Still refused on mixed meshes: shrink, subdivide, and clip of surfaces with polygons. Phase 6: the mixed-mesh refusal is gone. shrink and subdivide refuse a mixed mesh as they refuse polyhedra, for want of a reference element, as H1 above order 1 and L2 do; the sweeps hold every refusal to the all-polyhedra mesh's (in 2D, the all-polygons surface's), and the device sweeps run every call on a mixed mesh and on quads beside polygons, on a second device and on one without f64. On AhmedML a pipeline of gradients, quality, integration, surface, slice, contour, clip, probe and streamlines takes 3.7 s warm mixed against 4.2 s as all polyhedra, after reading in 8.7 s against 32 s, at a peak of 25 GB against 33; the answers agree. The mixed mesh costs what its parts cost alone.

Reading polygons (review of #389, d8e8b85e): vtk_to_dataset reads a grid or vtkPolyData with polygons as an UnstructuredTopology, the polygons beside the triangles and quads, so the reader breaks once, in 0.8, for polyhedra and polygons alike; polyhedral="all" gives the old PolygonalTopology. Whether PolygonalTopology and PolyhedralTopology are still needed is question 4 below, measured together.

Left open: shrink and subdivide on polyhedra, and clip of polygons, are polyhedra and polygon work, not mixed meshes'; advection's walk to face neighbours is off with polyhedra.

Point gradients at bad corners (done after phase 6). AhmedML's tangled and nearly flat hexahedra gave point gradients up to 2e9, on the hexahedra alone too. A cell folded over (a corner's scaled Jacobian below -1e-3) now gives every corner its Green-Gauss gradient over its faces, as a polyhedron's; a corner flat but for rounding (within 1e-3 of 0, f32 rounding of AhmedML's small cells having reached 2e-4) is left out of the point's average, as an exactly flat one was. The largest is now 1.7e3 (a valid corner of scaled Jacobian 0.008), against 132 read as all polyhedra, with a median difference of 1e-4.

Each phase is a pull request. After phase 1 a mixed mesh can be read, written, checked and subset, and every algorithm not yet ported refuses it by name, never silently computing on part of the cells.

Phase 3's joined outputs. Contour and slice produce triangles from the shapes and polygons from the polyhedra, so their output is a 2D mixed topology: triangles beside polygons. Clip produces shapes and polyhedra in one topology. That also opens the way for clip through mesh points (#309), whose 132 table patterns fit only as polyhedra.

Phase 4's cell lookup. Today the mesh class is chosen once for the whole topology (_FlatCells, _PolyhedralMesh, ...). A mixed mesh's kernel reads the cell's type and takes the fan path for a polyhedron. This must not slow the shaped meshes: the check is compiled in only when the topology has polyhedra, and benchmarks/polyhedra_compare.py and the advection benchmarks hold the shaped cost.

Reading and writing VTK

• vtk_to_dataset of an unstructured grid with polyhedra gives a mixed topology:
  • types, offsets and connectivity are shared with VTK where the ids allow;
  • VTK's per-cell face copies are matched once, among the polyhedra only (_match_rings), into the face arrays;
  • orient=True and merge_points= work as they do today, since AhmedML needs both (198,919 faces wound the same way by both cells, and 1.56M duplicate points).
• The old readings stay available: polyhedral="all" reads every cell as a polyhedron, as today's default does. polyhedral="split" is deprecated and removed after a release.
• dataset_to_vtk writes SetPolyhedralCells with each face copied per cell, reversed for side 1, as for a PolyhedralTopology now.
• warp-simdata comes in phase 3 or later, as its users ask.

The default reading of a file with polyhedra changes from a PolyhedralTopology to a mixed UnstructuredTopology. That is a breaking change in the CHANGELOG. polyhedral="all" gives the old result.

Deprecation

Tack has no way to deprecate something yet; it has only removed things outright. This work needs one for polyhedral="split", and later perhaps for PolyhedralTopology. Proposed:

• tack._deprecation (tack-core) with deprecated(what, *, since, removed_in, use). It warns once per call site with TackDeprecationWarning, a subclass of DeprecationWarning (as NumPy's are), at the caller's line.
• A test that fails once a deprecation is due. Every use is recorded in a registry; test_deprecations.py fails once removed_in is at or before the package's version, so a release cannot ship a deprecation past its date.
• Documented. The CHANGELOG gains a ### Deprecated section (tools/changelog.py accepts it), and the docs list current deprecations on one page generated from the registry.
• Policy. Deprecate in one minor release and remove in the next or later. Before 1.0, a breaking change may still go without deprecation where keeping the old behaviour would cost more than it saves; the CHANGELOG says so.

Tests and oracles

• A mixed mesh against the same mesh read two other ways: as all polyhedra (as_polyhedra of its shaped cells, polyhedral="all") and split by type. Per cell and per point, values agree to rounding with the all-polyhedra reading wherever the cell is a polyhedron in both, and with the shaped reading wherever it is a shape in both.
• VTK for what it computes on mixed grids: cell sizes, contours, clips and slices, surfaces.
• Small meshes built for the purpose:
  • a hexahedron next to a polyhedron sharing a quad face;
  • a wedge next to a polyhedron sharing a triangle;
  • two polyhedra sharing a pentagon;
  • a polyhedron whose face covers two hexahedra' faces (a T-junction: refused or one-sided until #319);
  • a 2D mesh of triangles, quads and polygons.
• AhmedML for scale: the memory and time of reading, extract_surface, slices and streamlines, against both of today's readings.
• Every device, through test_dataset_devices' fixtures, as for any dataset.

Questions to settle

1. Must a polyhedron's connectivity be sorted? PolyhedralTopology's derived cell points are sorted and unique. VTK's are only distinct. Taking VTK's as they are keeps them shareable. Does anything rely on the order? (merge_points, ghost matching and keys built from a cell's points need checking.) Answer: no. Per-cell loops take every point (centres, boxes, marks, thresholds); corner keys and point links take each corner's own point; faces are matched by their own points, sorted there (_face_keys in faces(), _match_rings for ghosts, merge_points and the VTK reader), never by a cell's; the locator's weights are made and read in the same order. The mixed-cells tests build their polyhedra in hexahedron-corner order, which is not sorted, and hold them to the sorted all-polyhedra reading. The views' docstring says so (#389's review).
2. Size buckets (SizeBuckets) group polyhedra by their face-point counts. Do they work over a selected group unchanged, or do they need the group's ids? Answer: unchanged. cell_keys keys every cell of the topology by its global id (a mixed topology's cell_face_offsets are per cell, a shape's range empty, so the smallest bucket), and launch_groups splits each group, selected ones included, by its cells' keys. Contour's mixed path (_lopez_polygons(only=)) runs this way.
3. Non-manifold faces (#308): still refused, as now: any face of more than two cells, whatever the cells are. A shaped cell's face on a face two polyhedra share was left a face of its own until #389's review; it is refused too.
4. Folding PolyhedralTopology in, after phase 6: is it worth one class fewer, given the deprecation it would take?