2. Python ecosystem for CAE

A CAE automation workflow is rarely limited by solver speed. More often, it is limited by how efficiently engineering data can be moved between the solver, the pre/post-processor, and the optimization driver. Python's scientific stack provides a mature, well-tested foundation for that data movement. This chapter maps the five core libraries — NumPy, SciPy, pandas, and matplotlib, plus the file I/O layer that connects them — onto the concrete artifacts a CAE engineer actually handles: mesh node tables, load case matrices, solver input decks, result files, and convergence histories.

NumPy: the numerical substrate

NumPy's ndarray is the shared memory format that every other library in the stack either consumes or produces. For CAE work, the most important property is that a node coordinate table, a stiffness matrix block, or a time-history vector can all be represented as a single contiguous array with a known dtype and shape. This matters because solver integration almost always involves reshaping and slicing: extracting a subset of nodes by ID, reordering DOFs, or broadcasting a load vector across time steps.

import numpy as np

# Node table: columns = [id, x, y, z]
nodes = np.array([
    [1, 0.0, 0.0, 0.0],
    [2, 1.0, 0.0, 0.0],
    [3, 1.0, 1.0, 0.0],
    [4, 0.0, 1.0, 0.0],
], dtype=float)

coords = nodes[:, 1:]              # (4, 3) coordinate block
centroid = coords.mean(axis=0)     # (3,) geometric center
dist = np.linalg.norm(coords - centroid, axis=1)  # (4,) radial distance

Two conventions deserve early attention because they propagate through the entire course. First, array indexing is zero-based while CAE node and element IDs are typically one-based; the mapping must be explicit and consistent. Second, NumPy arrays are row-major by default, whereas many solver result files are written column-major or in Fortran order. Reading such files with the wrong order silently transposes the data — a failure mode that produces plausible-looking but incorrect results.

SciPy: solvers, interpolation, and optimization primitives

SciPy extends NumPy with the numerical algorithms that CAE workflows repeatedly need: linear algebra decompositions, sparse matrix handling, interpolation, root finding, and optimization. In a solver-integration context, scipy.sparse is often the most consequential module, because assembled stiffness and mass matrices are sparse by nature. Storing them as dense arrays wastes memory and time; csr_matrix or csc_matrix preserves the sparsity pattern and enables efficient matrix-vector products.

from scipy import sparse
from scipy.sparse.linalg import spsolve

# Example: 1D bar stiffness assembly (illustrative, not a real solve)
n = 5
K = sparse.lil_matrix((n, n))
for i in range(n - 1):
    k = 100.0  # element stiffness, arbitrary units
    K[i, i] += k
    K[i, i+1] -= k
    K[i+1, i] -= k
    K[i+1, i+1] += k
K = K.tocsr()

F = np.zeros(n); F[-1] = 10.0
u = spsolve(K[1:, 1:], F[1:])  # apply boundary condition at node 0

SciPy also provides scipy.optimize.minimize, which is frequently used as a lightweight surrogate for a full CAE optimization loop during development. The same objective function signature can later be wrapped around a solver call, so the transition from prototype to production is mostly a matter of replacing the analytical model with a file-based one.

pandas: tabular engineering data

CAE data is often tabular: load case definitions, material property libraries, sensor channels, or DOE matrices. pandas provides labeled axes, mixed dtypes, and grouping operations that make these tables manageable without manual index bookkeeping. A typical use is reading a CSV of load cases, filtering by analysis type, and iterating over the resulting rows to generate solver input decks.

import pandas as pd

load_cases = pd.read_csv("load_cases.csv")
# columns: case_id, fx, fy, fz, analysis_type

static = load_cases[load_cases["analysis_type"] == "static"]
for _, row in static.iterrows():
    forces = row[["fx", "fy", "fz"]].to_numpy(dtype=float)
    # forces is now a (3,) array ready for the solver interface

The critical discipline with pandas is dtype awareness. A column that looks numeric may be read as object if it contains a stray unit string or a blank cell. Calling .to_numpy(dtype=float) will then raise, which is preferable to silently propagating bad values. Use pd.to_numeric(..., errors="coerce") and inspect isna() when ingesting externally produced tables.

matplotlib: convergence and result visualization

In an optimization loop, the most valuable plot is usually the convergence history: objective value and constraint violation versus iteration. matplotlib makes this a few lines, and the same figure can be regenerated automatically after every run.

import matplotlib.pyplot as plt

history = np.array([
    [0, 1.000, 0.250],
    [1, 0.812, 0.140],
    [2, 0.744, 0.061],
    [3, 0.731, 0.018],
    [4, 0.729, 0.004],
])  # columns: iteration, objective, max_constraint_violation

fig, ax1 = plt.subplots()
ax1.plot(history[:, 0], history[:, 1], "o-", color="#1f77b4", label="objective")
ax1.set_xlabel("Iteration"); ax1.set_ylabel("Objective")
ax2 = ax1.twinx()
ax2.plot(history[:, 0], history[:, 2], "s--", color="#d62728", label="max violation")
ax2.set_ylabel("Constraint violation")
fig.tight_layout(); fig.savefig("convergence.png", dpi=150)

For CAE, always set explicit axis labels and units, and prefer saving to PNG or SVG rather than relying on interactive display, so that reports are reproducible in headless environments.

File I/O: the connective tissue

The libraries above are only useful if engineering data can enter and leave them reliably. CAE file formats fall into three practical categories, each with a different Python strategy.

For solver-native text, a robust parser reads line by line, strips comments, and converts fixed-width or delimited fields explicitly. The following pattern handles a simple node block and is representative of the parsing logic used throughout the course.

import re

def parse_nodes(path):
    pattern = re.compile(r"^\s*(\d+)\s+([-\d.eE+]+)\s+([-\d.eE+]+)\s+([-\d.eE+]+)")
    rows = []
    with open(path, "r") as f:
        for line in f:
            if line.lstrip().startswith(("#", "$", "!")):
                continue
            m = pattern.match(line)
            if m:
                rows.append([int(m.group(1))] + [float(m.group(i)) for i in (2, 3, 4)])
    return np.array(rows, dtype=float)

Two failure symptoms are worth recognizing early. If the returned array has fewer rows than expected, the regex is likely rejecting lines with leading whitespace or scientific notation variants. If coordinates are swapped or negated, the column order or sign convention in the source file differs from the assumed one — always verify against a known node before trusting the parser.

The diagram above summarizes the data path that this chapter establishes. Solver and pre-processor files enter through a dedicated I/O layer, are normalized into NumPy arrays or pandas DataFrames, are processed by SciPy routines, and are visualized with matplotlib. The dashed feedback path is the optimization loop itself: results are evaluated against constraints, and the design variables are updated and written back to the solver. Every subsequent chapter in this course builds on this pipeline, so the conventions chosen here — explicit dtypes, one-based to zero-based ID mapping, a single I/O module, and reproducible figure output — should be treated as project standards rather than optional style choices.