meshio++

There are various mesh formats available for representing unstructured meshes. meshio++ can read and write all of the following and smoothly converts between them:

Native parity additions: XDMF2/3 and absolute DataItem references; Netgen name/periodic tables, two-line edges and zlib-enabled .vol.gz; PCD compressed writes via MIO_CODEC_LZF / codec="lzf"; and parameterized glTF export in C, Fortran, Julia, R and WASM. See formats, C API and WASM for options and indexing conventions.

Region-backed point/cell sets now convert to and from scalar integer data in the native core, C/Fortran/Julia/R/WASM bindings and native CLI convert -s/-d. Python's mutating Mesh methods preserve their ordering and metadata contract; pipelines expose SetsToData / DataToSets, and MCP exposes sets_data. See sets/data conversions.

Native C/Fortran/Julia/R pipelines also expose owning structured JSON reports (steps and warnings), without replacing their existing status-only entry points; see pipeline reports.

Exodus writers now preserve node/side sets and expose fixed-grid time-series writers across Python, C++, C, Fortran, Julia, R and WASM; sequence fan-in can write one Exodus file. Pipeline schema Version 2 adds file-backed Merge/Interpolate/UndoGreen and terminal Split/Partition with Output.Pattern. Version 1 remains unchanged. See Exodus and pipeline v2.

XDMF time-series writers preserve fixed point, cell and side regions alongside the shared topology. Python's array-based writer accepts regions= and its series reader exposes reader.regions; native writers use the regions already on the mesh. See XDMF time series.

Serial VTK XML readers (.vtu, .vtp, .vts, .vtr, .vti) accept multiple pieces and appended raw/base64 arrays in both the native core and Python reference, with UInt32/UInt64 headers and either byte order. Pieces concatenate without welding; VTP/VTS/VTR/VTI writers remain inline. Native legacy .vtk reads also support structured points, structured grids and rectilinear grids in ASCII and big-endian binary; see the format documentation.

Abaqus (.inp), Abaqus results file (.fil, ASCII and binary; every increment is a step, nodal and element results by output location; written too, the model and an increment per step), ANSYS Fluent mesh (.msh; cells rebuilt from faces, written face-based as Fluent reads it, zones as named regions), Ansys MAPDL coded database (.cdb, .inp; degenerate shapes resolved, components as named regions) and results (.rst, .rth; nodal solution, reactions and element nodal stresses, strains and forces, every result set is a step; distributed solves merged, the full rotor of a static cyclic model), AVS-UCD (.avs), CAE sample layout (.npz, what PhysicsNeMo's DoMINO/Transolver datapipes read), CGNS (.cgns), DOLFIN XML (.xml), DOLFINx VTX (.bp, an ADIOS2 directory, read and written; each step a mesh, with an optional ADIOS2 build or the adios2 package), COMSOL mesh (.mphtxt text and .mphbin binary; Selections as named regions, COMSOL's own node order), Exodus (.e, ), mesh directory (////, text or ElmerGrid's binary; bodies and boundaries as named regions, ElmerGrid partitions merged, written as , recognised by content), (geometry, /; C and Fortran binary), input (, spec 2.5/3.0/4.0 read, 4.0 written; sets and surfaces as named regions, MeshData read and written) and plot files (; every state is a step, remeshed runs included, surface variables on facet blocks), neutral file (, 4.41 to 2020.1 read, 8.2 mesh written; properties and groups as named regions, every output set is a step, a series written as one output set per step), (, named cell groups), (mesh , field ), (), (/, , ; writing via a vendored gidpost, reading is meshio++'s own code), (), native mesh (; / as named regions, Code_Aster's own node order), results (, ASCII or binary written too; every increment is a step, tensors keep the file order, opt-in von Mises and principal values), (, experimental, meshio++-specific), (, , read-only, through OpenCASCADE; triangulated, with its NURBS B-rep), (, ), (), (; tables, geometries, blocks and , including nested sub-model-part geometry/constraint membership, kept in a side channel on every surface; the Kratos IGA geometry beside a deck loaded with it), (; trimmed NURBS faces, trimming loops and coupling edges for Kratos IGA, sampled as Lagrange cells and kept exactly in a side channel on every surface, and the dialect), (, , also gzip/bzip2; active cells of refined meshes, subdomains, side, edge, shell-face and node sets as named regions, written back in the 1.8.0 layout), (keyword input , , ; and as named regions, followed), (the binary state database, and its ... family, read and write; every state is a step, parts are regions, deletion flags a mask; SPH, airbag particles, rigid bodies and roads), (the LSDA binary output; outputs are steps, the other databases field data; written too, a cloud per step with each database on its own cadence), (, ) and its solutions (, ; size maps and metrics for MMG, vertex and element fields attached to their mesh), mesh (, recognised by content; order-2 curved meshes in MFEM's own numbering, attribute sets as named regions, NURBS meshes with their patches kept and written back) and grid functions (), (; native named multi-mesh access and nodal/element profile expansion), (mesh , field ), and (bulk data, , , ; HyperMesh components and OptiStruct s as named regions), results (; every subcase, mode, time or frequency is a step; coordinate systems applied; corner, ply, station and grid point force values; written too, with nodal results and the centre stress and strain of shells, solids and rods, one step per call or a series), results (, MSC and NX, 32- and 64-bit; every subcase, mode, time or frequency is a step; coordinate systems applied; corner, ply, station and grid point force values), (, ), , (format versions 2.2, 4.0, and 4.1, ; native periodic reads across language bindings via format metadata, periodic writes in 2.2/4.1), (read as a general glTF 2.0 reader; written as the surface with per-vertex normals, fields as attributes, into ) (, ), (), (), (), (, , ), starter deck (; parts, subsets, groups, boxes and surfaces as named regions, units applied, followed; written with optional placeholder materials and properties), animation files (…; one state per file, a run's files a transient sequence) and time-history files (; each output a step of field data; written too, from the field data), input deck (, written by name) and formatted post file (; sets as named regions, every increment is a step), (, ; named components as regions), (Point Cloud Library point clouds, //) (), (, , , ), (, memory-mapped), (), (; a time-indexed collection of serial or parallel VTK XML files, read by ParaView), (; a parallel index over per-piece files) and (; a parallel index over per-piece files, ghost cells and output included), (a closed surface written as an AP214 solid that CAD tools open, one face per plane, cylinder, cone, sphere or torus; read through OpenCASCADE as a triangulation plus its NURBS B-rep) (, ), (), (, , and binary ; FE, ordered, polygonal and polyhedral zones; through a TecIO you supply), , , (2D drawing: 2D direct, 3D via skin projection; reads back only what meshio++ wrote, as planar) (), (LaTeX 2D drawing: 2D direct, 3D via skin projection; reads back only what meshio++ wrote, as planar) (), (), (), (VTK XML ImageData; a regular lattice) (), (), (; Kitware's HDF5-based VTK format: time, partitions and fields in one file, read by ParaView), (; a VTK XML multiblock index, one piece per cell block), (), (, VTK XML RectilinearGrid), (, VTK XML StructuredGrid), (; / writes raw appended binary, VTK's own default layout, about a quarter smaller than inline binary with no base64), () (), (, ), XYZ point clouds (headerless ASCII: , , , , , ), (, recognised by its file name, with its displacements and stresses), (, the PhysicsNeMo mesh layout).

Files only their vendor's software can read — Abaqus .odb, MSC Marc .t16, Femap .modfem, and Ansys results beyond the native .rst reader — have vendor-runtime routes: scripts under contrib/ that run in the vendor's own Python and write a series meshio++ reads.

One C++ core, six language surfaces, and the tools built on them; the architecture page is the map, with every box linked to the page that owns it.

Release v16.32.0 closes the performance roadmap: the browser viewer renders from typed arrays instead of a VTP file, more readers and operations read through shared cursors and reuse scratch, import meshioplusplus loads formats lazily, and transform works on ragged meshes. C++ ABI 23 is a bump (DistanceQuery layout); see the changelog for details.

meshio++ ships a C++20 core (built with pybind11 + scikit-build-core) that reads and writes most formats with zero-copy numpy at the I/O boundary, plus optional HDF5/netCDF acceleration and a selectable parallel backend (AUTO by default — prefers OpenMP, then STL+TBB, then sequential; override with -DMESHIOPLUSPLUS_PARALLEL_BACKEND=..., including a bring-your-own Kokkos host backend). Every format has a pure-Python fallback, so behaviour and file compatibility are identical whether or not the native libraries are present. For a standalone C++ build use build/configure.sh (Linux/macOS) or build/configure.bat (Windows). Full docs (install, data model, per-format options, CLI) live at the documentation site (sources under doc/).

Install with

pip install meshioplusplus[all]

([all] adds the h5py and netCDF4 runtime extras and the polyscope desktop viewer. By default, meshio++ uses numpy and rich.) The wheels carry the compiled core with OpenCASCADE linked in, so CAD output (fit_cad, to-cad) and STEP/IGES reading work out of the box. The wheels are licensed MIT AND (LGPL-2.1-only WITH OCCT-exception-1.0), because that kernel is LGPL-2.1 with the OCCT exception; see LICENSE-BINARIES.md and THIRD_PARTY_NOTICES.md. You can then use the command-line tool

meshioplusplus convert    input.msh output.vtk   # convert between two formats

meshioplusplus info       input.xdmf             # show some info about the mesh

meshioplusplus compress   input.vtu              # compress the mesh file
meshioplusplus decompress input.vtu              # decompress the mesh file

meshioplusplus binary     input.msh              # convert to binary format
meshioplusplus ascii      input.msh              # convert to ASCII format

meshioplusplus merge      a.vtu b.vtu out.vtu    # merge meshes (optional --weld)

meshioplusplus transform  in.vtu out.vtu --translate 1,2,3   # affine transform
meshioplusplus mirror     half.vtu whole.vtu --keep-original  # reflect, cells right side out
meshioplusplus extrude    profile.vtu slab.vtu --num-layers 4  # sweep a surface into a volume
meshioplusplus clean      in.vtu out.vtu --weld              # weld / prune / de-dup
meshioplusplus crop       in.vtu out.vtu --bbox 0,0,0,1,1,1  # subset by region
meshioplusplus crop       in.vtu out.vtu --where "d<0"       # ... or by a data predicate
meshioplusplus split      in.vtu 'out_{key}.vtu' --by type   # split by criterion
meshioplusplus stats      mesh.vtu                           # geometric statistics
meshioplusplus convert-cells in.msh out.vtu --mode simplexify  # hexes -> tetra
meshioplusplus subdivide  in.vtu out.vtu                      # polyhedral refinement, any 3D cell
meshioplusplus refine     in.vtu out.vtu --levels 2          # uniform subdivision
meshioplusplus refine     in.vtu out.vtu --where "q<0.3"     # adaptive, closed conformingly
meshioplusplus partition  in.vtu 'out_{part}.vtu' --nparts 4 # N balanced parts
meshioplusplus remesh-volume in.vtu out.vtu --cell-size 0.5  # retetrahedralize (isosurface stuffing)
meshioplusplus optimize-volume in.vtu out.vtu               # ODT-remesh a tet mesh (relocate + flips)
meshioplusplus smooth     in.vtu out.vtu --iterations 20     # relax node positions
meshioplusplus smooth     in.vtu out.vtu --method odt        # ODT smoothing, tet-only
meshioplusplus interpolate src.vtu tgt.vtu out.vtu           # transfer fields across meshes
meshioplusplus slice      in.vtu out.vtu --normal 0,0,1      # planar cross-section
meshioplusplus curvature  in.vtu out.vtu                     # per-vertex mean/Gaussian curvature
meshioplusplus normals    in.vtu out.vtu --split-angle 30    # point/cell normals, split at creases
meshioplusplus feature-edges in.stl edges.vtp --angle 45   # sharp/open/non-manifold edges as lines
meshioplusplus segment    in.stl patches.vtu --method primitive  # one patch per plane/cylinder/... as regions
meshioplusplus fit-primitives patches.vtu fits.vtu --json  # plane/sphere/cylinder/cone/torus per patch
meshioplusplus to-cad part.stl part.step                 # one STEP face per surface (any build; OpenCASCADE trims and checks)
meshioplusplus convert    in.stl part.step               # a closed skin as a CAD solid, coplanar facets merged
meshioplusplus convert    part.step part.cad.json        # a CAD model as Kratos IGA geometry (through OpenCASCADE)
meshioplusplus hausdorff  a.stl b.stl --max 1e-3             # worst-case distance; nonzero exit above it
meshioplusplus check      part.vtu -r "scaled_jacobian >= 0.2"  # quality gate: exit 1 on failure, 2 if unchecked
meshioplusplus regions    in.msh out.msh --union inlet=in_a,in_b --rename "Surface 3=Outlet"  # region algebra
meshioplusplus periodic   in.msh --slave inlet --master outlet --translate 2,0,0  # periodic node pairs
meshioplusplus repair     in.vtu out.vtu                     # orientation, holes, bowties
meshioplusplus shrinkwrap in.vtu scan.stl out.vtu           # project onto a target surface
meshioplusplus sobolev-deform in.vtu out.vtu --array d --length-scale 0.5  # filter a displacement
meshioplusplus isosurface in.vtu out.vtu --array T --values 350  # level set of a field
meshioplusplus sdf        skin.stl field.vti --resolution 128,128,128  # signed distance field
meshioplusplus data gradient in.vtu out.vtu --array T           # grad / div / curl of a field
meshioplusplus data hessian  in.vtu out.vtu --array T           # second derivative, gradient's companion
meshioplusplus data estimate-error in.vtu out.vtu --array T --marking dorfler --marking-value 0.6  # ZZ error indicator + marking
meshioplusplus data integrate in.vtu --array density            # total / mean, per region

meshioplusplus data info  mesh.vtu                           # summarize data arrays
meshioplusplus data calc  in.vtu out.vtu --point "s = norm(v)"   # derive a field
meshioplusplus data to-cell  in.vtu out.vtu --keys T         # point -> cell average
meshioplusplus data normalize in.vtu out.vtu --cell damage --to 0,1

meshioplusplus pipeline   settings.json                      # run a whole declarative
                                                             # read -> ops -> write chain

with any of the supported formats.

The same verbs are available as a standalone C++ binary that needs no Python: grab a ready-to-run, statically-linked build for Linux/macOS/Windows from the GitHub Releases page, or build it yourself with build/configure.sh --cli --build (or -DMESHIOPLUSPLUS_BUILD_CLI=ON). It links only the C++ core. Named regions — and so point/cell sets — are carried there since v8.1.0, so info lists them and diff compares them; convert -s/-d is still Python-only.

A whole chain of operations can be described declaratively in one settings.json (Input → Operations → Output) and run with meshioplusplus pipeline settings.json — or meshioplusplus.run_pipeline(...) from Python, and the same engine from C/Fortran/Julia/R/WASM. See the settings pipeline.

In Python, simply do

import meshioplusplus

mesh = meshioplusplus.read(
    filename,  # string, os.PathLike, or a buffer/open file
    # file_format="stl",  # optional if filename is a path; inferred from extension
    # see meshioplusplus convert --help for all possible formats
)
# mesh.points, mesh.cells, mesh.cells_dict, ...

# mesh.vtk.read() is also possible

to read a mesh. To write, do

import meshioplusplus

# two triangles and one quad
points = [
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 1.0],
    [1.0, 1.0],
    [2.0, 0.0],
    [2.0, 1.0],
]
cells = [
    ("triangle", [[0, 1, 2], [1, 3, 2]]),
    ("quad", [[1, 4, 5, 3]]),
]

mesh = meshioplusplus.Mesh(
    points,
    cells,
    # Optionally provide extra data on points, cells, etc.
    point_data={"T": [0.3, -1.2, 0.5, 0.7, 0.0, -3.0]},
    # Each item in cell data must match the cells array
    cell_data={"a": [[0.1, 0.2], [0.4]]},
)
mesh.write(
    "foo.vtk",  # str, os.PathLike, or buffer/open file
    # file_format="vtk",  # optional if first argument is a path; inferred from extension
)

# Alternative with the same options
meshioplusplus.write_points_cells("foo.vtk", points, cells)

For both input and output, you can optionally specify the exact file_format (in case you would like to enforce ASCII over binary VTK, for example).

Skin extraction

meshioplusplus.extract_skin derives the boundary surface of a 3D volume mesh (the Kratos SkinDetectionProcess face-hashing algorithm — faces occurring exactly once are boundary; points are compacted, point_data follows):

vol = meshioplusplus.read("part.msh")     # tetra/hexa/wedge/pyramid mesh
skin = meshioplusplus.extract_skin(vol)   # triangle/quad/... surface mesh

The STL and PLY writers do this automatically for volume meshes (pass skin=False for the legacy drop-volume-cells behavior), and the SVG/TikZ writers render 3D meshes by projecting the skin through an orthographic camera (azimuth/elevation/roll in degrees, default the classic CAD isometric view) with painter's-algorithm depth ordering — that is exactly how the Stanford-bunny logo above is drawn.

Publication-quality vector figures

The SVG and TikZ writers can colour each face by a data array, turning them into figures you can drop straight into a paper — resolution-independent, and with no extra dependency: the colormaps are built into the core.

annotated = meshioplusplus.attach_quality(mesh)
meshioplusplus.write(
    "quality.svg", annotated,
    color_by="quality:scaled_jacobian",   # or any point_data / cell_data name
    cmap="viridis",                       # viridis / coolwarm / turbo / magma / … (and *_r)
    colorbar=True,
)

The bundled bracket coloured by element quality — the same figure tools/gen_doc_images.py regenerates.

Point data colours a face by the mean of its corner values, cell data by its owning cell's value — for a volume mesh, tracked through the extracted skin's parent-cell provenance, so a per-cell material or metric lands on the right facet. Multi-component arrays reduce to a component or to their magnitude; vmin/vmax set the range (default: the drawn faces' finite range), and non-finite values take nan_color. From the command line:

meshioplusplus convert mesh.vtu figure.svg --color-by temperature --colorbar

Colouring is available from Python, from C++ directly, and from both CLIs; the flat C/Fortran/WebAssembly bindings reach these writers through the shared registry and always emit the default styling.

glTF export

meshioplusplus.write("model.glb", mesh) writes the surface of a mesh as glTF 2.0, which three.js, Blender and every web dashboard load directly: the skin of volume cells, 2-D cells, lines and point clouds, one named node per cell region, with normals split at creases and every point field exported raw as an _NAME attribute. color_by also bakes a field through a colormap into COLOR_0 with an unlit material. The output is Y-up and float32; the axis change, unit scale and recentring offset live on the root node rather than in the coordinates. See glTF.

meshioplusplus.write("model.glb", mesh)                                  # defaults
meshioplusplus.write("result.glb", mesh, color_by="temperature", cmap="turbo", scale=0.001)

meshioplusplus convert result.vtu result.glb --color-by temperature --split-angle 45

Surface extraction

meshioplusplus.extract_surface is the general form of skin extraction: it picks the dimension automatically (a volume mesh → boundary faces, a 2D surface mesh → boundary edges) and can record each facet's parent cell id (record_parent_ids=True). See the surface extraction docs (doc/extract_surface.md).

surf = meshioplusplus.extract_surface(vol)                  # faces (or edges for a 2D mesh)
edges = meshioplusplus.extract_surface(sheet, record_parent_ids=True)

Mesh quality

meshioplusplus.compute_quality scores every cell on a set of geometric quality metrics (area/volume, scaled Jacobian, aspect ratio, skewness, interior/dihedral angles, warpage) and flags inverted/degenerate cells; attach_quality writes them back as cell_data. See doc/mesh_quality.md.

report = meshioplusplus.compute_quality(mesh)
print(report["num_inverted"], "inverted cells")
annotated = meshioplusplus.attach_quality(mesh)   # metrics as cell_data

meshioplusplus.check_quality(mesh, "scaled_jacobian >= 0.2; aspect_ratio <= 5 @ 1%") turns the metrics into a pass/fail gate tested per cell, and the check CLI verb exits 1 when it fails and 2 when it could not run — what a CI job over meshes scripts; every report verb of both CLIs also takes --json with one shared, strictly valid shape. See doc/mesh_quality.md and doc/cli.md.

meshioplusplus.check_conformity(mesh) reports the hanging nodes, T-junctions and partly overlapping or unwelded facets of a mesh, flags the hanging points in conformity:hanging, and says which find_interface mode suits it; the check-conformity CLI verb exits 1 when the mesh is not conforming. See doc/check_conformity.md.

Reordering / renumbering

meshioplusplus.reorder renumbers nodes and elements to reduce sparse-matrix bandwidth (Reverse Cuthill–McKee) or improve cache locality (Morton / Hilbert space-filling curves). It is a pure permutation — geometry and all data preserved — and returns the applied node/cell permutations so external arrays can be remapped. compute_bandwidth measures the before/after connectivity bandwidth. See doc/reorder.md.

out = meshioplusplus.reorder(mesh, method="rcm")            # "morton" / "hilbert" too
out, node_perm, cell_perms = meshioplusplus.reorder(mesh, return_permutation=True)
print(meshioplusplus.compute_bandwidth(mesh), "->", meshioplusplus.compute_bandwidth(out))

Comparison (diff)

meshioplusplus.diff compares two meshes and reports whether they are equivalent within a tolerance (abs_err <= atol + rtol*|expected|), with a structured breakdown (points, cells, data, named sets) and an overall verdict (identical / equal within tolerance / different); meshes_equal is the boolean wrapper for test suites. An optional unordered=True mode matches nodes by spatial proximity, so a shuffled node order still compares equal. See doc/diff.md.

assert meshioplusplus.meshes_equal(a, b, atol=1e-8)         # ideal in a regression test
report = meshioplusplus.diff(a, b, unordered=True)          # tolerant to shuffled node order
print(report["verdict"])

The meshioplusplus diff a.vtu b.vtu CLI verb sets a nonzero exit code when meshes differ, for direct use in CI / Makefiles.

Hausdorff distance, feature edges, region algebra and periodic pairs

Four analysis and editing operations (v16.23.0), each on every surface — Python, C++, C, Fortran, Julia, R, WASM, both CLIs and MCP. meshioplusplus.hausdorff_distance(a, b) is how far apart two surfaces are at their worst, the scalar a remeshing or format round trip is asserted on (--max D turns the CLI verb into a CI gate) — doc/hausdorff.md. meshioplusplus.feature_edges(mesh) writes the sharp, open, non-manifold and inconsistently wound edges as a line mesh; it is the one crease test decimate and smooth pin nodes with — doc/feature_edges.md. meshioplusplus.edit_regions(mesh, edits) unions, intersects, subtracts, renames, retags and deletes named regions — doc/regions.md. meshioplusplus.match_periodic_nodes(mesh, slave, master, translate=...) pairs the nodes two boundary regions share under a translation or rotation, the input of a periodic boundary condition — doc/periodic.md.

d = meshioplusplus.hausdorff_distance(before, after, face_samples=4)["distance"]
edges = meshioplusplus.feature_edges(mesh, feature_angle=45)
mesh = meshioplusplus.edit_regions(mesh, [{"op": "union", "inputs": ["a", "b"], "output": "inlet"}])
pairs = meshioplusplus.match_periodic_nodes(mesh, "inlet", "outlet", translate=(2, 0, 0))

Skin to CAD: segment, fit primitives, STEP, CAD faces

The mesh-side half of reverse engineering a triangle skin into a CAD model, with no CAD kernel in the loop, on every surface — Python, C++, C, Fortran, Julia, R, WASM, both CLIs, MCP and the pipeline. meshioplusplus.segment(mesh) splits a surface into smooth patches by region growing, bounded by its feature edges, with a curvature test that separates a fillet from the planes it is tangent to, or, with method="primitive", one patch per plane, cylinder, cone, sphere or torus its facets lie on; the patches come back as segment:patch cell data, a patch-adjacency table and one cell region each. meshioplusplus.fit_primitives(mesh) fits a plane, sphere, cylinder, cone or torus to each patch, deterministically, and keeps the simplest within tolerance — a fit outside it is reported as none with its residual, never rounded to a neat primitive. The step writer hands a closed skin to CAD as an AP214 B-rep solid with no CAD kernel, validated against OpenCASCADE: each fitted patch one analytic face (plane, cylinder, cone, sphere, torus) bounded by lines and circles, each flat region one planar face with its holes, not one face per triangle (faces="facet" keeps that layout). meshioplusplus.fit_cad(mesh, path) runs the chain on a raw skin and builds the solid: without a CAD kernel through its native backend (the analytic faces above), and with OpenCASCADE (included in the PyPI wheels and the release binaries, on by default from source; LGPL-2.1 with the OCCT exception) doing the whole CAD half: each fitted patch becomes an analytic face trimmed by its boundary, a rolling-ball fillet (of one radius or a radius varying along it, also between freeform faces) or a helical thread face (straight or rounded, with or without a cylinder beside it) its own face, and any other patch no primitive explains a freeform B-spline face (split first into pieces a plane or an unrolled cylinder parametrises), sewn into a checked solid and written as STEP a CAD tool can edit (or as Kratos CAD JSON), with anything the kernel cannot build demoted to facets, never silently wrong; since v17.2.0 a chamfer, a conic blend, a tapered or variable-pitch thread and a corner where blends meet (a vertex blend) get their own faces too (every surface, WASM through its full artifact; to-cad in both CLIs). Since v17.11.0 the loop is iterative and tolerance-driven, each stage off by default: meshioplusplus.refine_patches(mesh) (or fit_cad(refine_segmentation=True)) refines a segmentation by variational shape approximation, splitting by seeded RANSAC, merging and absorbing until the patches stop moving; freeform_general=True fits a patch the height-field B-spline cannot hold over a conformal disc parametrisation with parameter correction and knot insertion; and , , and tighten the edges and measure each face's two-sided distance to the skin. (unreleased) turns what no primitive or feature explains into B-spline faces even without a CAD kernel, cutting non-discs along their topology, climbing a refine/re-parametrise/split ladder and giving every face it still demotes a reason, with the share of the skin left as facets in the report; lays such a region out in quads and fits them as B-spline patches that share their edge curves exactly, G1 across them, in both backends (a topology-optimised hook keeps well under one per cent of its area as facets, from a quarter), and the native solid comes back as trimmed NURBS for Kratos CAD JSON too; (experimental) refines with truncated hierarchical B-splines; (unreleased) fits related primitives together in both backends, so a turned part's cylinders, cones, tori and dome share one axis exactly, its shoulders are exactly perpendicular to it and nearly equal radii are one radius. See and .

seg = meshioplusplus.segment(skin, curvature_tolerance=1.0, min_patch_cells=20)
fit, report = meshioplusplus.fit_primitives(seg, tolerance=2e-3, return_report=True)
fit.field_data["fit:patch_params"]   # per patch: point, axis, radii, half-angle
meshioplusplus.write("part.step", fit)                      # analytic faces, no CAD kernel
solid = meshioplusplus.fit_cad(fit, "part_cad.step")         # trimmed and checked (OpenCASCADE)
meshioplusplus.fit_cad(meshioplusplus.read("bracket.stl"), "bracket.step")  # raw skin, one face per surface

Isogeometric layer: evaluate, refine, sample NURBS models

The NURBS model a reader keeps beside its mesh (mesh.info["nurbs"] for MFEM NURBS meshes, mesh.info["brep"] for Kratos CAD JSON, STEP and IGES) can be acted on directly, on every surface (Python, C++, C, Fortran, Julia, R, WASM, both CLIs, MCP): iga_refine inserts knots and elevates degrees exactly, with the geometry, the NURBS fields and a trimmed face's trimming loops unchanged; iga_evaluate returns points, first derivatives, normals and the trimming test at parametric coordinates; iga_sample samples a model into Lagrange cells again. A refined MFEM model writes back as an MFEM NURBS mesh. See isogeometric layer.

mesh = meshioplusplus.read("pipe-nurbs.mesh", file_format="mfem")
fine = meshioplusplus.iga_refine(mesh, subdivisions=2, elevate=1)   # exact: same geometry
meshioplusplus.write("fine.mesh", fine, file_format="mfem")         # patches intact
ev = meshioplusplus.iga_evaluate(fine, [[0.3, 0.7, 0.4]], derivatives=1)

Interfaces and contact

meshioplusplus.region_adjacency(mesh, regions=None) extracts conforming facets shared by Cell regions (or by cell blocks when fewer than two Cell regions exist). find_interface matches two parts by shared node ids or geometric proximity, contact_pairs projects a Point region onto master facets, and split_interface duplicates point fans along a Side region with optional cohesive elements. The API and tolerances are documented in doc/region_adjacency.md.

facets, report = meshioplusplus.find_interface(mesh, "part_a", "part_b", return_report=True)
contacts = meshioplusplus.contact_pairs(slave_mesh, "slave_nodes", "master_boundary")
split = meshioplusplus.split_interface(mesh, report["side_a"], add_cohesive=True)

The Python CLI provides region-adjacency, find-interface, contact-pairs, and split-interface; the native CLI provides region-adjacency.

Merge / combine

meshioplusplus.merge combines two or more meshes into one: it concatenates points (offsetting connectivity so indices stay valid), merges cell blocks by type, concatenates data (per a configurable data_policy), and tags each cell's origin. With weld=True it fuses coincident nodes across inputs within atol using a spatial hash (never O(N²)) — the standard way to stitch adjacent blocks into a watertight mesh. Overlapping set / field-data names are namespaced by source id. See doc/merge.md.

combined = meshioplusplus.merge([a, b, c])                 # concatenate
welded = meshioplusplus.merge([left, right], weld=True, atol=1e-8)  # fuse the shared interface

Editing (transform / clean / crop / split) and statistics

A bundle of dependency-free mesh-editing utilities:

• meshioplusplus.transform — apply an affine transform (translate / scale / rotate / 4×4 matrix / unit-scale) to the points; connectivity and data are carried through. See doc/transform.md.
• meshioplusplus.mirror — reflect across a plane with every cell right side out (a reflecting transform leaves them inverted); keep_original=True makes a whole model from a half, the points on the plane shared. See doc/mirror.md.
• meshioplusplus.clean — weld coincident points (spatial hash), drop degenerate and duplicate cells, and remove orphaned points, in one toggleable pass. See doc/clean.md.
• meshioplusplus.crop — extract the part of a mesh inside a bounding box or half-space, pruning unused points (mode="all"/"any"). See doc/crop.md.
• meshioplusplus.split — partition a mesh into several by cell type, connected component (flood-fill), region (cell_sets / integer tag), or one piece per named Cell region (by="regions", cross-binding, not a partition — overlapping regions overlap). See doc/split.md. meshioplusplus regions (CLI) / read_metadata(...)["regions"] lists a mesh's regions cheaply, without a full read where possible.
• meshioplusplus.compute_stats — geometric statistics (bounding box, centroid, per-type counts, area, signed/unsigned volume, inverted cells) — the geometric complement to info. See doc/stats.md.

out = meshioplusplus.transform(mesh, rotate=("z", 90))
out = meshioplusplus.clean(mesh, weld=True, atol=1e-8)
sub = meshioplusplus.crop(mesh, bbox=[0, 0, 0, 1, 1, 1])
pieces = meshioplusplus.split(mesh, by="type")             # {"triangle": ..., ...}
s = meshioplusplus.compute_stats(mesh)                     # dict of measures

Cell conversion (linearize / simplexify / elevate)

meshioplusplus.convert_cells converts a mesh's element representation — which cell types it is built from — while leaving the object it describes intact. See doc/convert_cells.md.

• mode="linearize" — every higher-order cell becomes its linear base (tetra10 → tetra, hexahedron27 → hexahedron), keeping the corner connectivity verbatim and pruning the nodes that become unreferenced.
• mode="simplexify" — every cell is decomposed into simplices of the same topological dimension (quad → 2 triangle, hexahedron → 6 tetra, wedge → 3, pyramid → 2, an n-gon into an (n−2)-triangle fan). No points are added, each parent's cell_data is replicated to its children, and every emitted simplex is positively oriented with volume conserved.
• mode="elevate" — every linear cell is promoted to its serendipity quadratic counterpart (triangle → triangle6, hexahedron → hexahedron20), adding one node per unique edge at the edge midpoint with point_data set to the endpoint mean.

linear = meshioplusplus.convert_cells(mesh, mode="linearize")
tets = meshioplusplus.convert_cells(mesh, mode="simplexify")   # hexes -> tetra
quadratic = meshioplusplus.convert_cells(mesh, mode="elevate")

Each mode is idempotent on cells it does not apply to, so it is safe on a mixed-order mesh, and output is byte-identical across mesh backends and thread counts.

Curved tessellation

meshioplusplus.tessellate isoparametrically subdivides the five higher-order cell types (triangle6/quad8/quad9/tetra10/hexahedron27, deliberately excluding hexahedron20 — see doc/tessellation.md) onto a refinement lattice, using each cell's own shape functions so a curved boundary tessellates onto the curve rather than being chopped straight. Unlike convert_cells(mode="simplexify"), it records a full provenance map (Tessellation.source_point/source_cell) from every synthetic point and tessellated cell back to the higher-order cell it was carved out of, so a prediction made on a tetrahedron is written onto the hexahedron it was carved out of via Tessellation.gather/.scatter/.aggregate. fields=True (default) interpolates the input's own point_data/cell_data onto the tessellated mesh. The kernel is in the C++ core, so the native CLI, every flat binding (C, Fortran, Julia, R, WASM) and the pipeline's Tessellate step run it too. See doc/tessellation.md.

tess = meshioplusplus.tessellate(mesh, levels=2)
mp.write("curved_bracket_tessellated.vtu", tess.mesh)

Polyhedral refinement (subdivide)

meshioplusplus.subdivide splits every eligible 3D cell into one polyhedral child per face, connected to a new interior point. refine and decimate both raise by name on a polyhedron, pointing at convert_cells(mode="simplexify") — both are built on fixed same-type subdivision templates, and an arbitrary polyhedron has none. subdivide needs no per-type table at all: tabulated types (reduced to corners for a quadratic variant) and existing polyhedron blocks are handled uniformly, so the same code covers every 3D cell type the mesh already supports. Automatically conforming (no closure, no hanging nodes), unlike refine. See doc/subdivide.md.

out = meshioplusplus.subdivide(mesh, record_parent_ids=True)

Polyhedral coarsening (agglomerate)

meshioplusplus.agglomerate merges groups of cells into single larger polyhedral cells — the many-to-one counterpart to subdivide. decimate raises by name on a polyhedron, pointing at convert_cells(mode="simplexify") — its fixed-template QEM edge collapse has no analogue for merging arbitrary polyhedral cells. agglomerate is a genuinely different algorithm: greedy seed-and-grow over the mesh's shared-face dual, absorbing face-adjacent neighbours by accumulated shared-face area until a target group size; each group emits one polyhedron whose faces are exactly its external boundary, conserving volume exactly. Non-volume blocks pass through unchanged; points are never pruned or renumbered (clean(mesh, remove_orphans=True) is the follow-up for a minimal point set). See doc/agglomerate.md.

coarse = meshioplusplus.agglomerate(mesh, target_group_size=8)
# fuse coplanar faces and keep groups compact
coarse, report = meshioplusplus.agglomerate(
    mesh, merge_coplanar_faces=True, min_sphericity=0.7, return_report=True
)

Refinement

meshioplusplus.refine subdivides cells into congruent children of the same cell type, increasing a mesh's resolution: line → 2, triangle → 4, quad → 4, tetra → 8, wedge → 8, hexahedron → 8, with levels=n applying the templates n times. Given a selection — a cell list, a region name, or a cell_data threshold — it refines only those cells and resolves the resulting hanging nodes, so the output is still conforming. See doc/refine.md.

New nodes sit at the midpoints of the parent's edges, quad faces and (hexahedron only) body, and carry the mean of that entity's corner values for every point_data array — so a linear field is interpolated exactly. Mid-edge and quad-face-centre nodes are shared between every cell touching the entity, so the refined mesh has no hanging nodes; each parent's cell_data row is replicated to its children.

record_hierarchy=True attaches two persistent cell_data arrays, refine:cell_id/refine:parent_id: a stable identity that survives across separate refine calls — a link between two meshes, not a tree inside one, resolved by a multigrid caller against the sequence of meshes it keeps.

fine = meshioplusplus.refine(mesh)                      # one level
finer = meshioplusplus.refine(mesh, levels=2)           # 64x the cells in 3D
tagged = meshioplusplus.refine(mesh, record_parent_ids=True)

# adaptive: refine the worst cells and close up conformingly around them
graded = meshioplusplus.refine(
    meshioplusplus.attach_quality(mesh),
    where="quality:scaled_jacobian < 0.3",
    record_levels=True,                                 # colour by refine:level
)

# multigrid: keep the coarse mesh, resolve the fine mesh's parent ids against it
fine = meshioplusplus.refine(coarse, cells=[4, 8], record_hierarchy=True)

Children inherit the parent's orientation (zero newly-inverted cells for a well-oriented input), and volume is conserved — exactly for tetra always, and for wedge/hexahedron when the parent is affine. Higher-order cells, pyramid, and ragged blocks have no same-type subdivision and raise by name.

Green-element undo

meshioplusplus.undo_green(coarse, fine) restores refine's transitional ("green") cells back to their original parent — the standard rule for selective refinement: before a new pass touches a region a prior pass already closed up, restore the transitional cell to its parent and re-split from scratch, rather than refining the transitional children directly (which degrades element quality without bound over repeated passes). It is a two-mesh operation, like interpolate: coarse is the mesh a prior refine(coarse, ..., record_hierarchy=True, record_levels=True) call was run on, fine is that call's output. Since refine's point map is always the identity, a green parent's exact connectivity and cell_data are already sitting, byte-for-byte, in coarse — so this is a lookup and substitution, not a reconstruction. See doc/undo_green.md.

fine = meshioplusplus.refine(coarse, cells=[4, 8], record_hierarchy=True, record_levels=True)
# ... later, decide to refine a different region ...
undone = meshioplusplus.undo_green(coarse, fine)
redone = meshioplusplus.refine(undone, cells=[12, 19], record_hierarchy=True, record_levels=True)

Decimation

meshioplusplus.decimate is refine's inverse: it reduces a surface mesh's face count by greedy quadric-error-metric (Garland–Heckbert) edge collapse, preserving shape, boundaries and features. Exactly one stopping criterion is given — ratio (fraction of faces to keep), target_faces, or max_error — and the output is all-triangle (quad/polygon blocks are triangulated first, block structure kept 1:1). See doc/decimate.md.

coarse = meshioplusplus.decimate(mesh, ratio=0.25)            # keep 25% of the faces
coarse = meshioplusplus.decimate(mesh, target_faces=5000)     # absolute face budget
coarse, report = meshioplusplus.decimate(mesh, max_error=1e-6, return_report=True)

Boundary vertices (once-used-edge test) and feature vertices (the endpoints of feature edges: an edge whose two faces' normals differ by more than feature_angle, default 30°) are pinned by default, so an open patch keeps its outline exactly and a cube keeps its corners; the link condition and a normal-flip guard reject any collapse that would change topology, create a non-manifold edge, or fold the surface. Float point_data blends along the collapsed edge; integer arrays keep the survivor's value. Volume meshes raise by name — run extract_surface first, then decimate the skin.

Volume decimation

meshioplusplus.decimate_volume is decimate's volume-mesh sibling — a separate operation, not a mode on it — reducing a tetrahedral mesh's cell count by greedy quadric-error-metric tet-edge collapse. Unlike decimate, boundary vertices participate by default (preserve_boundary=False): every vertex accumulates a quadric from its incident boundary-triangle planes only, so a purely interior vertex's quadric is exactly zero and interior-only edges are scored by squared length instead, always ranking behind boundary-touching collapses. See doc/decimate_volume.md.

coarse = meshioplusplus.decimate_volume(mesh, ratio=0.25)             # keep 25% of the tets
coarse = meshioplusplus.decimate_volume(mesh, target_cells=5000)      # absolute tet budget
coarse, report = meshioplusplus.decimate_volume(mesh, max_error=1e-6, return_report=True)

Validity is guarded by an exact vertex-link set-equality condition, a duplicate-tet check, and a tet-inversion guard, plus — for boundary-touching collapses — decimate's own ring/shared-face link condition and normal-flip check reused over the mesh's own outer skin. Tet-only: any non-tetra 3D block raises by name pointing at convert_cells(mode="simplexify").

Partitioning

meshioplusplus.partition decomposes a mesh into exactly N balanced pieces for domain decomposition — the count-driven complement to the criterion-driven split. See doc/partition.md.

• SFC (the default fallback, always available, dependency-free): cells are cut into contiguous ranges along a Hilbert space-filling curve of their centroids — equal-weight part sizes differ by at most one cell, weights=<cell_data> balances a per-cell cost instead, and the assignment is deterministic and byte-identical across mesh backends and thread counts.
• KaHIP (the optional quality path): the shared-face dual graph goes through KaHIP's serial kaffpa(), which actively minimizes the edge cut. Configure imbalance (default 3%), mode (fast/eco/strong, default eco — eco/strong carry the quality) and seed. KaHIP is MIT-licensed like meshio++ itself, so enabling it changes nothing about licensing; it is bring-your-own (-DMESHIOPLUSPLUS_WITH_KAHIP=ON + KAHIP_ROOT, Conan with_kahip, vcpkg feature kahip — next to the HDF5/zstd-style optional deps), links only the serial interface (no MPI), and pip install meshioplusplus[kahip] gives pure-Python installs the same quality path via the MIT kahip wheel. Requesting it where absent fails by name — never a silent downgrade.

pieces = meshioplusplus.partition(mesh, 4)                    # list of 4 meshes
labels = meshioplusplus.partition_labels(mesh, 4)             # per-block Int64 part ids
quality = meshioplusplus.partition(mesh, 16, method="kahip", mode="strong")

Pieces keep the input's block structure 1:1, so they recombine into the input: every cell lands in exactly one piece. ghost_layers=N instead grows each piece by N shared-node layers of its neighbours' cells (an MPI-style halo), tagged partition:ghost.

Surface remeshing (ACVD clustering)

meshioplusplus.remesh replaces a surface mesh's own triangulation with a new, near-uniformly-sized, well-shaped one at a caller-chosen vertex count, by approximated centroidal Voronoi diagram (ACVD) clustering — the one resolution-changing operation that does not inherit the input's element shapes, so it can raise quality at every target count rather than only remove elements. See doc/remesh.md.

out = meshioplusplus.remesh(mesh, num_clusters=5000)
out = meshioplusplus.remesh(mesh, num_clusters=5000, metric="quadric")       # preserves sharp edges/corners
out = meshioplusplus.remesh(mesh, num_clusters=5000, metric="anisotropic", max_anisotropy=4.0)

metric="isotropic" (default) is fast and rounds sharp features; "quadric" (Garland-Heckbert error) pins clusters onto edges/corners instead; "anisotropic" shapes clusters with a local curvature tensor, elongating elements along low-curvature directions. The output has entirely new points and connectivity — no point/cell map — so point_data/cell_data/named regions are dropped and field_data carries through; transfer a field onto the result with interpolate/conservative_interpolate.

Adaptation with MMG (optional backend)

meshioplusplus.adapt remeshes a simplicial mesh to a constant size, a size field or an anisotropic metric through MMG: mmg3d for tetrahedra, mmg2d for planar triangles, mmgs for surfaces. Cell and side regions come back under their names, required points and regions survive exactly, and metric_from_hessian turns a hessian into the metric. discretize_levelset embeds a level set's isovalue conformally between an inside and an outside region, and mesh_domain triangulates the planar domain a line mesh bounds, holes left out by a point inside each, and move moves a mesh by a displacement of part of its boundary with MMG's Lagrangian mode (an MMG built with ELAS, GPL-3.0). MMG is LGPL-3.0 and bring-your-own (-DMESHIOPLUSPLUS_WITH_MMG=ON, build/build-mmg.sh), so the published wheels do not carry it. See doc/mmg.md.

Volumetric remeshing (isosurface stuffing)

meshioplusplus.remesh_volume is remesh's volumetric sibling — the tet-mesh counterpart decimate_volume needed but could not itself provide, since QEM edge collapse can only remove elements. It retetrahedralizes a volume mesh (or a closed surface, unlike remesh) at a caller-chosen resolution by isosurface stuffing over a body-centered cubic (BCC) lattice, generating an entirely new mesh rather than working on the input's own cells. See doc/remesh_volume.md.

out = meshioplusplus.remesh_volume(mesh, cell_size=0.5)
out = meshioplusplus.remesh_volume(mesh, resolution=(64, 64, 64))
out, report = meshioplusplus.remesh_volume(mesh, cell_size=0.5, return_report=True)
print(report["num_tets"], report["num_non_manifold_edges"])

Every uncut lattice tet has a dihedral angle from a small, mesh-size-independent fixed set. warp_fraction (default 0.35) moves lattice vertices near the surface onto it, trading a small, measured chance of non-manifold boundary edges (reported in num_non_manifold_edges) for substantially better boundary tet quality; 0 disables warping and gives an exactly watertight but lower-quality boundary. Implemented from the published description of Labelle & Shewchuk's isosurface stuffing (SIGGRAPH 2007) only — no predicate library needed, unlike literal Delaunay/CVD tetrahedralization.

Volume meshes on an octree (graded tetrahedra fitted to a skin)

meshioplusplus.octree_mesh meshes the inside of any skin (STL, OBJ, a volume mesh's boundary; holes and self-intersections are tolerated through the fast winding number) on a 2:1-balanced octree whose leaves follow a size field (curvature, sharp features, thin gaps, a user array). The default output is a conforming graded tetrahedral mesh: a graded BCC lattice cut to the skin by isosurface stuffing, sharp corners and edges snapped, then optimised under a guard that never worsens the extreme dihedral angles; skin_grading="uniform" gives Labelle and Shewchuk's 10.7°–164.8° bound. output="castellated" gives the inside leaves as hexahedra with the surrogate_boundary side region and sdf:distance/sdf:closest_point, for embedded and Shifted Boundary Method solvers. output="hex" gives hexahedra only, on one grid at the skin's size, or with hex_grading="graded" the all-hexahedral dual of the graded octree, conforming through the transition schemes of Livesu, Pitzalis and Cherchi (ported from CinoLib). The same operation is in C++, C, Fortran, Julia, R, WASM, both CLIs (octree-mesh), the pipeline (OctreeMesh) and MCP. See doc/octree_mesh.md.

tets, report = meshioplusplus.octree_mesh(skin, min_size=0.02, max_size=0.2, return_report=True)
print(report["num_tetrahedra"], report["min_dihedral"], report["warnings"])
cells = meshioplusplus.octree_mesh(skin, output="castellated", min_size=0.02)

Surface meshes from CAD (conforming triangulation of a B-rep)

meshioplusplus.mesh_cad_surface triangulates the faces of a CAD model (a Kratos CAD JSON read, or STEP and IGES through OpenCASCADE) into one conforming surface mesh at a target edge size, with no CAD kernel: every B-rep edge is sampled once and shared by the faces that use it, so a closed solid comes out watertight, and each face is a constrained Delaunay triangulation of its parameters on Shewchuk's exact predicates, lifted onto the surface. A face that cannot be meshed is reported by index and reason, never meshed wrong, and the output is identical across backends and thread counts. min_angle refines each face by Delaunay refinement until its triangles' smallest angles on the surface reach the bound (guaranteed to 20.7 degrees away from small input angles and poles, which are reported), and smooth_iterations smooths the inside under a guard that never lowers a star's smallest angle. order=2 or 3 writes curved triangle6 or triangle10 elements whose nodes lie on the exact surfaces and edges, each proved valid everywhere inside by a Bézier bound on its Jacobian or kept straight-sided and reported. curvature_deviation and cells_across_gaps size the mesh by the model's curvature and its thin parts through one gradation-limited field (grading), pole_rings grades the triangles around a pole, and a triangle that would lift folded is refined away. Sizes can be given per face (face_sizes), per edge (edge_sizes) or by a background mesh, with the report naming the edges that hold a face finer than asked; virtual_edges skips the slivers of edges a model splits its rims into, and metric="varying" meshes each face in its surface's own metric, so stretched faces and the sphere's poles meet min_angle too. A face whose parameters are unusable (a side collapsed onto a point, a pole) can be meshed in a flattening of its own coarse mesh (reparam_condition_number), a sliver face with the neighbour it meets tangentially (virtual_faces), and the gaps a STEP or IGES model declares edge by edge are accepted with model_tolerances; lop_corners lops acute corners (Pav and Walkington 2005) and states, on planar faces, the angle bound that proves; metric_background sizes the triangles by an anisotropic metric from a field's Hessian (R |Λ| Rᵀ / (8 metric_error), Loseille and Alauzet 2011), so they stretch along a layer. untangle moves the nodes of folded curved elements on the surface and along the CAD edges to keep them curved and valid (Toulorge et al. 2013), and quads="dominant" recombines triangle pairs into quads by a maximum-weight matching (Remacle et al. 2012; pairs best first) while leaves no triangle at all, curved at and at , sized by the same sources; curved meshes work on flattened and merged faces too, and takes a field itself, its Hessian computed on the background. exposes the planar triangulator, with the same . A settings pipeline carries the model a CAD read returns to a step, Gmsh 4.1 output gives every CAD face its own physical surface, and fills a closed result with tetrahedra whose dihedral angles keep inside the isosurface-stuffing bound, or in one call (its octree backend, or Netgen keeping the triangles). The same operation is in C++, C, Fortran, Julia, R, WASM, both CLIs () and MCP, and a scheduled job measures it on real STEP models beside Gmsh's own meshes (). See .

surface, report = meshioplusplus.mesh_cad_surface(meshioplusplus.read("part.cad.json"), 0.1, return_report=True)
print(report["num_failed"], report["num_open_edges"], report["min_angle_deg"])

meshioplusplus.project_to_cad moves a mesh's points back onto the CAD model it was made from, by the cad:face and cad:edge cell data mesh_cad_surface writes, so a smoothing or an adaptation that knows no CAD can be put back on the model; points that land outside their face's trims are reported, never clamped. See doc/project_to_cad.md.

meshioplusplus.boundary_layer grows prismatic layers from a triangle and quad surface: num_layers layers of wedges and hexahedra along the most normal normals (Aubry and Löhner 2009), the first first_height thick and growing by growth, stopped short of the surface across a gap and where the last layer's front would cross itself or the surface (front_check), and halved, then collapsed into pyramids and tetrahedra (or, with collapse="layers", grown fewer layers) where a cell would fold. See doc/boundary_layer.md.

meshioplusplus.extrude and meshioplusplus.revolve sweep a surface into a volume, or a curve into a surface, along a vector (graded layers) or around an axis (the cells on the axis collapsing into pyramids, tetrahedra and wedges): the swept cells only, with bottom/top (start/end) and side regions, the input's regions and data carried; extruding rectangle gives box exactly. See doc/extrude.md.

meshioplusplus.tetrahedralize fills a closed triangle surface (a mesh_cad_surface result, a boundary layer's cap) with tetrahedra, nested shells becoming regions and cad:face side regions: through its native octree backend in every build (isosurface stuffing, the boundary fitted within a reported distance), or through Netgen's nglib, an optional backend that keeps the surface's triangles exactly (LGPL-2.1, linked dynamically when built with -DMESHIOPLUSPLUS_WITH_NETGEN=ON; backend="auto" takes it when present). See doc/tetrahedralize.md.

meshioplusplus.error_size turns estimate_error's per-cell indicator into a per-point size by Zienkiewicz and Zhu's adaptive procedure, a background for mesh_cad_surface, so an adaptation re-meshes from the CAD model. See doc/error_size.md.

layers, report = meshioplusplus.boundary_layer(surface, 1e-3, num_layers=8, growth=1.25, return_report=True)
slab = meshioplusplus.extrude(meshioplusplus.rectangle(dims=(8, 4)), (0.0, 0.0, 0.5), num_layers=6, ratio=1.2)
body = meshioplusplus.revolve(meshioplusplus.rectangle(dims=(4, 8)), angle=270.0, num_layers=12)
back = meshioplusplus.project_to_cad(meshioplusplus.smooth(surface, iterations=10), meshioplusplus.read("part.cad.json"))
volume, report = meshioplusplus.tetrahedralize(surface, max_size=0.1, return_report=True)  # Netgen if built, else octree

ODT remeshing (relocate + flip connectivity)

meshioplusplus.optimize_volume raises a tetrahedral mesh's worst element quality by ODT remeshing — relocating vertices AND changing connectivity. It is the genuine "ODT remeshing" and the third member of a trio whose other two each do half the job: smooth(method="odt") moves points on fixed connectivity, remesh_volume discards the input's tets for a fresh lattice mesh. It alternates the ODT vertex relocation with quality-improving topological flips (2-3 and 3-2), predicate-free — a flip is applied only when a pure signed-volume test finds the local configuration convex and the minimum scaled Jacobian strictly improves (no in-sphere/Delaunay predicate). See doc/optimize_volume.md.

out = meshioplusplus.optimize_volume(mesh)                       # a tetrahedral mesh
out, report = meshioplusplus.optimize_volume(mesh, return_report=True)
print(report["num_flips"], report["min_quality_before"], report["min_quality_after"])

The flips touch only interior faces/edges, so the boundary surface is preserved exactly (watertight in ⇒ watertight out); the point set is invariant, so point_data/field_data and named Point regions carry through while cell_data/Cell/Side regions are dropped. Tet-only, and (like remesh_volume) C++-core only.

Smoothing

meshioplusplus.smooth relaxes point coordinates to improve element shape, leaving topology and every data value alone: only the points move. See doc/smooth.md.

Three methods. Laplacian (x <- x + lambda*L(x), the edge-neighbour centroid displacement) smooths strongly per pass but shrinks — over 40 iterations on a jittered 8×8 quad grid it contracts the bounding box by 57%. Taubin (the default) follows each +lambda pass with a larger-magnitude -mu pass that deliberately un-shrinks, leaving the same grid 3.6% smaller. ODT (optimal-Delaunay-triangulation smoothing, tet-only) instead moves each free interior vertex to the closed-form volume-weighted average of its incident tets' circumcenters — the "ODT" half of the volumetric-remeshing roadmap item, closed as smoothing on existing connectivity rather than remeshing. Neighbours for Laplacian/Taubin are the nodes joined by an actual cell edge, not the element clique, so a structured hex block is a fixed point rather than being bevelled toward a sphere. Boundary nodes, feature nodes (the endpoints of feature edges — on a surface mesh its own creases as well as the corners of its boundary), an optional frozen mask, and the nodes of blocks whose edge topology is unknown are all pinned by default, and the inversion guard rejects any move that would turn a valid cell inverted.

relaxed = meshioplusplus.smooth(mesh)                          # 10 Taubin iterations
harder = meshioplusplus.smooth(mesh, iterations=40)            # shrink-free even so
lap = meshioplusplus.smooth(mesh, method="laplacian", lambda_=0.4)  # note the underscore
odt = meshioplusplus.smooth(tet_mesh, method="odt", iterations=10)  # tet-only, C++-core only
out, report = meshioplusplus.smooth(mesh, return_report=True)  # nodes moved, max displacement

lambda_ carries a trailing underscore because lambda is a Python keyword, and a negative value means "this method's own default" (0.5 Laplacian, 0.33 Taubin). Point and cell counts, connectivity, cell_data, field_data, point_data values and the points array's dtype all come through unchanged, and output is byte-identical across mesh backends and thread counts.

Field transfer (interpolation)

meshioplusplus.interpolate samples a source mesh's data arrays onto a target mesh — the first cross-mesh operation that transfers data (diff compares, merge concatenates). The result is a copy of the target, its own geometry/data/sets preserved exactly, with the requested source arrays sampled on: source point_data at the target's points, source cell_data by nearest source-cell centroid. method="nearest" (default) copies the nearest source point's value bit-for-bit; method="barycentric" simplexifies the source first and interpolates linearly — exact on a linear field — with default_value/extrapolate deciding what happens outside the source domain. Both search grids are bucket-grid spatial hashes (never O(N²)), and output is byte-identical across backends, thread counts and the C++/numpy boundary. See doc/interpolate.md.

coarse = meshioplusplus.read("solution.vtu")   # carries point_data "T"
fine = meshioplusplus.read("remeshed.vtu")
mapped = meshioplusplus.interpolate(coarse, fine, method="barycentric")
meshioplusplus.write("mapped.vtu", mapped)

meshioplusplus.conservative_interpolate is a separate, mass-preserving sibling: over the region the two meshes share, sum(target value * target measure) equals sum(source value * source measure), a property interpolate's pointwise sampling does not have. Both meshes are simplexified first (accepting ragged/polyhedron blocks for free), overlapping simplex pairs are measured with an exact geometric clip (Sutherland-Hodgman in 2D, a bounded tetrahedron-tetrahedron clip in 3D), and a target cell's value is the overlap-measure-weighted mean of every source cell it intersects. Unlike interpolate, an unset arrays transfers every source point_data and cell_data array. point_data is transferred by composition (point_data_to_cell_data → the same clip algorithm → cell_data_to_point_data), a layered approximation rather than exact nodal conservation. C++-core only, with no pure-Python fallback (the 3D clip's discrete branches could disagree near a degenerate overlap). See doc/conservative_interpolate.md.

Slicing / cross-sections

meshioplusplus.slice computes the planar cross-section of a mesh — the actual intersection of the mesh with a plane, one topological dimension below the cut cells: a 3D volume mesh yields a triangle/quad surface, a 2D surface mesh a line mesh. Unlike crop (plane mode), which keeps whole cells on one side, slice computes the intersection and lowers the dimension. It uses robust marching tetrahedra (the input is simplexified first, so each cell's cross-section is a well-defined convex primitive), deduping crossing points on shared edges into single nodes so the section is watertight, and winding every face consistently toward the +normal side. Each section cell inherits its parent's cell_data; record_parent_ids=True attaches slice:parent_cell, and point_data is interpolated at the cut. Output is byte-identical across backends, thread counts and the C++/numpy boundary. See doc/slice.md.

vol = meshioplusplus.read("part.vtu")                          # a tetra/hex/wedge mesh
section = meshioplusplus.slice(vol, origin=(0, 0, 0.5), normal=(0, 0, 1))
meshioplusplus.write("section.vtu", section)                  # a triangle/quad surface at z=0.5

(slice shadows the Python built-in only as a module attribute — meshioplusplus.slice is intended.)

Isosurfaces / contours

meshioplusplus.isosurface computes the level set of a scalar field — the locus where a point_data array equals a given isovalue, as a mesh one topological dimension below the cut cells: a 3D volume mesh yields a triangle/quad surface, a 2D surface mesh a line contour. It is the data-driven sibling of slice (which cuts where the distance to a plane is zero) and shares its marching-tetrahedra cutter, so contours are watertight in the same way. The field must be point_data — cell_data is piecewise constant and has no level set, so naming one raises and points at cell_data_to_point_data. Several isovalues land in one mesh, cut in ascending order and tagged per cell with a Float64 iso:value and an Int64 iso:index (the ordinal — the integer tag split(by="region", tag=…) needs). The contoured field reads back as exactly the isovalue on the cut points, faces are wound toward increasing field, and an out-of-range isovalue is an empty contour rather than an error. Output is byte-identical across backends, thread counts and the C++/numpy boundary. See doc/isosurface.md.

vol = meshioplusplus.read("part.vtu")                          # carrying point_data["T"]
shells = meshioplusplus.isosurface(vol, "T", [300.0, 350.0, 400.0])
meshioplusplus.write("shells.vtu", shells)                     # three tagged contour surfaces

Field derivatives (gradient / divergence / curl)

meshioplusplus.gradient differentiates a point_data field: its gradient, divergence or curl. meshio++ could already transform, transfer, summarize and contour a field — this is what lets it differentiate one, which is the missing input for contouring a derived quantity (|∇T|, vorticity) and for the gradient-based error indicators that drive the selective refine. Two methods: Green-Gauss (the default) applies the divergence theorem over the cell, fanning each face into triangles about its corner average, which is exact for a linear field on any cell — planar faces or not, because the fan surface is closed; least-squares fits over the cells sharing a node and falls back to Green-Gauss (counted, never silently wrong) on a degenerate neighbourhood. An nc-component input yields 3·nc gradient components laid out [component][derivative], so a scalar gives (n, 3) and a 3-vector (n, 9); divergence gives 1 and curl 3. Output is Float64, named <input>:gradient / :divergence / :curl, at either the cell (default) or point location. A cell_data input raises by name — a piecewise-constant field has no derivative — and cells that cannot be differentiated yield NaN and are counted rather than approximated. Geometry, regions and existing data pass through bit-identically; output is byte-identical across backends, thread counts and the C++/numpy boundary. See doc/gradient.md.

import numpy as np

vol = meshioplusplus.read("solution.vtu")                      # carrying point_data["T"]
g = meshioplusplus.gradient(vol, "T", location="point")        # point_data["T:gradient"], (n, 3)

grad = np.asarray(g.point_data["T:gradient"])
g.point_data["gradT"] = np.sqrt((grad**2).sum(axis=1))
shells = meshioplusplus.isosurface(g, "gradT", [2.0])          # contour where T changes fastest

These operations are exposed across every binding surface (Python, C API, Fortran, WASM) and as the CLI verbs meshioplusplus quality, meshioplusplus extract-surface, meshioplusplus reorder, meshioplusplus diff, meshioplusplus merge, meshioplusplus transform, meshioplusplus clean, meshioplusplus crop, meshioplusplus slice, meshioplusplus split, meshioplusplus stats, meshioplusplus convert-cells, meshioplusplus tessellate, meshioplusplus subdivide, meshioplusplus agglomerate, meshioplusplus refine, meshioplusplus undo-green, meshioplusplus partition, meshioplusplus remesh, meshioplusplus remesh-volume, meshioplusplus optimize-volume, meshioplusplus smooth, meshioplusplus interpolate, meshioplusplus conservative-interpolate, meshioplusplus isosurface, meshioplusplus curvature, meshioplusplus normals, meshioplusplus repair, meshioplusplus shrinkwrap and meshioplusplus sobolev-deform (plus meshioplusplus data gradient, meshioplusplus data hessian, meshioplusplus data estimate-error and meshioplusplus data integrate, mesh operations grouped under data because that is where a user looks for them).

Second derivatives (Hessian)

meshioplusplus.hessian computes the Hessian (second derivative) of a scalar point_data field — gradient's companion one order further, for curvature-based adaptive refinement. A composition of two gradient calls, not a new numerical kernel: the field is differentiated once (point location), then that (n, 3) gradient is differentiated again with the default gradient operator, producing (n, 9) — the flattened row-major 3x3 Hessian. method forwards to both internal passes. A field that is at most linear has an exactly zero Hessian everywhere — the one mesh-shape-independent guarantee; a genuinely quadratic field's composition is exact on a structured/symmetric mesh away from its own boundary and a good, standard, but genuinely approximate curvature estimate on an irregular mesh. Input must have exactly one component. A curvature-driven refinement indicator needs no new code: data_calc's norm(...) on the 9-component output is exactly its Frobenius norm, ready for refine's where selector. See doc/hessian.md.

h = meshioplusplus.hessian(vol, "T")                            # cell_data["T:hessian"], (n, 9)
curv = meshioplusplus.data_calc(h, "norm(`T:hessian`)", location="cell", output="curv")
adapted = meshioplusplus.refine(curv, where="curv > 3.0")       # refine where curvature is largest

Error estimation (Zienkiewicz-Zhu recovery + marking)

meshioplusplus.estimate_error estimates the per-cell recovered-gradient error of a point_data field, and can mark cells for refinement — the piece that closes the adaptive loop: gradient differentiates and selective refine's --where consumes any scalar cell_data predicate; this produces one, so estimate → mark → refine is three verbs and no new numerical code. It is a composition, not a new kernel: gradient (Green-Gauss, cell location) → the measure-weighted point↔cell averaging round trip recovers a smoothed gradient, and eta_K = sqrt(|measure| · sum((recovered − raw)²)) per cell is the standard ZZ indicator. marking turns it into a boolean error:marked array: "absolute" (threshold), "fraction" (top N by indicator), or "dorfler" (the smallest indicator-descending prefix covering a bulk fraction of the total — the usual AMR criterion). Cells that cannot be evaluated read NaN in the indicator and 0 (never NaN) in the marking array, counted and excluded from the global error. See doc/error.md.

out, report = meshioplusplus.estimate_error(
    vol, "T", marking="dorfler", marking_value=0.6, return_report=True,
)
print(report["global_error"], report["num_marked"])
adapted = meshioplusplus.refine(out, where="error:marked > 0.5")

Field integration (total / mean, per-region)

meshioplusplus.data_integrate computes a cell-measure-weighted total and mean of one or more cell_data arrays — gradient's integration counterpart (gradient differentiates a field, this integrates one), for a density field's total mass, a heat-flux field's total power, or an occupied volume. Every sum is weighted by |measure(cell)| (the cell's own length/area/volume); a cell whose measure is not computable, or a component whose value is non-finite, is excluded from that component's numerator and denominator — never given a fallback weight of 1, unlike cell_data_to_point_data's own measure weighting, since a silent unit-weight substitution would corrupt a physical total. Reported for the whole mesh and independently for every named Cell region — regions are not a partition, so a cell in two regions contributes fully to both. A point_data-only name raises by name, pointing at point_data_to_cell_data. The mesh is never modified; this is a read-only report, like data_info and compute_stats. See doc/field_integration.md.

report = meshioplusplus.data_integrate(mesh, arrays=["density"])
report[0]["domain"]["total_per_component"]   # sum(value * |measure|) -- total mass
report[0]["regions"]                         # the same, independently, per named Cell region

Geometry helpers (curvature, normals, repair, shrinkwrap, Sobolev deformation, grids and SDF)

Six operations built for the physics-ML surface, each with its own CLI verb, MCP tool, and Python-API and doc page: meshioplusplus.compute_curvature (per-vertex mean/Gaussian curvature by the angle defect and cotangent Laplace-Beltrami operator — doc/curvature.md), meshioplusplus.compute_normals (point/cell normals with an optional crease split — doc/normals.md), meshioplusplus.repair (fixes a surface's orientation, holes and pinched vertices, the defects clean does not touch — doc/repair.md), meshioplusplus.shrinkwrap (projects a mesh's points onto a target surface in one pass — doc/shrinkwrap.md), and meshioplusplus.sobolev_deform (moves points by a raw displacement field filtered through the mesh's own P1 operators — doc/sobolev_deform.md). meshioplusplus.voxelize/compute_sdf build a regular grid around or over a surface, filled by occupancy or signed distance — doc/grids.md and doc/sdf.md.

Parametric primitives (grids, disks, spheres, tori, cones, balls, shells, curves, point clouds, reference cells)

Thirty-four generators build a shape from parameters, so a test, a notebook or a demo needs no fixture file: meshioplusplus.rectangle, box, disk, uv_sphere, icosa_sphere, cylinder, the planar domains reference_triangle, ngon, annulus (or a sector of one, with a seam pair), lshape, plate_with_hole and honeycomb (a ragged polygon block of hexagons), reference_cells (one cell of each type in its reference coordinates), the curves and point clouds line, circle, helix, sphere_points (a Fibonacci lattice) and box_points (seeded, the same cloud on every platform), the surfaces tetra_sphere, octa_sphere, cube_sphere, geo_sphere, torus, moebius (non-orientable), tube, cone, capsule, box_surface and platonic, and the solids ball (the seven-block butterfly), spherical_shell and pipe (graded walls), fichera and solid_torus. Surfaces come out oriented outward, volumes positive, the boundary as named side regions (xmin … zmax, rim, side/top/bottom, inner/outer, hole), and order=2 places the mid-edge nodes on the exact geometry, so a quadratic sphere or cylinder is genuinely curved; complete=True gives the full-Lagrange quad9, hexahedron27 and wedge18, and order 3 to 10 the VTK_LAGRANGE_* family, every node on the geometry too. The same generators are on every surface — C, Fortran, Julia, R, WASM, both CLIs (meshioplusplus generate box block.vtu --dims 8,8,8 --cell-type tetra), MCP and the pipeline ({"Op": "Generate", ...}, which needs no Input). See parametric primitives.

block = meshioplusplus.box(dims=(4, 4, 4), cell_type="tetra", order=2)   # tetra10, six side regions
ball = meshioplusplus.uv_sphere(radius=2.0, nu=24, nv=12, order=2)      # triangle6 on the exact sphere

Data operations (rename / average / calc / condition / invariants / summarize)

A second bundle operates on the data arrays a mesh carries (point_data / cell_data / field_data) rather than on its geometry, which none of them ever modifies:

• meshioplusplus.data_rename / data_drop / data_keep — rewrite which arrays a mesh carries and under what names; values, dtypes and shapes are copied verbatim. See doc/data_manage.md.
• meshioplusplus.point_data_to_cell_data / cell_data_to_point_data — move data between locations by averaging, optionally weighted by cell area/volume. See doc/data_average.md.
• meshioplusplus.data_calc — derive a new array from an elementwise expression (+ - * /, parentheses, abs/sqrt/min/max/norm) evaluated by a hand-written parser — no external parser library, no arbitrary-code path. See doc/data_calc.md.
• meshioplusplus.data_condition — clamp, normalize to a target range, or standardize to zero mean / unit standard deviation, per component or by row magnitude. See doc/data_condition.md.
• meshioplusplus.tensor_invariants — von Mises, principal, hydrostatic and deviatoric fields of a symmetric or general 3x3 tensor array. See doc/tensor_invariants.md.
• meshioplusplus.data_info — a read-only per-array summary (dtype, shape, components, min/max/mean, NaN/inf counts) — the data-side complement to info and compute_stats. See doc/data_info.md.

out = meshioplusplus.data_calc(mesh, "norm(velocity)", location="point", output="speed")
out = meshioplusplus.point_data_to_cell_data(out, keys=["speed"], suffix="_c")
out = meshioplusplus.data_condition(out, "cell", ["speed_c"], mode="normalize")
out = meshioplusplus.data_rename(out, "point", "T", "temperature")
arrays = meshioplusplus.data_info(out)                     # list of per-array dicts

These are likewise exposed across every binding surface, and as the ten CLI verbs under the meshioplusplus data group (info, rename, drop, keep, to-cell, to-point, calc, clamp, normalize, invariants). See doc/data_operations.md. A second nested group, meshioplusplus dataset (add, list, split, tag, annotate), curates the hand-editable dataset manifests used for ML training collections (Python CLI only, like data export).

Time series

The XDMF format supports time series with a shared mesh. You can write times series data using meshio++ with

with meshioplusplus.xdmf.TimeSeriesWriter(filename) as writer:
    writer.write_points_cells(points, cells)
    for t in [0.0, 0.1, 0.21]:
        writer.write_data(t, point_data={"phi": data})

and read it with

with meshioplusplus.xdmf.TimeSeriesReader(filename) as reader:
    points, cells = reader.read_points_cells()
    for k in range(reader.num_steps):
        t, point_data, cell_data = reader.read_data(k)

Transient / multi-file datasets

Most formats cannot express time at all, so transient output usually arrives as a set of files (out_0000.vtu … out_0500.vtu). meshio++ treats such a set as one ordered dataset:

meshioplusplus convert 'out_*.vtu' series.xdmf         # fan-in (quote the glob!)
meshioplusplus convert series.xdmf 'step_{step}.vtu'   # fan-out
meshioplusplus pipeline transient.json                 # chain applied per step

for time, mesh in meshioplusplus.read_sequence("out_*.vtu"):  # lazy
    ...
meshioplusplus.write_sequence("series.xdmf",
                              meshioplusplus.read_sequence("out_*.vtu"))

Ordering is natural-numeric, so out_10.vtu follows out_9.vtu; each step's time comes from an explicit list, the file, its filename or its index, and which one applied is reported. Fan-in and fan-out stream — one mesh is alive at a time, whatever the step count — and a multi-step input aimed at a single-step output is an error naming {step}, never a silent write of step 0. meshioplusplus.resample_sequence("out_*.vtu", "u_{index}.vtu", times="0:2:0.1") puts a run on new times (a linear blend, nearest or previous step, or another run's own times with times_from) holding at most two steps in memory, which aligns two solvers before a diff. Available from Python, both CLIs, C, Fortran, Julia and R. See doc/sequences.md.

Interactive viewer

The browser viewer — reading, rendering and converting entirely client-side.

One call, two backends:

import meshioplusplus

mesh = meshioplusplus.read("part.msh")
meshioplusplus.view(mesh)                      # pick a backend automatically
meshioplusplus.view(mesh, backend="polyscope") # a native desktop window
meshioplusplus.view(mesh, backend="browser")   # vtk.js, in a browser or notebook

The desktop backend is Polyscope, an optional Python-only extra (pip install meshioplusplus[viewer]). It draws solids you can slice into, colours by any point or cell array, and renders headless screenshots for CI and docs:

meshioplusplus.screenshot(mesh, "part.png", color_by="temperature")

The bundled example.msh bracket coloured by element quality — this image is generated by screenshot() itself.

The browser backend needs nothing extra. The same app is hosted as a live demo: drag in any supported format, colour by point or cell data, and convert and download to another format — all client-side, with no server and no upload. Since it runs the WebAssembly build, every format meshio++ reads works there too.

A second page, the dataset manager, curates the dataset manifests used for ML training collections visually — an overview of every manifest in a directory (cards with split balance, health badges and thumbnails, plus a manifest diff view, and — with an optional local companion process — server-side scans and PhysicsNeMo training launched, followed, compared and previewed from the page; see the dashboard page), directory picking with in-place manifest save (Chromium), per-entry previews with a time-series scrubber, and NaN/Inf/quality scanning — against the same hand-editable JSON the CLI and Python API use:

Splits, tags and groups on the left; a live preview and per-array summary on the right — all against the manifest JSON the CLI and Python API read and write.

From the command line:

meshioplusplus view part.msh
meshioplusplus screenshot part.msh part.png --size 1600 1200

See the viewer docs for how volume meshes are handled and what each backend can and cannot do.

Terminal rendering

Where there is no window, browser or GPU — an SSH session on a login node, a container, a CI log — snapshot draws the mesh anyway. A dependency-free software rasterizer (depth buffer, flat or smooth shading, element or feature edges, a mapped data field with a colour bar) turns it into pixels, and the pixels into Unicode block characters with terminal colours, a PNG, an HTML page or an asciinema orbit. The output is the same bytes on every platform, thread count and CLI.

meshioplusplus snapshot part.vtu -                                    # draw it in this terminal
meshioplusplus snapshot result.vtu - --color-by temperature --colorbar --encoding braille
meshioplusplus snapshot part.vtu part.png --edges feature --axes      # a PNG, no display needed
meshioplusplus tui part.vtu                                           # orbit and zoom it right here

print(meshioplusplus.render_text(mesh, color_by="temperature"))      # ANSI text for this terminal
image = meshioplusplus.render_image(mesh, 800, 600, supersample=2)   # an (H, W, 4) uint8 array
meshioplusplus.snapshot(mesh, "part.png", edges="all")               # .png .txt .ansi .html .cast
meshioplusplus.tui(mesh)                                             # the interactive viewer

Field rendering colours by an array, a data_calc expression or a tensor invariant on a linear, log or symlog scale with percentile clipping and a legend, draws categories, regions, isolines, vector arrows, streamlines that follow a vector field and a warp with its undeformed outline, and offers quality, orientation and free-edge diagnostic views. The same options reach C, Fortran, Julia, R, WebAssembly and the settings pipeline's Snapshot step.

tui is the interactive viewer for a machine with no display: drag to orbit, wheel to zoom, click to probe a cell, x/y/z to cut the model open, [ ] to step through a time series (--follow watches a run that is still writing), --compare to put two meshes side by side, a : command line that takes the flags of snapshot, and sessions saved as JSON; only the changed cells are sent, and the terminal is restored however it ends. view() picks it by itself when there is no display but a terminal. For fun, --theme synthwave puts a banded sunset and a grid floor behind the model in neon, and --music plays a generated loop on whichever audio player the system has (off unless asked for, and silent in CI and over SSH).

It is a preview, not a replacement for the viewers above; screenshot() falls back to it when Polyscope is not installed, and the MCP server's render_mesh tool lets an agent look at a mesh. See doc/tui.md.

Interoperability

Hand a mesh straight to the tools you reach for next — no file round-trip, and the numpy buffers are shared, not copied, wherever the target accepts them as they are:

import meshioplusplus

mesh = meshioplusplus.read("bracket.msh")

grid = meshioplusplus.to_pyvista(mesh)      # a pyvista.UnstructuredGrid
tm = meshioplusplus.to_trimesh(mesh)        # a trimesh.Trimesh (triangles only)

Both directions exist (from_pyvista, from_trimesh). Mixed-type meshes are the normal case, and named regions ride along as region:<name> mask arrays plus a metadata sidecar, so even a gmsh physical group's integer tag survives a PyVista round-trip. zero_copy_only=True turns any step that would copy into a named error instead of a silent one.

Data arrays also export to Apache Arrow and Parquet for the analytics stack, or straight to a pandas / polars frame with no Parquet detour:

df = meshioplusplus.to_pandas(mesh, location="cell")    # a pandas.DataFrame
pf = meshioplusplus.to_polars(mesh, location="point")   # a polars.DataFrame

meshioplusplus.write_parquet(mesh, "cells.parquet", location="cell")

import pandas
pandas.read_parquet("cells.parquet").head()

meshioplusplus data export bracket.msh cells.parquet --location cell

Multi-component arrays keep their shape wherever the target can hold it — Arrow fixed_size_list columns, polars pl.Array columns — while to_pandas (whose columns are one-dimensional) flattens them into suffixed columns v_0/v_1/v_2 with the grouping recorded in df.attrs. The mesh's counts, cell types and region names travel in the Arrow schema metadata / df.attrs. This is a data export, not a mesh format — it does not round-trip geometry and is deliberately not in the format registry.

PyVista, trimesh, pyarrow, pandas and polars are Python-only optional extras (pip install meshioplusplus[interop], or one at a time with [pyvista] / [trimesh] / [arrow] / [pandas] / [polars]). They are kept out of [all], which means "the optional dependencies the formats need". None of them reaches the C++/WebAssembly/C/Fortran core, which stays dependency-free.

Meshes also hand off to the GPU through the standard exchange protocols — DLPack as the primary export (which covers host arrays too), __cuda_array_interface__ consumed on the way back:

gpu = meshioplusplus.to_cupy(mesh)            # one host→device transfer per array
gpu.point_data["T"] *= 2.0                    # any CuPy / RAPIDS kernel
meshioplusplus.write("out.vtu", meshioplusplus.from_cupy(gpu))

The host→device move is always a bus transfer — what this removes is the file round-trip and every extra copy around it. to_dlpack(mesh) exports host arrays any DLPack consumer (PyTorch, JAX, Numba, …) adopts in place. There is deliberately no [gpu] extra: CuPy wheels are CUDA-version-specific, so install the one matching your toolkit (e.g. pip install cupy-cuda13x for CUDA 13.x, cupy-cuda12x for 12.x).

And for machine-learning pipelines, meshes become graphs, feature matrices, datasets and framework tensors directly:

ei = meshioplusplus.edge_index(mesh)                  # (2, E) int64 — PyG/DGL layout
fm = meshioplusplus.feature_matrix(mesh, "point")     # (N, F) float64 + recorded columns
meshioplusplus.write_dataset("out_*.vtu", "dataset/") # mesh_id-keyed partitioned Parquet
t = meshioplusplus.to_torch(mesh)                     # torch tensors, adopted zero-copy

feature_matrix's column order is a stated, versioned contract recorded in the returned schema, so training and inference cannot silently disagree; write_dataset streams a glob / directory / transient series into hive-partitioned Parquet (or chunked zarr/hdf5 groups, [zarr]/h5py) with a strict shared schema and a JSON manifest; to_torch/to_jax adopt the DLPack payload per framework (no [torch]/[jax] extra, deliberately — the CuPy precedent). See the ML docs.

A collection of solution outputs is catalogued by a hand-editable dataset manifest (DatasetManifest — sources, train/valid/test splits, tags, groups, notes; curated in Python, by the meshioplusplus dataset CLI group, over MCP, or with a text editor, all against the same JSON), and the PhysicsNeMo adapter (meshioplusplus.physicsnemo) trains straight off it: graph_sample builds the MeshGraphNet tensor set per mesh, field_stats/edge_stats stream normalization stats in PhysicsNeMo's own convention, make_dataset yields a PyTorch Geometric dataset and make_reader a Gen-2 Reader, and run_training/predict (also python -m meshioplusplus.physicsnemo.train --spec, and the dashboard's Start training) run the loop and write predictions back onto the mesh — for five model families (MeshGraphNet on the mesh graph; SRResNet, FNO and AFNO on regular grids, the operators 2-D through a thin-axis squeeze; and DeepONet for parameters-in/field-out over each entry's Metadata) — with worked, GPU-executed end-to-end examples in example/physicsnemo/.

See the interoperability docs for the full mapping tables, the zero-copy contract, and the Open3D/DOLFINx design sketch, and the GPU docs for the device handoff.

For a grid-shaped model — a 3-D CNN, or a superresolution network taking a coarse grid to a fine one — the shape a mesh has to become is a dense array, and sample_grid is that step:

spec  = mio.GridSpec.from_mesh(mesh, resolution=(64, 64, 64))
array = mio.sample_grid(mesh, spec, fields=["T", "vel"])   # (C, D, H, W), channels first then z,y,x
array.coverage                                             # how much of the grid is actually inside the mesh
result = mio.scatter_grid(prediction, mesh)                # and back onto the mesh, trilinearly

Grids are ordinary hexahedron meshes, so view, crop and every writer work on one; cache them as .vti, which stores the lattice exactly. power_spectrum reports whether a super-resolved field carries the right small-scale content, which a pointwise error cannot see. See mesh and regular grids.

For a point-cloud model — Transolver, FLARE, DoMINO, anything whose cost is quadratic in its token count — the shape is a fixed number of points, and select_points is that step:

budget = mio.select_points(mesh, 4096)                     # farthest-point sampling; method="grid" scales to millions
tokens = budget.take(mio.feature_matrix(mesh).matrix)      # (4096, F), the column contract intact
cloud  = mio.subsample_points(mesh, budget)                # a vertex-block point cloud every writer accepts

The selection is kept in selection order, so one farthest-point budget at 8192 serves 4096 and 2048 by slicing. See point-cloud budgets.

For a particle method there is no connectivity to build a graph from at all — a smoothed-particle or discrete-element state is a cloud of positions, and what makes two particles interact is the interaction radius:

edges = mio.proximity_graph(cloud, radius=0.015)                    # or method="knn", max_neighbors=16
edges = mio.proximity_graph(cloud, radius=0.015, box_size=[2, 2, 2])  # periodic: minimum image
attrs = mio.edge_vectors(cloud.points, edges)                       # displacement, then its norm

For a transient surrogate the shape is a window — the last K states of every node, oldest first — and iter_windows builds them in three schemes over a manifest, while rollout feeds a model its own predictions back and reports how the error grows, which a one-step validation loss cannot show. Augmentation randomizes each case's pose per epoch, drawn deterministically from (seed, epoch, index) and replayed coherently across a paired input and target.

A trained model answers any mesh it is given, and the answer for a part unlike anything it saw is finite, plausible and wrong. GeometryGuard is the check it cannot make for itself — fit a description of the shapes a dataset contains, score a new one against it, and get back the descriptors that put it there. The descriptors are deliberately not scale- or position-invariant, because a scaled part is a different part.

Once a model is trained, applying it to one mesh needs no manifest at all — meshioplusplus predict model.mdlus part.vtu part_pred.vtu, or predict_file from Python. Everything comes from the checkpoint's own model card, and a mesh carrying no truth predicts anyway rather than reporting an error against itself.

The same call adds world edges beside a mesh's own — contact between surfaces that are near but not connected — and bistride_hierarchy coarsens a graph so a signal crosses the mesh in logarithmically many message-passing steps. See proximity graphs.

If PhysicsNeMo itself is unfamiliar, PhysicsNeMo basics is a fourteen-page map of the framework — what a Module and a .mdlus checkpoint are, which of the 25 architecture families fits the shape of your data, how simulation output becomes batched tensors, and where meshio++ ends and the framework begins. Each page closes by naming what meshio++ supplies, or by saying plainly that nothing does.

MCP server

Every operation in this README is also exposed to AI agents as a tool over the Model Context Protocol — reading/writing all the formats, conversion, and the full mesh- and data-operation suite:

pip install "meshioplusplus[mcp]"     # the mcp SDK needs Python >= 3.10
claude mcp add meshioplusplus -- meshioplusplus-mcp

Then ask the agent to convert, inspect, slice, partition, … and it drives the 87 tools itself. Tools are file-path based (optionally sandboxed with --root DIR), parsed inputs are cached between calls and re-read when a file changes, and every report is strict JSON. meshioplusplus-mcp --http (pip install "meshioplusplus[dashboard]") serves the same tools over HTTP — MCP over streamable HTTP for agents, plus the JSON API the browser dataset dashboard uses as its local companion process. See the MCP docs for the tool table and client setup.

Blender add-on

Blender ships Python and reads almost no FEA formats. The add-on puts all 43 of meshio++'s behind File > Import.

It is a Blender 4.2+ extension, so the meshio++ wheel travels inside the zip — no pip step, no network at install time, nothing to configure. Download meshioplusplus-<version>-<platform>.zip from the latest release and drag it into a Blender window, or use Edit / Preferences / Add-ons / Install from Disk.

Volume meshes arrive as their boundary surface — Blender has no tetrahedron — with each cell_data array carried through to the faces of the cell that owned it, so a solid can still be coloured by its material. Quads and n-gons are kept, not triangulated. point_data, cell_data and named regions land as Blender attributes on the right domain.

The same bridge is two public functions, usable from Blender's scripting console: to_blender(mesh) and from_blender(obj). See the Blender docs.

ParaView plugin

If you have downloaded a binary version of ParaView, you may proceed as follows.

• Install meshio++ for the Python major version that ParaView uses (check pvpython --version)
• Open ParaView
• Find the file paraview-meshioplusplus-plugin.py of your meshio++ installation (on Linux: ~/.local/share/paraview-5.9/plugins/) and load it under Tools / Manage Plugins / Load New
• Optional: Activate Auto Load

You can now open all meshio++-supported files in ParaView.

Benchmarks

How much does the C++ core help? The benchmark/ folder times read/write conversions against the original pure-Python meshio on the formats both support (same in-memory mesh, same machine). The headline input is the bundled example.msh — a real Gmsh bracket (~52k nodes, ~293k cells).

meshio++'s biggest wins are the parallel and text paths: VTU binary+zlib ~16× write (the zlib blocks run across cores via an OpenMP backend with dynamic scheduling — hybrid P+E-core CPUs load-balance too), VTU ASCII ~7× write / ~5× read, and mixed-topology XDMF read ~10×. The binary and HDF5 formats that used to be slower — VTK/Gmsh binary, UGRID, and MED — are now at or above parity after an optimisation pass (bulk-buffered binary I/O, single-instruction bswap endianness conversion, a real parallel backend, an Eigen-backed MED transpose, zero-copy cell reconstruction that moves the connectivity buffer straight into the mesh, and uninitialised reader buffers + thread-parallel block copies so nothing is written twice); binary reads now match or beat numpy's fromfile — Gmsh ~1.7×, single-type VTK ~1.45×, and even mixed-topology VTK ~1.1×. Output stays byte-identical throughout.

The speedup is per-element: text/parallel formats climb out of the small-mesh regime and plateau (large meshes realise the full speedup):

Full methodology and a reproducible notebook are on the Benchmarks doc page (source: benchmark/01_benchmark.ipynb). The benchmark trends page keeps weekly history with raw measurements and environment metadata; timings remain informational. See fuzzing and sanitizers for the POSIX fuzz tooling, instrumented HDF5/netCDF campaign and the XDMF-bundle and optional-runtime (ADIOS2, CGNS MLL, TecIO) production-path campaigns, including fault-injection checks for mutator isolation. The OSS-Fuzz integration includes the confirmed maintainer contact and reproducible container checks for both declared sanitizers; upstream acceptance remains a separate onboarding step.

Contributors can audit the native core with tools/performance-tidy.sh --jobs 4: it runs Clang-Tidy's built-in performance-* checks separately from the include-hygiene gate, records deduplicated findings and actual build capabilities, and fails on incomplete scans without applying fixes. See the performance-audit workflow for ABI review and before/after validation.

The I/O sweep includes the single-type MFM tetrahedral input: python benchmark/bench.py --sizes M,L --formats mfm --repeats 7. Its native reader now uses bounded token views, with matching before/after file and mesh digests and substantially fewer allocation requests; see the measured comparison. MFM's mixed-cell conformance rejection and writer bytes are unchanged.

The measured text-I/O batches extend bounded token/source-line views across point clouds, scientific mesh formats and solver decks, preserving writer bytes and parsed data across SEQ/OpenMP/TBB. Evidence includes separate timing/allocation measurements, both SEQ measurement orders and retained unfavorable observations. Ownership tests cover source replacement, nested Marc includes and mapped/buffered paths. DEX's sweep input supplies its required scalar nodal field without mutating the shared input; Tecplot's removed unused helper has no runtime-speedup claim.

The shared-card pass retains the installed owning API while native readers use private views and bounded Fortran-number parsing. Fixed/free-record scope and neutral results are documented alongside the improvements; C++ ABI 22 is unchanged.

Eight format-private adapters also share bounded byte/record cursors, preserving grammar and binary positioning. Evidence explicitly exercises native MDPA (whose public Python reader remains the reference) and Gmsh 4.1's tested single-type input, without changing mixed-cell conformance.

The remaining deck-token pass extends views to Abaqus comma/data rows, Nastran logical cards, Netgen whitespace fields and FLAC3D ASCII tokens, with stable ownership for noncontiguous joins and include/gzip sources. Both SEQ orders confirm reader gains; writer bytes and ABI 22 are unchanged.

DEX normalization and source-copy removal preserve temporary-string ownership and add sanitizer-tested cursor guards. The evidence distinguishes measured text paths from PLY's binary timing input and retains unfavorable writer observations rather than claiming universal speedups.

Reading only what you need

import meshioplusplus

mesh = meshioplusplus.read("big.vtu", points_only=True)   # geometry, no data arrays
mesh = meshioplusplus.read("big.vtu", arrays=["u", "p"])  # only these arrays
meta = meshioplusplus.read_metadata("big.vtu")            # counts/names, no heavy arrays

mesh = meshioplusplus.read("run.exo", time_step=-1)       # the last step of a time series
meta["time_values"]                                       # how many steps there are

VTU, VTP, XDMF and Gmsh skip the unwanted array bodies outright; other formats are read in full and filtered, and meta["fell_back_to_full_read"] says which happened. time_step picks one step of a multi-step file (0 = the first, negative counts from the end); out of range is an error naming the available count rather than a silent fallback to step 0. Currently honoured by Exodus. A lenient option downgrades "this reader cannot represent construct X" errors to a warning plus a skip (currently MDPA's Table/Geometries/Mesh/Constraints blocks) — not "ignore all errors": a malformed file still fails. Large files can also be memory-mapped (automatic above 16 MiB), which roughly halves peak memory during a read. See selective reads and memory-mapped reading.

VTK XML output can additionally use lz4 (ParaView-readable) or zstd (a meshio++ extension) instead of zlib, when built with -DMESHIOPLUSPLUS_WITH_LZ4=ON / -DMESHIOPLUSPLUS_WITH_ZSTD=ON. zlib remains the default. See compression codecs.

The native reads, writes and operations release the GIL, so a batch of files converts in a thread pool, even on the sequential Linux wheels:

from concurrent.futures import ThreadPoolExecutor

def convert(path):
    meshioplusplus.write(path.replace(".msh", ".vtu"), meshioplusplus.read(path))

with ThreadPoolExecutor() as pool:
    list(pool.map(convert, paths))  # list() re-raises any worker's exception

The HDF5, netCDF, ADIOS2, TecIO and gidpost formats release it too, but those libraries are not thread-safe, so the core serialises them on one lock: their threads queue instead of running in parallel. See threads and the GIL.

Installation

meshio++ is available from the Python Package Index, so simply run

pip install meshioplusplus

to install.

Additional dependencies (netcdf4, h5py) are required for some of the output formats and can be pulled in by

pip install meshioplusplus[all]

For JavaScript / browser use, the C++ core also ships as a WebAssembly npm package covering every format above — including the HDF5- and netCDF-backed ones (CGNS, H5M, HMF, MED, Exodus), as of v8.0.0:

npm install @meshioplusplus/wasm

See the WebAssembly / JavaScript doc page for usage and the format-support table. Beside that lean artifact (MIT; Eigen is MPL-2.0) the package carries a full one (v17.2.0), loaded with loadMeshioPlusPlus({}, {full: true}): it links OpenCASCADE (STEP and IGES reading, fitCad), zstd, lz4 and bzip2, and is LGPL-2.1 with the OCCT exception for its OpenCASCADE part (lean and full builds).

meshio++ is also on Spack, upstream in spack/spack-packages — no checkout of this repo required:

spack install py-meshioplusplus +hdf5 +netcdf +zlib

See Installation → Spack for the standalone C API package and the full variant list.

C++ API

The full C++ core installs as a normal CMake package — the real Mesh, the format registry, every mesh/data operation, and the header-only Kratos bridge, none of which fit through a C ABI:

cmake -S . -B build -DMESHIOPLUSPLUS_BUILD_PYTHON=OFF -DMESHIOPLUSPLUS_INSTALL_CPP=ON
cmake --build build && cmake --install build --prefix /opt/meshioplusplus

find_package(meshioplusplus 17.23.0 EXACT CONFIG REQUIRED COMPONENTS CXX)
target_link_libraries(my_solver PRIVATE meshioplusplus::core)

EXACT is the conservative pin: the C++ API makes no ABI promise (Mesh, ModelPart and GeometricalEntity are header-defined types whose layout moves with the headers), so library and consumer must be built from compatible headers. The finer pin is MESHIOPLUSPLUS_ABI_VERSION, which moves only when a change really would break an already-compiled consumer — so a release that cannot affect you costs no rebuild. Either way a mismatch now fails at link time rather than corrupting memory, and the C API is the stable one — pin find_package(meshioplusplus 16 … COMPONENTS C) there. See ABI compatibility.

#include "meshioplusplus/registry.hpp"
#include "meshioplusplus/operations/partition.hpp"

auto mesh  = meshioplusplus::registry_read("bracket.msh", "", {});
auto parts = meshioplusplus::partition(mesh, {.mNParts = 8});

All three mesh backends install side by side (meshioplusplus::core_meshio, ::core_native, ::core_kratos), so one prefix serves consumers that disagree about the backend; each carries its own backend macro, making a mismatch a compile or link error rather than silent UB. Kratos consumers get the whole ModelPart surface: application entity names such as SmallDisplacementElement3D4N are preserved end to end (file → Mesh → ModelPart → file), material data crosses as Properties key/value pairs, and nested SubModelParts round-trip as parent/child region names. meshio++ itself is serial — there is no MPI anywhere in the API; partition(mesh, {nparts, ghost_layers}) produces the shared-node halo an MPI assembly needs and each rank takes its own piece. See the C++ API doc page.

C / Fortran API

For HPC codes written in C or Fortran, the C++ core also builds as an installable shared library (libmeshioplusplus, pure-C99 header, pkg-config + find_package support) with a modern OO Fortran 2008 module on top:

./build/configure.sh --fortran --tests --build     # --c-api for the C API alone
cmake --install build/cpp-release --prefix /opt/meshioplusplus

mio_mesh* m = mio_read("in.msh", NULL);
printf("%lld points\n", (long long)mio_mesh_num_points(m));
mio_write("out.vtu", m, NULL);
mio_mesh_free(m);

use meshioplusplus
type(mio_mesh) :: m
call m%read("in.msh")
call m%write("out.vtu")
call m%free()

The C API is also packaged for Conan (root conanfile.py) and vcpkg (overlay port under packages/vcpkg/meshioplusplus/), both driving the same install/find_package path, plus Spack (upstream in spack/spack-packages, no checkout needed):

conan create . -o meshioplusplus/*:with_hdf5=True
vcpkg install meshioplusplus --overlay-ports=ports
spack install meshioplusplus +fortran +hdf5

Full mesh access (build meshes from raw arrays, zero-copy readback) is covered on the C API and Fortran doc pages.

Julia / R bindings

The same installed C library also carries bindings for Julia and R, the two remaining languages of the scientific-computing audience. Both are layered on libmeshioplusplus exactly as the Fortran module is — no new C++, and the core stays untouched:

import MeshioPlusPlus as mio
m = mio.read("bracket.msh")
mio.write(mio.extract_surface(m), "surface.vtu")

library(meshioplusplus)
m <- mio_read("bracket.msh")
mio_write(mio_extract_surface(m), "surface.vtu")

Julia and R are both column-major, so — as in Fortran — points shaped (dim, n) and connectivity (nodes_per_cell, n) are the same memory as the C API's row-major shapes, and nothing is ever transposed. Node indices are 1-based, with the shift applied inside the copying accessors only; Julia additionally exposes genuine zero-copy borrows (points_ptr, connectivity_ptr) whose validity window is enforced rather than merely documented. R is copy-only — R vectors are R-managed, so a borrow cannot survive into R — and says so plainly instead of implying parity.

[!IMPORTANT] The Julia binding is not MIT. bindings/julia/ is released under the GNU General Public License, version 3 (GPL-3.0) — a copyleft license, not a permission-required one: anyone may use, modify or sell it commercially with no permission needed, but distributing it or a modified version of it must be under GPL-3.0 too, with source available; purely private use carries no obligation. Everything else in this repository, including the C API it calls and the R binding, remains MIT.

See the Julia and R doc pages.

Single-header C++

The whole C++ core is also amalgamated into one self-contained, STB-style header — src/single_include/meshioplusplus/meshioplusplus.hpp — with pugixml bundled and no external dependencies by default. Drop it in, no CMake or linking required:

// in exactly ONE .cpp:
#define MESHIOPLUSPLUS_IMPLEMENTATION
#include "meshioplusplus/meshioplusplus.hpp"
// elsewhere: just #include it (declarations only)

g++ -std=c++20 -I src/single_include main.cpp

It is generated by ./tools/amalgamate.sh and kept in sync by CI. See the single-header doc page (optional HDF5/netCDF/zlib formats via MESHIOPLUSPLUS_HAS_* macros).

example/cpp/ is the C++ counterpart of example/python/: the same tour of meshio++, called directly against this single header instead of the Python bindings, on the xeus-cpp Jupyter kernel — no PyVista either, renders go through meshio++'s own SVG writer. example/julia/ and example/r/ are the same tour again, called through the Julia and R bindings on their own Jupyter kernels (IJulia / IRkernel) — since those flat-ABI bindings can't drive the SVG writer's data-driven colouring either, quality/field renders are small charts instead of colour.

C++ mesh backends

Standalone C++ builds (no Python) can swap the in-memory mesh structure at compile time via MESHIOPLUSPLUS_MESH_BACKEND — every format works identically under each backend:

• MESHIO (default; the Python extension and PyPI wheels always use it) — mirrors the Python meshio.Mesh;
• NATIVE — the fastest pure-C++ structure (canonical Float64/Int64 storage, cell-type enum, CSR ragged blocks); the WebAssembly build uses it;
• KRATOS — a Kratos Multiphysics-style ModelPart (Nodes/Elements/Conditions/SubModelParts) plus a header-only templated bridge that populates a real Kratos::ModelPart with no Kratos build dependency.

./build/configure.sh --mesh-backend NATIVE --tests --build

All three also install side by side from a single prefix (MESHIOPLUSPLUS_INSTALL_CPP=ON), as meshioplusplus::core_meshio / ::core_native / ::core_kratos — see the C++ API page. See the C++ mesh backends doc page.

Testing

To run just the Python suite, check out this repository, install it with the test extras, and type

pytest tests/python/

To run everything — the Python suite alongside the C++ GoogleTest suite (and, when configured, the Fortran API test) — with a single command, pip install -e . first so meshioplusplus and pytest are importable, then configure with -DMESHIOPLUSPLUS_BUILD_TESTS=ON and let CTest discover both:

pip install -e .
cmake -S . -B build -DMESHIOPLUSPLUS_BUILD_TESTS=ON
cmake --build build
ctest --test-dir build --output-on-failure

ctest picks up the GoogleTest cases (via gtest_discover_tests) plus a python_suite entry that runs pytest tests/python/ — the latter is only added when meshioplusplus/pytest are importable by the Python CMake found, so a tree built without the editable install still runs the C++ suite alone.

License

meshio++ is published under the MIT license, with one exception: the Julia binding in bindings/julia/ is released under the GNU General Public License, version 3 (GPL-3.0). GPL-3.0 is copyleft, not permission-required: anyone may use, modify or sell it commercially without asking, but distributing it (or a modified version) must be under GPL-3.0 too, with source available; purely private/internal use carries no obligation at all. Nothing else is affected at the source level: the C++ core, the C API that binding calls, and the R binding are all MIT. The PyPI wheels, the release command-line binaries, the Blender extension zip (which carries the wheel) and the Conan and vcpkg packages built with OpenCASCADE, which is their default, contain LGPL-2.1 code, so they are licensed MIT AND (LGPL-2.1-only WITH OCCT-exception-1.0), as LICENSE-BINARIES.md and THIRD_PARTY_NOTICES.md state. A build from source without OpenCASCADE is MIT.

Acknowledgements

meshio++ is a fork of meshio by Nico Schlömer and its contributors (MIT). meshio is where the Mesh data model, the cell-type naming, and the great majority of the format readers and writers come from; this fork adds a C++20 core, an operations layer, and the C/Fortran/WebAssembly bindings on top of that foundation. The original copyright is retained in LICENSE.

Some code and test fixtures also come from Simvia's meshlane fork of meshio (MIT), credited per-change in CITATION.cff and CHANGELOG.md.

The C++ core has no mandatory dependency. Everything below is either optional, bundled, or confined to one binding or tool; the PyPI wheels and the release CLI binaries link OpenCASCADE, zlib, gidpost, pugixml and nlohmann/json statically (THIRD_PARTY_NOTICES.md).

Runtime dependencies (Python)

Project	Used for	License
NumPy	the array type the whole data model is built on	BSD-3-Clause
Rich	CLI output formatting	MIT

Optional dependencies

Project	Extra	Used for	License
Polyscope	[viewer], [all]	the desktop viewer and screenshot()	MIT
PyVista	[pyvista] / [interop]	to_pyvista() / from_pyvista()	BSD-3-Clause
trimesh	[trimesh] / [interop]	to_trimesh() / from_trimesh()	MIT
pyarrow	[arrow] / [interop]	the Arrow/Parquet data export	Apache-2.0
pandas	[pandas] / [interop]	to_pandas()	BSD-3-Clause
polars	[polars] / [interop]	to_polars()	MIT
zarr	[zarr]	the write_dataset(format="zarr") chunked-dataset backend	MIT
CuPy	— (wheels are CUDA-version-specific, e.g. cupy-cuda13x — see the GPU docs)	to_cupy() / from_cupy()	MIT
MCP Python SDK

Bundled and build-time

Project	Used for	License
pugixml (vendored, 1.14)	XML parsing in the C++ core (vendored in src/cpp/third_party/); statically linked into the wheels and CLI binaries	MIT
gidpost (vendored, 2.14)	the GiD postprocess writer; statically linked into the wheels, the CLI binaries and both WASM pairs	BSD-2-Clause-Views
nlohmann/json (git submodule, v3.12.0)	the JSON pipeline parser; header-only, statically linked into the wheels and CLI binaries	MIT
Eigen	the MED Fortran↔C transpose (git submodule, optional); compiled only where HDF5 is built, so only into the lean WASM build	MPL-2.0
pybind11	the Python bindings	BSD-3-Clause
scikit-build-core	the CMake-driven build backend	Apache-2.0
Emscripten	the WebAssembly build	MIT / NCSA
GoogleTest	the C++ test suite	BSD-3-Clause
zlib (vendored, 1.3.1)	statically linked into the wheels and CLI binaries, where it backs VTU compression and the GiD writer; a source build uses the system zlib	Zlib
Zstandard, LZ4

Browser viewer (isolated under src/viewer/)

Project	Used for	License
vtk.js	all rendering, colour maps and the scalar bar	BSD-3-Clause
Vite	the app build	MIT
vite-plugin-singlefile	the self-contained wheel-bundled build	MIT
Playwright	the end-to-end tests and the documentation screenshots	Apache-2.0

Also

Kratos Multiphysics (BSD-3-Clause) is the source of the skin-detection algorithm, the EnSight writer logic, the KaHIP partitioning approach and the FindKaHIP.cmake module, all these 3 implementation are from the original author of this project as well. Also thanks to its ModelPart design informs the KRATOS mesh backend. VTK and Verdict (both BSD-3-Clause) define the mesh-quality formulas. The documentation is built with VitePress (MIT) and Doxygen (GPL-2.0, used as a tool only). The logo renders the Stanford Bunny ("Stanford Bunny — Digitized!" by MakerBot, CC-BY).

Thank you to all of them.