2. Python ecosystem for CAE post-processing

CAE post-processing is fundamentally a data problem disguised as a visualization problem. A solver writes displacement, stress, pressure, or temperature fields to disk; the engineer must read those fields, associate them with mesh topology, compute derived quantities, and present them in a form that supports a decision. Python's scientific stack covers this entire chain, but each library occupies a distinct layer. Understanding which layer does what prevents the common failure mode of forcing one tool to do another's job — for example, using pandas to hold a 50-million-cell field array, or using PyVista to compute a group-by aggregation.

The Layered Stack

The ecosystem organizes naturally into four layers, each with a clear responsibility boundary:

• Array layer — NumPy. Owns dense, homogeneous, N-dimensional numeric data. A nodal displacement array of shape (n_nodes, 3) and a connectivity array of shape (n_cells, 8) are NumPy arrays. NumPy provides the vectorized arithmetic, indexing, and reduction operations that make field manipulation fast.
• Tabular layer — pandas. Owns heterogeneous, labeled, row-oriented data. Time histories, probe-point extractions, load-step summaries, and result tables with mixed dtypes belong here. pandas gives you named columns, alignment, grouping, and export to CSV or Excel.
• Mesh I/O layer — meshio. Owns format translation. It reads and writes VTK, Abaqus, Nastran, Gmsh, Exodus, CGNS, and dozens of other mesh formats, exposing a uniform Mesh object with points, cells, and point_data/cell_data dictionaries.
• Visualization layer — VTK and PyVista. Owns rendering and geometric processing. VTK is the underlying C++ engine; PyVista is a Pythonic wrapper that exposes VTK's data model through NumPy-aware classes such as PolyData and UnstructuredGrid.

The dependency direction matters: pandas and PyVista both build on NumPy; PyVista builds on VTK; meshio is largely independent and serves as the bridge between solver output and the visualization layer. A well-structured post-processing script respects this direction and does not, for instance, route mesh geometry through a DataFrame.

NumPy: The Field Arithmetic Engine

Every downstream library ultimately hands you NumPy arrays, so fluency with array semantics is the highest-leverage skill. Consider a typical operation: extracting the von Mises stress from a six-component symmetric stress tensor stored as a (n_cells, 6) array in Voigt order [sxx, syy, szz, sxy, syz, sxz]. The computation is a direct vectorized expression:

import numpy as np

# stress: (n_cells, 6) in Voigt order [sxx, syy, szz, sxy, syz, sxz]
sxx, syy, szz = stress[:, 0], stress[:, 1], stress[:, 2]
sxy, syz, sxz = stress[:, 3], stress[:, 4], stress[:, 5]

vm = np.sqrt(0.5 * ((sxx - syy)**2 + (syy - szz)**2 + (szz - sxx)**2)
             + 3.0 * (sxy**2 + syz**2 + sxz**2))

Two details deserve attention. First, the Voigt ordering convention is not universal — some solvers store shear components as [sxy, sxz, syz] or apply engineering shear strains. Always confirm the ordering from the solver documentation before writing the extraction. Second, the expression is fully vectorized: no Python loop over cells, so a 10-million-cell array completes in well under a second. Replacing this with a per-cell loop would inflate runtime by two to three orders of magnitude.

NumPy also handles the geometric operations that precede visualization: computing cell centroids via connectivity-indexed point averaging, calculating bounding boxes with points.min(axis=0) and points.max(axis=0), and masking regions of interest with boolean index arrays. These operations produce the derived fields that later get attached to the mesh for rendering.

pandas: Structured Result Tables

Where NumPy excels at dense fields, pandas excels at the sparse, labeled, mixed-type data that surrounds them. A transient analysis with 200 load steps and 15 probe points yields 3,000 rows — small enough for a DataFrame, and naturally expressed with columns step, time, probe_id, displacement_mag, von_mises_max. The idiomatic workflow reads solver-exported CSV, coerces types, and aggregates:

import pandas as pd

df = pd.read_csv("probe_history.csv", parse_dates=["time"])
summary = (df.groupby("probe_id")
             .agg(peak_vm=("von_mises_max", "max"),
                  peak_step=("von_mises_max", "idxmax"),
                  mean_disp=("displacement_mag", "mean"))
             .reset_index())

The critical discipline is to keep pandas at the tabular layer. A DataFrame with one row per cell and a column per field component is a legitimate result table; a DataFrame holding the full connectivity array is an anti-pattern that destroys performance and obscures intent. When a field must move from NumPy to pandas, it should be reduced first — for example, extracting the top 1% of cells by stress rather than the entire field.

meshio: Format-Agnostic Mesh I/O

meshio's value is uniformity. Whether the source is an Abaqus .inp, a Gmsh .msh, an Exodus .exo, or a VTK .vtu, the result is a Mesh object with the same attribute names:

import meshio

mesh = meshio.read("bracket_result.vtu")
print(mesh.points.shape)          # (n_nodes, 3)
print(mesh.cells[0].type)         # e.g. 'hexahedron'
print(mesh.point_data.keys())     # dict_keys(['displacement', 'temperature'])
print(mesh.cell_data.keys())      # dict_keys(['von_mises', 'eq_plastic_strain'])

This uniformity means a conversion pipeline can be written once and applied across solver formats. A common pattern is to read a native solver format, attach a derived field computed in NumPy, and write a .vtu for PyVista to consume. Note that meshio distinguishes point_data (nodal fields) from cell_data (element fields); mixing them up produces either interpolation artifacts or outright shape errors, so verify array lengths against len(mesh.points) and the total cell count respectively.

PyVista and VTK: Geometry, Processing, and Rendering

PyVista wraps VTK's data model in classes that behave like NumPy containers. An UnstructuredGrid built from meshio output supports the full range of engineering visualization operations — contouring, clipping, warping by a displacement vector, thresholding by a scalar range — and renders either interactively or off-screen for batch report generation:

import pyvista as pv

grid = pv.UnstructuredGrid("bracket_result.vtu")
warped = grid.warp_by_vector("displacement", factor=50.0)
contour = warped.contour(isosurfaces=8, scalars="von_mises")

pl = pv.Plotter(off_screen=True)
pl.add_mesh(contour, cmap="turbo", show_scalar_bar=True)
pl.add_mesh(warped, opacity=0.15, color="gray")
pl.camera_position = "xy"
pl.screenshot("stress_contour.png", window_size=[1600, 1000])

The factor=50.0 is a display exaggeration, not a physical scaling — a point worth stating explicitly in any report caption, since exaggerated deformation is easily misread as actual displacement. Off-screen rendering with off_screen=True is what makes automated report generation feasible: no display server is required, and the same script runs in CI or on a headless compute node.

Choosing the Right Tool

The practical takeaway is that these libraries are complementary, not competing. A robust post-processing pipeline reads with meshio, computes derived fields with NumPy, summarizes with pandas, and renders with PyVista — each layer doing what it was designed for, with NumPy arrays as the common currency flowing between them.