01 — Mesh I/O¶

How to get meshes in and out of imatools, and what you get back.

By the end you will be able to:

• build a mesh without needing any data files
• save and load meshes with io.mesh_io
• read a MeshContract and know what is actually inside it
• drop to raw VTK and to NumPy when you need to
• round-trip a mesh through the CARP .pts/.elem format

This notebook is self-contained. It generates its own mesh and writes to a temporary directory, so it runs anywhere with no external data.

Setup¶

Paths use pathlib.Path throughout — that is the convention for new code.

In [ ]:

import tempfile
from pathlib import Path

import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy

from imatools.io.mesh_io import load_mesh, save_mesh
from imatools.io import carp_io

work = Path(tempfile.mkdtemp())
work

import tempfile from pathlib import Path import numpy as np import vtk from vtk.util.numpy_support import vtk_to_numpy from imatools.io.mesh_io import load_mesh, save_mesh from imatools.io import carp_io work = Path(tempfile.mkdtemp()) work

A mesh from nothing¶

VTK ships procedural sources, so a tutorial never needs a data file. Here is a sphere as a triangulated surface (vtkPolyData) — the mesh type most of imatools works with.

In [ ]:

src = vtk.vtkSphereSource()
src.SetRadius(10.0)
src.SetThetaResolution(12)
src.SetPhiResolution(12)
src.Update()

sphere = src.GetOutput()
print(f"{sphere.GetNumberOfPoints()} points, {sphere.GetNumberOfCells()} cells")

src = vtk.vtkSphereSource() src.SetRadius(10.0) src.SetThetaResolution(12) src.SetPhiResolution(12) src.Update() sphere = src.GetOutput() print(f"{sphere.GetNumberOfPoints()} points, {sphere.GetNumberOfCells()} cells")

Saving¶

save_mesh takes either a raw VTK object or a MeshContract. overwrite=False is the default and it will refuse to clobber an existing file — pass overwrite=True deliberately.

In [ ]:

mesh_path = work / "sphere.vtk"
save_mesh(sphere, mesh_path, mesh_type="polydata", overwrite=True)

print(mesh_path.name, "->", mesh_path.exists(), f"({mesh_path.stat().st_size} bytes)")

mesh_path = work / "sphere.vtk" save_mesh(sphere, mesh_path, mesh_type="polydata", overwrite=True) print(mesh_path.name, "->", mesh_path.exists(), f"({mesh_path.stat().st_size} bytes)")

Loading: the contract¶

load_mesh returns a MeshContract by default. A contract is the typed object that moves between layers: it carries the mesh plus metadata computed at load time, so callers do not each re-derive point counts and bounds.

In [ ]:

contract = load_mesh(mesh_path, mesh_type="polydata")

print("type:      ", type(contract).__name__)
print("points:    ", contract.metadata.num_points)
print("cells:     ", contract.metadata.num_cells)
print("mesh_type: ", contract.metadata.mesh_type)
print("bounds:    ", tuple(round(b, 1) for b in contract.metadata.bounds))
print("path:      ", contract.path.name)

contract = load_mesh(mesh_path, mesh_type="polydata") print("type: ", type(contract).__name__) print("points: ", contract.metadata.num_points) print("cells: ", contract.metadata.num_cells) print("mesh_type: ", contract.metadata.mesh_type) print("bounds: ", tuple(round(b, 1) for b in contract.metadata.bounds)) print("path: ", contract.path.name)

Gotcha. MeshContract has both polydata and points fields, but load_mesh only fills the VTK object — contract.points is None. The NumPy fields are for contracts built from arrays. To get NumPy from a loaded mesh, go through the VTK object (below).

In [ ]:

print("polydata:", type(contract.polydata).__name__)
print("points:  ", contract.points)  # None — not populated by load_mesh

print("polydata:", type(contract.polydata).__name__) print("points: ", contract.points) # None — not populated by load_mesh

Escape hatch: raw VTK¶

Contracts are the default, not a cage. Pass return_contract=False to get the plain VTK object — useful when handing off to VTK's own filters.

In [ ]:

raw = load_mesh(mesh_path, mesh_type="polydata", return_contract=False)
print(type(raw).__name__, "with", raw.GetNumberOfPoints(), "points")

raw = load_mesh(mesh_path, mesh_type="polydata", return_contract=False) print(type(raw).__name__, "with", raw.GetNumberOfPoints(), "points")

Mesh to NumPy¶

vtk_to_numpy turns VTK arrays into NumPy views — this is how you get from a mesh to something you can do maths on.

In [ ]:

points = vtk_to_numpy(raw.GetPoints().GetData())

print("shape:", points.shape, "dtype:", points.dtype)
print("centroid:", points.mean(axis=0).round(3))

points = vtk_to_numpy(raw.GetPoints().GetData()) print("shape:", points.shape, "dtype:", points.dtype) print("centroid:", points.mean(axis=0).round(3))

CARP round trip¶

CARP stores a mesh as two text files: .pts (vertices) and .elem (connectivity). saveToCarpTxt takes a points array and an element array and writes both, given a basename with no extension.

In [ ]:

elements = []
for i in range(raw.GetNumberOfCells()):
    ids = raw.GetCell(i).GetPointIds()
    elements.append([ids.GetId(j) for j in range(ids.GetNumberOfIds())])
elements = np.array(elements)

print("element array:", elements.shape, "(triangles)")

elements = [] for i in range(raw.GetNumberOfCells()): ids = raw.GetCell(i).GetPointIds() elements.append([ids.GetId(j) for j in range(ids.GetNumberOfIds())]) elements = np.array(elements) print("element array:", elements.shape, "(triangles)")

In [ ]:

basename = str(work / "sphere")  # no extension
carp_io.saveToCarpTxt(points, elements, basename)

print(sorted(p.name for p in work.glob("sphere.*")))

basename = str(work / "sphere") # no extension carp_io.saveToCarpTxt(points, elements, basename) print(sorted(p.name for p in work.glob("sphere.*")))

Read it back and confirm the round trip preserved the geometry.

In [ ]:

pts_back = np.asarray(carp_io.read_pts(basename + ".pts"))

print("read back:", pts_back.shape)
print("matches original:", np.allclose(pts_back, points, atol=1e-5))

pts_back = np.asarray(carp_io.read_pts(basename + ".pts")) print("read back:", pts_back.shape) print("matches original:", np.allclose(pts_back, points, atol=1e-5))

Recap¶

Task	Call
Load a mesh (typed)	load_mesh(path) -> MeshContract
Load a mesh (raw VTK)	load_mesh(path, return_contract=False)
Save a mesh	save_mesh(obj, path, overwrite=True)
Mesh -> NumPy	vtk_to_numpy(obj.GetPoints().GetData())
Write CARP	carp_io.saveToCarpTxt(pts, el, basename)
Read CARP points	carp_io.read_pts(basename + ".pts")

Next: 02_scalar_fields — point vs cell data, and projecting scalars onto a mesh.