meshio++ roadmap

Status at time of writing: v16.11.0 — 72 core formats plus four Python-only physics-ML ones, thirty-nine mesh operations + six data operations, six language surfaces (Python / C / Fortran / Julia / R / WASM), two viewers plus a browser dataset manager, a Blender add-on, a ParaView plugin, an MCP server, a settings-driven pipeline engine, a dataset-manifest layer with a PhysicsNeMo adapter, and a versioned ABI (MESHIOPLUSPLUS_ABI_VERSION 16).

This document lists what is not built. Nothing here duplicates shipped functionality; where a feature partially exists, the shipped half is named and the gap is stated explicitly. Release history lives in CHANGELOG.md, not here.

Effort key: S = days, M = a couple of weeks, L = a month or more, XL = a project in its own right. An item names a probe — the failing test that proves the gap — wherever one is cheap.

How this file works

• Sections are ordered, and the order is the recommendation. Each section states what belongs in it; an item that sounds exciting does not move up for that reason. An empty section is removed, not kept as a placeholder.
• A closed item is removed, not struck through. Its history is the CHANGELOG.md entry and the feature's own doc/ page; a partly closed item is narrowed to what remains (the AGENTS.md change checklist rule).
• Defect-shaped items go in a Correctness debts section at the top, regardless of size — behaviour that loses data, mis-orients cells, does not terminate on valid input or fails silently is not a feature request, even when the fix and the feature are the same work. None is open today, so there is no such section; the first defect found opens one as §1 and renumbers the rest.
• An item estimated from a doc, a .d.ts or a changelog alone says so ("verify first") and names its probe; the code has repeatedly been more or less capable than its description.
• Non-goals record decisions already taken, with their reasons, so they are not re-proposed as gaps.

The map

---

1. Format reach

Admission: a format a simulation or physics-ML workflow actually exchanges, or the missing half of a shipped one. "Exchanges" means a file that crosses a tool boundary in a real pipeline (mesher → solver, solver → post-processor, solver → training set), not a format that merely exists. Every item names the consumer on the other side of the file.

Sizes: S ≈ days, M ≈ two weeks, L ≈ a month or more, including tests, docs page and CLI wiring. Sizes assume the shared infrastructure in §1.8 exists; the first format that needs a component pays for it.

Link legend: unmarked links were opened or returned by a search while this section was written (September 2026); links marked † are quoted from memory and must be checked before relying on them. Vendor documentation portals move often — when a link is dead, search the document title.

---

A. FEM interchange, solver inputs and result files — native implementations

Python meshio covers none of §1.1–§1.5; these are greenfield. Cross-validation partners are named in each block.

1.1 Tecplot against preplot — S

v16.11.0 reads and writes polygonal and polyhedral zones (Tecplot) and closed Elmer's halo layer and .xplt edge variables. What is left: every .plt fixture, the polyhedral ones included, is written by TecIO rather than by Tecplot's own preplot (no licence was available). Probe: a preplot-converted FE and polyhedral file reading identically to its ASCII source.

1.2 Femap and MFEM, remaining halves — S–M

v16.11.0 writes Femap output sets, reads Femap 2401's property records, checks Patran against real P3/PATRAN and CUBIT exports, and reads MFEM meshes of any order, non-conforming meshes and parallel runs (Patran, Femap, MFEM). What was left out:

• Femap, newer element and group records. The Femap 2401 files found hold no elements, so element (404) and group (408) records from Femap 10 and from 12 on are still unverified; the 13-node pyramid (topology 19) has no documented slot layout; a written file has not been imported into Femap. Probe: a brick20 model exported by two Femap versions reading identically. S
• MFEM NURBS meshes. MFEM NURBS mesh files are refused: each patch would be evaluated per knot span to VTK Lagrange cells of the NURBS order, with rational basis functions and the patches' global control points. Parallel non-conforming meshes, grid functions on non-conforming meshes (their space-filling-curve numbering) and Bernstein or serendipity high-order spaces (read at the vertices) are still out. M
• Patran, other exporters and the rest of the file. No ANSA or HyperMesh export has been read; loads and boundary conditions (packets 06–08, 10) and the .nod/.els/.dis result files stay out until a consumer asks. S

1.3 Abaqus .fil, Radioss, Z88 and libMesh, remaining halves — S

v16.11.0 writes libMesh and reads its compressed files and edge and shell-face sets, reads and writes the whole Z88 deck with every element family's results and Z88Aurora's sets, and applies Radioss units, resolves its boxes, generators and surfaces and reads its animation files (Abaqus .fil, Radioss, Z88, libMesh). What was left out:

• Abaqus .fil against a real binary run. Eleven real ASCII files (pybaqus's tests and example) read, but no public binary .fil with a matching .dat printout exists. The original done-when, a Standard static case whose U and von Mises from S match the .dat, still needs one (ideally with shells and C3D20R). Explicit-only records, contact (15xx), modal and element-matrix records stay out until a consumer asks. S
• Radioss, old decks and the rest. Decks in the 4.x fixed formats (input versions 41 and 44) from real tools have not been seen; boxes in a skew system, per-card unit systems (/UNIT), /SURF/BOX, /PLANE and /ELLIPS, the engine deck, the T01 time history and animation files of older layouts (other magic numbers) are not read. S
• Z88, unrun element types. The layered shells 21 and 22 have not been run through Z88R, and a flat deck of 8-node shells (type 23) fails Z88R's own Jacobian check for reasons not found; the 16-node plate (type 19, read as corner quads) cannot be written back; surface loads (z88i5.txt) are neither read nor written, and Z88Aurora's sets are not written. S
• libMesh, refinement trees and native compression. The writer writes a refined mesh flat, as its active cells; bzip2 files and compressed output go through Python only (the native reader inflates gzip). S

1.4 Marc and Ansys .rst, remaining halves — S each

v16.8.0 reads Marc input decks and formatted post files (Marc) and completes the .rst reader with element results, reaction forces, distributed solves and the full rotor of static cyclic models (Ansys results). What was left out:

• Marc against a real run, and the rest of the deck. No Marc licence was available: the deck reader is checked on real Marc Mentat 2020 decks (DAMASK's element library, AGPL, not committed) and the post-file reader on one real .t19. Edge and face sets (elem:face members), INCLUDE files, the Herrmann, interface, gasket and rebar element types and the 15-node pentahedron (whose element number was not confirmed) are not read; remeshing increments and post files before revision 9 are refused; the binary .t16 stays on the PyPost route (§1.7). S
• .rst compressed records. /FCOMP,RST,1 and higher are refused: detail/zlib_inflate exists, the record framing around the stream needs a sample. S
• Modal cyclic expansion. ansys_rst_cyclic expands static analyses; a modal one needs each mode combined with its pair (or the duplicate sector) by harmonic index, as pymapdl-reader's CyclicResult does. Probe: pymapdl-reader's academic rotor, whose full-rotor mode shapes meshio++ cannot produce today. S
• Other element records (energies ENG, fluxes and gradients EFX/EGR, nonlinear data ENL, contact ECT) as consumers ask.

1.5 LS-DYNA d3plot and Nastran OP2, remaining halves — S–M each

v16.9.0 reads the d3plot family and OP2 files (LS-DYNA d3plot, Nastran OP2). What was left out:

• d3plot against LS-DYNA itself. No LS-DYNA licence was available: the reader is checked on lasso-python's real test families and on families lasso-python's writer produced. No real file with element deletion events, rigid bodies, rigid roads, SPH or airbags has been read — those blocks are sized from the manual and lasso-python and skipped, not read — and it is unknown whether LS-DYNA writes the ten-node solids' extra nodes once (the 2006 manual) or twice (lasso-python). Probe: a public crash model's family that deletes elements, read against LS-PrePost. S
• d3plot, the rest. Adaptive remeshing (NADAPT, a geometry per adaptive state; each state would carry its own mesh through the sequence engine), the 20/27-node and cubic solids, SPH and airbag particles as point clouds, rigid-body motion, and the d3part, d3thdt and binout (LSDA) siblings; femzip needs its closed library and stays out. M
• OP2, other tables. Complex (frequency response, complex modes), random and SORT2 tables; element forces OEF, strain energies ONRGY; springs and NX's newer solid types (300–303); PARAM,POST,-2 files, which are read only by extension (no header to sniff). A file with neither GEOM1 nor a sibling deck could take its points from BGPDT/EQEXIN. S–M

B. Formats that need a vendor runtime or a heavy optional dependency

1.6 DOLFINx ADIOS2 .bp (read; optional ADIOS2 build) — M

VTXWriter output: an ADIOS2 BP4/BP5 directory with a vtk.xml schema attribute and step-wise geometry, connectivity, types, NumberOfNodes, NumberOfEntities plus field variables. ADIOS2 is a heavy dependency for one consumer, so this is an opt-in build flag, never a default. adios4dolfinx checkpoints and Fides output are different layouts — out of scope.

• Done when. A DOLFINx demo's .bp reads as a sequence matching ParaView's VTX reader.
• References. DOLFINx † · ADIOS2 documentation † · VTK ADIOS2 module (VTX schema reader) · VTX reader changes in VTK 9.4 · adios4dolfinx †

1.7 Vendor-runtime routes (exporter scripts and optional plugins, no native parser) — M in total

For these, the deliverable is a documented, tested route — a script that runs inside the vendor's own Python and writes VTKHDF/XDMF, or a plugin compiled against an SDK the user already owns. meshio++ never redistributes vendor libraries.

• Abaqus .odb. No public specification; readable only through the ODB API on a licensed install, version-locked (abaqus upgrade -odb). Ship an abaqus python exporter script (instances → regions, steps/frames → sequence, field outputs by position) and document it; an optional C++ plugin against the ODB API comes only if a user asks. Where no licence is available, the .fil results file (Abaqus .fil, v16.7.0) is the fallback. — ODB2VTK (reference exporter) · abqpy †
• MSC Marc .t16. Read through PyPost (py_post), shipped with Marc/Mentat. Ship an exporter script; the formatted .t19 is read natively (Marc, v16.8.0).
• Femap .modfem. COM API on Windows only. The route is "export a neutral file", which meshio++ reads (Femap).
• Tecplot .szplt. Undocumented; TecIO only. Optional link against a user-installed TecIO, or "re-save as .plt", which is read natively (Tecplot, v16.10.0). — TecIO †
• ANSYS results beyond the native reader (§1.4). Document the PyDPF export route.
• Done when. Each route has a docs page, a script under contrib/, and one manual test recorded with the vendor version used.

---

1.8 Shared infrastructure worth building once

• Fixed-width card tokenizer, remaining adoption. detail/keyword_card.hpp (v15.2.0, built for LS-DYNA) splits a card into standard, long, I10 and free layouts per line, parses real fields with Fortran spellings and, since v16.3.0, parses a Fortran format line such as (3i9,6e21.13e3) into the fields split_fixed cuts (the ANSYS .cdb reader's). The Nastran/OptiStruct (v16.1.0) and Patran (v16.5.0) readers use it too, and the Radioss (v16.7.0) and Marc (v16.8.0) readers its field converters (their column widths follow the deck: Radioss's input version, Marc's EXTENDED); UNV should adopt it instead of growing its own (the .frd reader of v15.3.0 slices its own columns: the layout is fixed by the record key, not by a format line).
• Node-ordering permutation registry, remaining half. detail/node_order.hpp and _node_order.py (v16.0.0, built for Code_Aster .mail) hold the MED, Code_Aster, .frd, UNV, COMSOL, Elmer, FEBio, FLUX, Patran, libMesh, Radioss and Z88 tables in both directions, with a self-test that maps a reference element through every table and checks midpoints and a positive Jacobian (node ordering). What remains: the gmsh, CGNS, GiD, Exodus and Kratos tables still live inside their readers and should move in when those formats are next touched (Marc, v16.8.0, needs no entry: every Marc element type read is in meshio++'s node order; ANSYS's degenerate shapes and Femap's 20-slot brick layout are slot maps, not permutations, and live in their readers; MFEM needs none, since its order-2 nodes are placed by the entity they sit on).
• Fortran unformatted-record reader, remaining adoption. detail/fortran_records.hpp (v16.7.0, built for Abaqus .fil) sniffs 4/8-byte markers in either byte order and splits a file into its records, and detail/binary_stream.hpp reads values in a chosen byte order (libMesh .xdr). The OP2 reader (v16.9.0) builds on it, and the Python twins share _fortran_records.py; EnSight Fortran binary should too. The ANSYS .rst reader addresses records by pointer and keeps its own view, and d3plot (v16.9.0) is word-addressed, not record-framed.
• Directory-as-format support, remaining half. Reading sniffs a directory by its files since v16.2.0 (Elmer, OpenFOAM). What remains: sequence globs keep only suffixed directories (.pmsh, .zarr), so an Elmer directory must be listed explicitly, and ADIOS2 .bp (§1.6) will need its own rule.
• Optional-plugin mechanism and a contrib/ script convention — §1.6, §1.7.

1.9 Suggested order

1. Next FEM wave: §1.1, §1.2, §1.3, §1.4 and §1.5 as Tecplot, Femap, MFEM, Patran, Abaqus, Radioss, Z88, libMesh, Marc and Ansys users ask (the .fil check against a real binary Abaqus run first, and a real Marc run, as soon as one is in hand; MFEM NURBS when an isogeometric user asks).
2. Heavy binaries, on demand: the rest of §1.5 when a d3plot or OP2 user brings a file.
3. Routes: §1.7 scripts as users ask; §1.6 when a FEniCSx user asks.

Considered, not queued

• Graphics, CAD, VFX or lab-specific rather than simulation interchange: 3MF, X3D, DXF, Alembic, OpenVDB, Silo, LAS/LAZ, E57, OBJ materials (.mtl).
• Closed FEM databases with no route worth shipping: Altair .h3d (undisclosed format, SDK not redistributable — Altair's own advice for OpenRadioss is ANIM → VTK, see the OpenRadioss discussion); HyperMesh .hm; Abaqus .sim.
• No native format to implement: deal.II (reads and writes UCD, VTK, Gmsh, Abaqus, Exodus — all reachable today); MOOSE (Exodus).
• Nastran OP4 matrices and punch .pch: matrix and card dumps, not mesh/field exchange; pyNastran serves them.
• Non-MSC Nastran HDF5 dialects: incompatible schemas; revisit with a file in hand.
• Partitioning niche: gmsh $PartitionedEntities (Kratos partitions through MDPA).
• Operators and tooling, not formats (listed only because they came up in the same review): ICP registration; mesh duals; cylindrical/spherical coordinate transforms; shell completions.

Revisit any of them when a consumer asks with a file in hand — a real deck or result file, the tool version that wrote it, and the tool that needs the converted output.

Open verification items before coding

• The Marc element number of the 15-node pentahedron, and the node counts of the Herrmann element types, against Volume B (§1.4).
• Every link marked †.

---

2. Spack package upkeep

Admission: a package that already exists but lives in someone else's repository and has fallen behind the release. Not a feature, and ordered straight after format reach because it is how HPC users get meshio++ at all (issue #3: "Interesting for HPC"; the maintainer's follow-up: done, but it has to be kept current with every release).

Both recipes are upstream in spack/spack-packages — meshioplusplus (a CMakePackage: C API, Fortran, installable C++ API, CLI) and py-meshioplusplus (a PythonPackage), added in PR #5624. Nothing in this repository builds or tests them, so nothing notices when they drift.

• Verify first, then bump the versions. Read on 2026-09-21, the newest tagged version() in both recipes is 9.10.0 (plus master) against 15.2.0 here, and the meshioplusplus recipe's url still points at the v9.10.0 tarball. Probe: spack versions meshioplusplus py-meshioplusplus against git tag. Add a version(..., sha256=...) line per release worth keeping (spack checksum computes it); spack install py-meshioplusplus@15.2.0 and meshioplusplus@15.2.0 must then concretize and build. S
• Audit the recipes against six major versions of drift. Every when="@X:" guard, the depends_on floors (py-scikit-build-core@0.8:, py-pybind11@2.11:, python@3.8:, cmake@3.15:) and the conflicts("%gcc@:9") must be re-checked against pyproject.toml and CMakeLists.txt at the new version, and the docstring's "~40 unstructured mesh formats" is now 52. The +cxx_api libraries install as libmeshioplusplus_core_<backend>.so.<abi> and MESHIOPLUSPLUS_ABI_VERSION is now 15, while the C library keeps SOVERSION 0 (ABI policy) — confirm the recipe finds both, and that a find_package(meshioplusplus X.Y.Z EXACT CONFIG) consumer against a Spack-installed +cxx_api still resolves (C++ API). S
• Variants that lag the CMake options. The recipes expose hdf5, netcdf, zlib, zstd, lz4, kahip, fortran, cli, parallel, mesh_backend, cxx_api and cxx_api_backends. CMakeLists.txt also defines MESHIOPLUSPLUS_WITH_CGNSLIB, MESHIOPLUSPLUS_WITH_GIDPOST, MESHIOPLUSPLUS_WITH_POLYSCOPE, MESHIOPLUSPLUS_WITH_EIGEN and MESHIOPLUSPLUS_WITH_JSON. Verify which of them post-date 9.10.0 and which have a Spack package to depend on (cgns does); expose the ones an HPC build would choose, leave the rest at their defaults, and keep every variant named as the Conan option and vcpkg feature are (C API). S–M
• Test the matrix, not one install. Build both recipes with spack install --test=root for the default variants, +fortran, +cxx_api (each cxx_api_backends), parallel=openmp|tbb|kokkos and mesh_backend=kratos, and run spack style and spack audit before opening the PR. None of this is in CI here today; a smoke job that installs from a spack-packages checkout on a schedule would catch the next drift without a person remembering to look. S–M
• Make the bump a release step. Installation → Spack and the C API page both say a new release "needs no action here" because a version(...) line is added upstream after each tag; the drift above is what that sentence produces when it is nobody's job. Add the upstream PR to the release checklist in AGENTS.md (after the tag: spack checksum, one version() line per recipe, spack style) and replace the sentence in both docs. Opening the PR from the release workflow needs a token on a fork and is a follow-up only if the manual step is skipped again. S
• Done when. Both recipes list the current release, spack install py-meshioplusplus@<current> +hdf5 +netcdf +zlib and spack install meshioplusplus@<current> +fortran +cxx_api succeed on a clean Spack, the docs no longer claim the step is automatic, and issue #3 is closed.
• References. Spack packaging guide † · spack checksum † · PR #5624, the original submission · Installation → Spack

---

3. Quality of implementation

Admission: work that makes every other item safer to land. None of it is a feature, so none of it competes for the same attention — it can run in parallel with everything.

• A sanitizer CI leg — ASan and UBSan over the existing cpp-tests job. No workflow passes -fsanitize today, and a fuzzer that finds a crash without one reports a symptom rather than the out-of-range read behind it. The precondition for the next item. S
• Fuzzing the readers (libFuzzer, then OSS-Fuzz if the project is accepted). 43 mostly hand-rolled parsers are reachable from a C ABI, a browser, a VS Code extension and an MCP server — untrusted input reaches them by design. One fuzz target per registry_readers() entry, seeded from tests/python/meshes/. The highest-value non-feature item in this document. M
• A format conformance matrix — one canonical mesh written to and read back from every writable format, asserting per format what survives (points, each cell type, point/cell/field data and their dtypes, each region kind) against a declared expectation. tests/python/test_region_roundtrip.py already does this for regions over Gmsh/Abaqus/MED, and tests/cpp/test_sequence.cpp's WriteSupportsTimeAgreesWithReality is the registry-iterating shape to generalise it to. The declared expectations become a lossiness column in the format table, which today has only Read/Write/dependencies, with lossiness scattered across its notes and fifty-five per-format quirks sections. The canonical mesh should be a primitive constructor from §6. M
• Property-based testing (Hypothesis) over the invariants the docs already articulate: partition-of-unity, volume conservation, conformity, byte-identical determinism, map composition. M
• A benchmark harness that covers what ships, with a CI leg. The suite exists (benchmark/, up to ~1M synthetic tets in Python, 257k in the C++ backend benchmark, which is off by default) but no CI job runs any of it, so a performance regression is found by a user. It is also narrow: benchmark/bench.py times 6 format labels of the 43 the core reads, and benchmarks has no numbers for any operation or for the parallel backends — src/cpp/benchmark/bench_backends.cpp compares mesh backends only. Widen bench.py to every registry format; add a bench_ops.cpp (extract_surface, smooth, refine, merge, clean, compute_sdf, decimate, partition, reorder) over a size sweep and SEQ/OpenMP/TBB/Kokkos; add a 10M+ cell tier; run it on a schedule that records rather than gates. Every §4 item is gated on this showing its before/after, and it decides whether the scale items in §8 matter at all. S–M
• Test and install the ParaView plugin. tools/paraview-meshioplusplus-plugin.py ships as a reader and writer, but nothing tests it, and the data_files entry that would install it is commented out in pyproject.toml, so its page describes a plugin path nothing writes. A pvpython smoke step plus the install fix. S
• Finish the fallback narrowing. The per-format shims route every decline through core_declined (_fallback.py), but two halves of the same defect remain. First, 44 broad except Exception handlers in 36 package-root operation files (_clean.py, _data_average.py, _curvature.py, …) still wrap the C++ core and, unlike a format fallback, silently substitute a different algorithm; each operation needs its own ruling on which exceptions mean "unsupported" and which mean "bad input", with _error.py's (ValueError, TypeError) re-raise as the model. Second, the C++ core has one ReadError for both "malformed file" and "construct I deliberately decline", and about a hundred std::stoi/stoll/.at() call sites in the format readers leak ValueError/IndexError instead (tests/python/meshes/tecplot/quad_zone_space.tec makes the C++ Tecplot reader throw std::stoull, and now logs a warning before the Python twin reads it); a ReadError subclass, or a per-entry-point wrapper in _core.cpp (not the global translator, which would remap the operations' std::invalid_argument), makes the "recognised decline" contract total. Running the suite once under MESHIOPLUSPLUS_STRICT_CORE=1 sizes both. M

---

4. Performance

Admission: a measured or code-verified slowdown in a path a user hits, with the shape of the fix named. Nothing here is scheduled before the §3 harness can show its before/after.

Two findings frame the section. First, the serial phases below are deliberate: each is documented in the code as a determinism pin, not an oversight — output is byte-identical across parallel backends and thread counts, and the reference-file tests enforce it — so every fix must keep that guarantee and prove it with a SEQ-versus-OpenMP diff, not assert it. Second, every parallel item is conditional on the backend: a SEQ build (and the stl fallback without TBB) runs parallel_for sequentially, so each change must also show that SEQ does not get slower — a parallel sort is O(n log n) where the hash map it replaces is O(n).

Measured regressions (verify first: benchmark/results.csv is one run on one machine). Gmsh binary write at 0.68× legacy meshio and MED read at 0.81–0.93×, both reported as parity or better in benchmarks until this change; XDMF (HDF5) write at 0.92–0.97×; XDMF read at 1.0× on a single-block mesh against 10× on mixed topology, which suggests the post-read conversions — xdmf.cpp has no parallel_for at all. Re-measure on the widened harness, then profile. S each to size

Text I/O.

• A shared tokenizer and number path. 29 readers split each line with std::istringstream and >> into a std::vector<std::string> — one stream and N heap allocations per line (su2.cpp, vtk_read.cpp, mdpa.cpp, avsucd.cpp, tecplot.cpp, unv.cpp, flac3d.cpp and 22 more) — while gmsh.cpp's GmshCursor, a string_view cursor calling strtod straight on the buffer, is the in-repo model to copy. Built on the shipped detail/fast_number.hpp (locale-independent parse_double/snprintf_c) and migrated reader by reader against the reference files; every one of those streams is already pinned to the classic locale by detail/classic_stream.hpp, so what remains here is performance only. The worst cases: XDMF ASCII DataItems (a std::string and a dtype switch per scalar), VTU ASCII arrays (push_back with no reserve, then a second pass with a dtype switch per element) and mdpa.cpp, which materializes one std::string per line of the whole file before parsing. M
• Hoist the per-element dtype switch. detail::dispatch_dtype (detail/value_io.hpp) exists to move a DType switch out of a hot loop, and is used in 14 source files against nearly 500 per-element read_double/read_int/read_point call sites. In I/O the hot ones are the VTU ASCII writer and vtu_to_int64 (detail/vtk_xml.cpp) and the Exodus reader's index shift, coordinate transpose and column_stack (exodus.cpp), which are also fully serial — the exact treatment MED already had, which took it from 0.2–0.6× to parity. S each
• Parallel row formatting in ASCII writers. abaqus.cpp and ansysinp.cpp format rows in parallel into one string per row and stream them in order, with byte-identical output; vtu, vtp, vtk, medit and tecplot still format serially with a locale-aware ostream << int per index, and the OpenFOAM writer formats every coordinate through ostream << std::setprecision(16), the slowest formatting route in the standard library. The reference files are the gate, not the claim — setprecision(16)'s exact output has to be reproduced. S per writer

Binary I/O (VTU).

• b64decode is serial and pushes one byte at a time (detail/vtu_binary.cpp), while b64encode beside it is parallel; decode is the VTU read hot path. Its inverse table is guarded by a hand-rolled static bool init, which is a data race as soon as decoding is parallel or the GIL is released — replace it with a magic static in the same change. A branchless resize-then-index loop first, then parallel chunks after a whitespace pre-scan (chunks are not independent without one). S
• The VTU binary read copies each payload five or six times: vtu_strip copies the base64 text and then substrs it, vtu_parse_binary decodes into a std::vector and memcpys that into the array, vtu_decode_uncompressed builds another vector to strip the header, and the codec returns one std::vector per 32 KiB block which is then copied again (on write, the blocks are concatenated with no reserve). Decode straight into the destination array; independent <DataArray>s can then decode in parallel. It roughly halves peak memory too. M
• Raw <AppendedData> is not written. Both engines read it in .vtu since v16.6.0 (with base64 appended data and BigEndian files); no writer emits it. It is the VTU encoding with no base64 at all, so files meshio++ writes could skip both items above; a format-reach item as much as a performance one. M

Memory and allocation.

• Small, mechanical: five readers slurp the file by hand instead of going through detail/file_source.hpp (ply, medit, ansys, wkt, stl), and several accumulate without reserve although the count is in the header (gid_read.cpp, openfoam.cpp, ansysinp.cpp, unv.cpp). S
• Per-cell heap allocations in every polyhedral reader. AddPolygonBlock/AddPolyhedronBlock take nested std::vectors, so the CGNS, EnSight, MED, OpenFOAM, UNV, FLAC3D and VTU readers allocate one vector per cell and one per face; the VTU path additionally sorts each cell's nodes and buckets cells through a std::map. A CSR ingestion overload — flat node ids plus face and cell offsets, the shape the WASM binding already crosses with — removes all of it at once. It changes the installed API, so it pairs with an ABI bump. M–L

Serial phases inside parallel operations.

• One shared, deterministic facet and edge table. Seven places fill facet or edge keys in parallel and then deduplicate them through a single-threaded unordered_map (surface.cpp, smooth.cpp, partition.cpp, convert_cells.cpp, refine.cpp, detail/marching.cpp, detail/surface_distance.cpp) — on a 10M-tet mesh, ~40M hash operations on one core while the rest idle; it is the dominant phase of extract_surface, extract_skin, smooth, partition, refine and elevate. Nine independent key and hash types exist for it, and detail/face_mesh.hpp has only four consumers. refine.cpp already states the property a parallel version needs — numbering is a pure function of (block, cell, slot) — so a parallel sort of (key, slot) followed by a segmented first-occurrence scan reproduces it exactly there; the other six sites need that argument made before the claim. parallel.hpp has only parallel_for, so the first step is a parallel_sort/parallel_reduce primitive; the second is caching the table so an N-step pipeline stops rebuilding it N times. L
• optimize_volume rebuilds everything on every sweep — a fresh Mesh with new point and connectivity arrays, smooth's node adjacency and boundary hash, and three single-threaded unordered_maps for the 2-3 and 3-2 flips — with one parallel_for in the whole file. Hoist the adjacency out of the sweep and maintain the face and edge maps across flips. S–M
• Welding is single-threaded. clean keeps its own point grid rather than detail::SpatialGrid (its header says so), reads each coordinate through the dtype switch, and deduplicates cells with a std::string built per cell as the hash key, alongside three per-cell vectors; merge scans 27 buckets per point on one core. Move to SpatialGrid, fixed-size keys and a chunked scan merged keep-first. M
• Distance-kernel construction is serial — the triangle soup, edge table, pseudonormals and grid insertion (detail/surface_distance.cpp) — while its queries are parallel; compute_sdf, shrinkwrap, voxelize, compute_curvature and repair all pay it. A BVH is deliberately not the fix: detail/surface_distance.hpp rejects one because its visit order would make the equidistant-triangle tiebreak observable. M
• reorder sorts an index vector through an indirect comparator after a serial bounding-box scan, and partition's space-filling-curve path does the same; a radix sort on the 64-bit curve key replaces both. S
• Hoist the dtype switch in operations, as in I/O: refine (29 call sites), interpolate (21), convert_cells (13), repair and diff (7 each), partition (6). S each
• Operations with no parallel phase at all: remesh (serial clustering), remesh_volume (serial cut), agglomerate (a std::set and an unordered_map allocated per seed), split, undo_green, and hessian (two gradient passes plus deep copies). compute_quality (chunked partials merged in chunk order) and sobolev_deform (a gather-form sparse product) are the in-repo models. L in total
• Radius and k-nearest search in the core. proximity_graph is numpy-only — its own docstring measures 200k points at 6.4 s for a radius graph and 23 s for k=16 on one core — while detail::SpatialGrid already has the needed primitives. Expose the neighbour search as a core operation the Python layer calls, keeping graph assembly in Python per Non-goals. M

Boundaries and startup.

• The Python bindings never release the GIL. Nothing in bindings/python uses gil_scoped_release or a call_guard, so a multi-gigabyte read or a 10M-cell operation blocks every other Python thread, and no caller can convert files in a thread pool. Release after the numpy→Mesh conversion and re-acquire before the Mesh→numpy one, once the b64decode table race above is fixed, auditing that no released region touches a Python object. S–M
• import meshioplusplus takes ~180–210 ms, ~160–190 ms of it the CLI. The package __init__ imports _cli, which imports every verb module, and _common.py imports rich for library use; every MCP process and CLI call pays it. Load _cli lazily through a module __getattr__ and move rich into the CLI. The heavy optional dependencies (h5py, netCDF4, torch, pxr, vtk) are already imported inside functions. S
• A declined C++ read costs a full parse before the Python one starts. vtu_read.cpp builds the whole XML DOM, base64 bodies included, before rejecting lzma; gmsh.cpp rejects $Periodic only after $Nodes and $Elements are parsed; with an ambiguous extension one .msh can be parsed up to six times. A cheap pre-flight — attribute and section-header scans — before the expensive parse, with only early rejections falling back. Every decline is now logged at DEBUG (and raises under MESHIOPLUSPLUS_STRICT_CORE=1), which is what shows how often it happens; the durable fix is §5's core parity. S–M
• The MCP server re-reads the input file on every tool call — 65 call sites in mcp/_tools.py go through an uncached _load(), so an agent's info → clean → decimate → convert parses one file four times. A bounded cache keyed on (path, mtime_ns, size); anything weaker manufactures a stale-read bug. S
• pybind11 per-call overheads, together: the Mesh→numpy conversion re-imports the Mesh class on every call, the contiguity check does a Python attribute lookup per array, operations clone connectivity they never change (transform), and polygon/polyhedron blocks cross the boundary one node id and one face at a time where the WASM binding already uses a CSR triple. S–M
• Flat-binding accessor copies, together: R copies the points twice and shifts connectivity to 1-based with a scalar loop, and Julia's safe accessors copy the borrowed view; both are documented, and both are fixable behind the same accessor names. S
• The browser viewer round-trips every operation result through a VTP file — written into MEMFS, copied out with .slice(), then parsed again by vtk.js on the main thread. A typed-array mesh channel between worker and renderer removes three full passes. L

Deliberately not, and recorded so it is not re-proposed:

• Explicit SIMD intrinsics or -march flags — portability across wheels, WASM and the release binaries is worth more than the scalar kernels cost; revisit only if the harness shows a kernel dominating.
• Kokkos device execution — every parallel_for body captures host pointers (parallel.hpp); the GPU route is the DLPack/CuPy handoff (GPU handoff).
• A BVH in place of the uniform grid — the tiebreak argument above.
• Tuning the pure-Python fallback readers — the fix is making the C++ path accept the file (§5), not a faster fallback; _decimate.py's heap-based twin is deleted once decimate accepts its inputs, not optimised.

Recommended posture: the harness and the measured regressions first; then the small isolated wins — b64decode, the lazy CLI import, optimize_volume — then the shared facet table, the largest total win; the boundary items as their consumers ask.

---

5. Core parity across surfaces

Admission: something the Python layer can do that the C++ core cannot, or that the core can do and a binding cannot reach. A construct that forces the Python fallback is not "slower from C" — it is unreadable from C, Fortran, Julia, R, WASM and the native CLI, none of which has a fallback. Ordered by this project's own consumers, Kratos first.

• MDPA beyond mesh-level blocks. The C++ core reads and writes Nodes/Elements/Conditions/SubModelParts, but Begin Table, Begin Geometries, Begin Mesh <id>, Begin Constraints and non-numeric ModelPartData throw — or, under a lenient read, are skipped and listed in MdpaInfo, which no flat binding exposes (MDPA). M
• Gmsh $Periodic and format 4.0 in the C++ core. $Periodic, both directions: a periodic 4.1 file is unreadable from every flat binding today (Gmsh). The C++ reader also accepts only versions 2.2 and 4.1 (gmsh.cpp), so a 4.0 file, which the Python reader reads, is unreadable from them too. Pairs with periodic node matching in §6. S–M
• VTK-family constructs the C++ readers refuse: <AppendedData> and multiple <Piece>s in .vtp, .vts, .vtr and .vti (the .vtu reader handles both since v16.6.0; its appended-data decoder in vtu_read.cpp is the one to share), and legacy .vtk structured points, structured grid and rectilinear grid (VTU, VTK). S–M
• XDMF 2 and XPath references — the C++ core implements XDMF 3 only, and Reference="XML" DataItems not at all (XDMF). M
• MED multi-mesh files and profiles, which are Python-only and not reachable even under a lenient C++ read (MED). M
• Netgen extras — periodic identifications, materials/bcnames/cd2names/cd3names, edgesegmentsgi2 and the .vol.gz container (Netgen). S–M
• An Exodus writer that carries sets and steps. It writes element blocks but no node sets or side sets, so only element-block regions round-trip, and it writes one step per file; a multi-step writer is a stateful object of the XdmfTimeSeriesWriter shape (Exodus). M
• Sets → regions, phase 2. XDMF Sets are not mapped (UNV groups became regions in v15.6.0 and Ansys components in v16.3.0), and VTU/VTP need a documented region convention (Named regions). M
• Named Side regions surviving operations. subdivide, agglomerate, undo_green and convert_cells(simplexify) drop them through their parent-cell remap, and the cutters (slice, isosurface, extract_surface, extract_skin) drop them outright (Named regions); for a Kratos model, Side regions are where the boundary conditions live. M
• A structured pipeline report on the flat ABI. C, Fortran, Julia and R receive status plus mio_last_error() only; a caller-buffer JSON accessor is recorded as a follow-up, as are the v2 multi-mesh steps (Inputs: for Merge/Interpolate/UndoGreen, an Output.Pattern for Split/partition) (pipelines). S–M
• PCD binary_compressed on the flat ABI. The C++ API and Python write it (write_pcd(..., PcdData::BinaryCompressed), data="binary_compressed"), but WriteOptions/mio_write_opts only carry the VTK block codecs, so C, Fortran, Julia, R and WASM can read it and cannot write it; the fix is an appended MIO_CODEC_LZF enumerator routed to write_pcd (PCD). S
• glTF write options on the flat ABI. mio_write("x.glb") and writeMesh write the defaults; the colour field, colormap and range, split angle, up axis and unit scale are reachable from Python, C++, both CLIs and MCP only, so C, Fortran, Julia, R and WASM cannot colour a .glb or change its axis (glTF). The fix is a mio_gltf_opts struct with a reserved tail, following the mio_curvature_opts shape. S
• Point/cell sets beyond regions in the core, so the convert -s/-d sets↔data conversions work in the native CLI and flat bindings (Julia). S–M

---

6. Operations

Admission: a new operation, or a public face for machinery that already exists privately inside one.

Generation. Almost every operation transforms a mesh you already have. The exceptions all start from something else — grid from a lattice (detail/grid_lattice.hpp), voxelize/compute_sdf from a surface's bounding box, remesh_volume from a closed surface — and nothing builds a shape from parameters, sweeps one, or triangulates a domain.

• Primitive constructors — box, sphere, cylinder, disk, in their own operations/primitives.hpp beside grid. Dependency-free, and it removes the fixture-file dependency from tests, docs, notebooks, the browser demo and the MCP server; it is also the canonical mesh the §3 conformance matrix needs. Highest leverage per line of code in this document. S
• extrude — 2-D → 3-D sweep (triangle → wedge, quad → hexahedron) with nlayers and per-layer offsets, carrying regions to side and cap regions. The most-requested generation primitive. M
• revolve — extrude's rotational sibling around an axis, sharing its layer machinery; degenerate cells on the axis are the only new work. M
• Delaunay / constrained 2-D meshing — genuinely useful, but robust geometric predicates are where dependency-free stops paying. Better as an optional Triangle or Gmsh backend, off by default, following the KaHIP pattern. L

Analysis and editing.

• Feature edges as a line mesh. The feature-angle crease test exists three times (decimate, smooth, remesh each carry mFeatureAngleDeg) and only ever pins nodes; one public op emitting line cells serves inspection, boundary-condition picking and those three in one place. S
• Hausdorff distance between two meshes — a symmetric max-reduction over the shipped distance_to_surface, returning the scalar the remesh/decimate tests and the conformance matrix want to assert on. S
• Periodic node-pair matching — given two boundary regions and a transform, return the matched node pairs. $Periodic already round-trips as metadata and proximity_graph already does minimum-image search; Kratos periodic conditions are the consumer. S–M
• A quality gate — pass/fail thresholds over the metrics and histograms compute_quality already produces, a check CLI verb that exits non-zero, and uniform --json output (today only a handful of verbs take it — not info, quality, diff or convert). What a CI pipeline over meshes actually scripts. S–M
• Region set algebra — union, intersection, difference, rename and retag of named regions. Regions are a first-class layer with add/replace/enumerate only; every consumer that builds boundary conditions from them hand-rolls this. S
• Time-axis resampling of sequences. The sequence engine reads time values and drives N→N, fan-in and fan-out, but never resamples; aligning two solvers' timelines (or a solver and a surrogate's) is the missing step before a pairwise diff or a training pair. S–M
• agglomerate follow-ups — coplanar boundary-face merging (fusing adjacent group faces on one plane into a single polygon) and a shape-quality absorption gate, both deferred when it shipped (agglomerate). S–M

---

7. Ecosystem reach

Admission: getting what exists to the people who would use it.

• Registration and distribution. All pure logistics and all blocking real adoption; the effort is small per target but the calendar time is set by each registry's review, so start early. S each, calendar-bound

Target	Status today
conda-forge	Not started.
CRAN	R CMD check --as-cran already runs in CI, so the gate exists; only submission remains (R).
Julia General	Not registered, so Pkg.add("MeshioPlusPlus") does not work yet (Julia).
Blender Extensions Platform	Extension zips are built and attached to every release; only the listing remains (Blender).
ConanCenter and the vcpkg registry	Recipes are self-hosted and CI-validated under packages/; neither is submitted (C API).
ParaView	The plugin exists; installing it with the wheel is the §3 item.

• Rust bindings over the C API — the next language by scientific adoption after Julia and R, and the ABI/SOVERSION work makes it cheap. M
• Interop phase 2 — the Open3D and DOLFINx bridges exist as named stubs that raise NotImplementedError (interop); the pinned-memory staging for CuPy is wired in C++ (v8.5.0) but not from Python (GPU handoff). S–M each
• A CLI chatbot — meshioplusplus chat, a natural-language entry point from the terminal. The MCP server already exposes the whole Python surface to an agent, but only to a host that speaks MCP (Claude Desktop, an IDE); a user with an LLM API key cannot ask the CLI "why does this file fail to convert" or "clean this mesh and tell me what changed". A thin client over existing work, not new mesh functionality. S–M
  • Reuse, not reimplementation. The tool-calling loop drives mcp/_tools.py's TOOL_REGISTRY in-process — never a second copy of its dispatch, path sandbox or result sanitizing — so every tool added there reaches the chatbot for free, under the parity guarantee test_every_operation_has_a_tool already gives MCP. The pattern is proven, not hypothetical: mcp/_http.py already dispatches the dashboard's /api/tools/<name> through the same registry.
  • v1 scope. A stateless-per-turn loop (history + tool schemas out, tool calls executed, results fed back until a plain-text reply) over files under the CWD or --root, reusing the MCP sandbox; context lives only in the terminal session — no persistence, no server.
  • Packaging. A chat extra carrying the model provider's SDK, kept out of [all] and imported only by the new _cli/_chat.py, so _tools.py, the default CI matrix and every other surface are untouched; without the extra or an API key the verb fails by name, following the meshioplusplus-mcp entry-point precedent. Provider-agnostic wiring and persisted history are follow-ups, not v1.
  • Docs. A doc/chat.md page modelled on MCP, plus the README and tool-table updates the AGENTS.md checklist requires.

---

8. Long run (spike first)

Admission: work whose shape is unknown until an investigation writes it down. Findings before code.

Scale. Memory-mapped reads (mmap, via ReadOptions and the C ABI but not Python's read()) roughly halve the peak footprint of a large read, and the XDMF series appender and Python's chunked write_dataset write a series or a dataset without holding it — but nothing writes one mesh larger than memory, and no operation streams. Run the §3 benchmark tier first; it decides whether either item below matters.

• Streaming / chunked writes of one mesh, the counterpart to selective and memory-mapped reads. L
• Out-of-core operations for the ops that are already block-local. XL

NURBS and higher-order geometry. The data model is strictly linear/Lagrange polytopes: a CellBlock is a cell-type string plus a node-index array. NURBS is a genuinely different object — control points, weights, knot vectors, and a parametric mapping — and CAD/IGA formats (STEP, IGES, Rhino 3dm, .iga) express geometry that no current cell type can hold. This is the most architecturally invasive item on the list and should be approached as a research spike, not a feature.

• Spike: how far can the current model stretch? Higher-order Lagrange cells already exist (hexahedron27, VTK-Lagrange types); a rational Bézier/NURBS patch needs weights and a knot vector, which have nowhere to live. Determine whether a side-channel struct (the MedInfo/GmshInfo precedent) suffices, or whether the Mesh needs a genuine second entity kind. Write the finding up before committing. M
• Read-only CAD ingestion first: a NURBS surface tessellated to a triangle mesh at a requested tolerance, with the parametric data carried out-of-band. This delivers most of the practical value (getting CAD into the mesh world) without touching the data model, and is the natural first release. L
• A real IGA data model — patches, control nets, weights, knots, trimming curves — plus formats and evaluation. This is XL, likely a separate library or a major version, and should only be attempted if the spike shows real demand.
• Dependency reality: robust STEP/IGES parsing effectively means OpenCASCADE, which is a heavyweight LGPL dependency. If ingestion goes ahead, it must follow the KaHIP/Polyscope pattern — strictly optional, off by default, never in the core, licence implications documented. Findings before code.
• No C++ tessellator exists yet — tessellate is Python-only — so read-only CAD ingestion reaching every surface also carries that port.

Recommended posture: spike and document; do not schedule implementation until the spike says what shape it takes.

---

Non-goals and decisions taken

Recorded so they are not re-proposed as gaps.

• MPI in the library — none planned (C++ API); partition's ghost layers produce the halo an MPI assembly in the owning application needs.
• Solver-coupled physics-ML — assembled solver residuals, adjoints and Sobolev training, co-simulation, active-learning labeling, adaptive remeshing driven by a surrogate, and MPI model-part gathering. Every one needs a live solver (an assembly routine, its tangent, its communicator) and meshio++ has no notion of a discrete system; they belong in the application that owns the solver, the division Symbolic and physics describes.
• The Python-only layers stay Python — pmsh, zarr, cae and usd, and the physics-ML surface (tessellate, grids, point budgets, proximity graphs, datasets, training) are export targets and tooling for a training pipeline, registered in Python rather than the shared C++ registry. A core kernel they call (the neighbour search in §4) does not change that.
• Polyscope in the release CLI binaries — excluded deliberately, so view/screenshot there report the build flag rather than opening a window (viewer).
• General meshing algorithms as operations — hex and hex-dominant meshing, boolean/CSG, boundary-layer inflation, quadrangulation and geodesic distance are each a library in their own right; the answer is an optional backend (the KaHIP pattern), not a native implementation.
• KaHIP in the WASM build — no Emscripten port, and a graph partitioner would bloat every consumer's bundle; "auto" resolving to the SFC method is the answer, and "kahip" throwing by name is the contract.
• zstd/lz4 codecs, memory mapping and the Kokkos backend under Emscripten — the codecs have no Emscripten port and are compiled out, CMake refuses the Kokkos backend under Emscripten, and there is nothing to map inside MEMFS.
• Single-file output for DOLFIN, TetGen and EnSight from WASM, or lifting DOLFIN's simplicial restriction — both are facts of the formats; writeMesh's written-paths return value (v11.2.0) is the fix for the bookkeeping they cause.
• Polyhedron blocks as a WASM gap — vtu, ensight, cgns, med and openfoam all write them from WASM; a consumer that cannot is constrained by its own data model.

---

Suggested sequencing

Open work only; what shipped is in CHANGELOG.md.

1. Format reach (§1) — the remaining halves as users ask.
2. Spack package upkeep (§2) — a small, mechanical catch-up (recipes, variants, one release-checklist line) that unblocks HPC users on the current release; then a checklist step, so it stays current.
3. Sanitizer leg, then fuzzing (§3) — a parallel track from day one; it does not compete for the same attention as features.
4. Performance (§4) — the harness and the measured regressions first, then the isolated wins (b64decode, the lazy CLI import, optimize_volume), then the shared facet table.
5. Primitive constructors (§6) — a few days, and a prerequisite of the conformance matrix, every demo surface and extrude/revolve.
6. Core parity (§5) — MDPA first, then Side-region survival and sets → regions, then the rest by consumer demand.
7. Registration (§7) — calendar-bound, so start the submissions early and let them run alongside everything else.
8. Long-run spikes (§8) — the benchmark tier decides the scale items; the NURBS spike is scheduled independently of the rest.