PDE Toolbox — Full FEA Workflow

End-to-end finite element analysis: geometry → model setup → solve → post-process. Uses the modern femodel workflow (R2025a+).

When to Use

• Building geometry from primitives (multicuboid, multicylinder, multisphere), STL/STEP import, or 2-D decsg
• Setting up femodel with boundary conditions, loads, materials, and initial conditions
• Solving thermal, structural, or electromagnetic problems (steady, transient, modal, frequency, conduction)
• Post-processing FE results: interpolation, derived quantities, visualization

When Not to Use

• General-equation PDE (createpde(N)) — that legacy workflow is not covered here
• System-level simulation (Simulink/Simscape) — use product-specific skills
• Mesh-only tasks with no PDE solve (e.g., surface meshing for visualization)

Workflow Overview

1. Geometry — Create with primitives, decsg, file import, or boolean ops → wrap in fegeometry
2. Model — femodel(AnalysisType=..., Geometry=gm) → material, BCs, loads, ICs
3. Mesh — generateMesh(model) (default first, refine if needed)
4. Solve — result = solve(model) or solve(model, tlist)
5. Post-process — Extract fields, interpolate, compute derived quantities, visualize

Phase 1: Geometry

fegeometry — The Hub

gm = fegeometry(multicuboid(1, 1, 1));       % From primitives
gm = fegeometry("model.stl");                % From STL/STEP file
gm = fegeometry(decsg(gd, sf, ns));          % From 2-D CSG
gm = fegeometry(nodes, elements);            % From mesh data

fegeometry is for the femodel workflow only. Do NOT use fegeometry with createpde(N) — that legacy workflow uses geometryFromEdges (2-D) or importGeometry (3-D) instead. This skill covers femodel exclusively.

Key properties: NumCells, NumFaces, NumEdges, Vertices

3-D Primitives

Nested cells (vectors), stacked layers (ZOffset), hollow (Void=[true,false]):

gm = fegeometry(multicylinder([0.3, 0.5], 1, Void=[true, false]));  % hollow pipe
gm = fegeometry(multicuboid([1, 1], [1, 1], [0.3, 0.7], ZOffset=[0, 0.3]));  % stacked

See references/primitives-and-import.md for full options and file import details.

Boolean Operations

gmCombined = union(gm1, gm2);                         % Merge into 1 cell
gmCombined = union(gm1, gm2, KeepBoundaries=true);    % Preserve cells (multi-material)
gmCombined = union(gm1, gm2, KeepBoundaries=[true, false]);  % Selective per shape
gmResult = subtract(gm1, gm2);                        % Cut gm2 from gm1
gmResult = intersect(gm1, gm2);                       % Keep only overlapping region

KeepBoundaries: Use true when shapes get different materials (preserves internal faces as cell boundaries). Omit or use false to merge into a single cell.

Cell modification after boolean ops:

gm = mergeCells(gm);              % Merge ALL cells into one
gm = mergeCells(gm, [2, 3]);      % Merge specific cells (must be connected)
gm = deleteCell(gm, cellIDs);     % Remove unwanted cells

Rules: Union first, subtract last. The function is subtract — NOT subtractgeom. Never assemble pre-hollowed pieces. Call mergeCells only ONCE at the end.

See references/boolean-and-cell-ops.md for full strategy (Sculpt+Carve, Void flags, addCell, addVoid, face imprinting).

2-D Geometry with decsg

Each shape is a column vector in the geometry matrix. First entry identifies the type:

All columns must have the same row count — pad shorter ones with zeros. Set formula: + (union), - (subtract), * (intersect).

R1 = [3; 4; 0; 1; 1; 0; 0; 0; 0.5; 0.5];
C1 = [1; 0.5; 0.25; 0.15; 0; 0; 0; 0; 0; 0];
gd = [R1, C1]; sf = '(R1-C1)'; ns = char('R1', 'C1')';
gm = fegeometry(decsg(gd, sf, ns));

See references/decsg-and-2d-geometry.md for polygon vertices, polar-coordinate shapes, extrude, namespace rules.

Entity Identification

topFace = nearestFace(gm, [0, 0, 1]);       % Single point: row vector
faceIDs = nearestFace(gm, [0 0 1; 0 0 0]);  % Multiple points: N×3 matrix
frontEdge = nearestEdge(gm, [0.5, 0, 0.5]);
cellID = findCell(model.Geometry, [x, y, z]); % fegeometry only (not DiscreteGeometry)
facesOfCell = cellFaces(gm, 1);              % All faces of cell 1
facesOfCell = cellFaces(gm, 1, "external");  % Only outer boundary faces
edgesOfCell = cellEdges(gm, 1);
edgeIDs = faceEdges(gm, faceID);
fIDs = facesAttachedToEdges(gm, edgeID);              % Faces sharing an edge
fIDs = facesAttachedToEdges(gm, edgeID, "internal");  % Only internal faces (3-D)

Identify faces/edges BEFORE meshing. nearestVertex does not exist — use gm.Vertices + distance calculation.

Transforms

gm = translate(gm, [dx, dy, dz]);
gm = rotate(gm, angle);                          % angle° about z through origin
gm = rotate(gm, angle, [cx cy cz]);              % about z through point [cx,cy,cz]
gm = rotate(gm, angle, [x1 y1 z1], [x2 y2 z2]); % about LINE from pt1 to pt2
gm = scale(gm, [1, 1, -1]);                      % reflect across z=0 (use -1 on axis to flip)

rotate 4-arg: Both args are points defining the axis line, not direction+origin. E.g., about y through origin: rotate(gm, 90, [0 0 0], [0 1 0]). The 3-arg form ONLY rotates about z.

Extrude

gm3d = extrude(gm2d, [0.1, 0.3, 0.1]);   % 2-D → 3-D stacked layers along z
gm = extrude(gm, faceID, 0.2);            % 3-D face extrusion along outward normal

Mesh Generation

Only mesh when needed. Escalate: default → global → local.

gm = generateMesh(gm);                            % Default (try first)
gm = generateMesh(gm, Hmax=0.1, Hmin=0.01);      % Global control
gm = generateMesh(gm, HFace={[3,5], 0.02});       % Local face refinement (cell array!)
gm = generateMesh(gm, HEdge={edgeIDs, 0.005}, Hmax=0.2, Hgrad=1.5);
gm = generateMesh(gm, HVertex={vtxID, 0.005}, Hgrad=1.5);  % Stress concentrations

GeometricOrder: Default "quadratic". Use "linear" for faster solves. Important: 3-D magnetostatic/magneticHarmonic/electricHarmonic require linear mesh (Nedelec elements) — generateMesh handles this automatically.

Mesh properties: mesh.Nodes (3×N), mesh.Elements (connectivity), meshQuality(mesh) (0–1), area(mesh), volume(mesh).

Phase 2: Model Setup

Units: PDE Toolbox is unit-less. All inputs are pure numbers — the user must keep dimensions, properties, loads, and constants in a self-consistent unit system. Do not assume SI.

Create Model

model = femodel(AnalysisType="thermalSteady", Geometry=gm);

AnalysisType Reference

Pattern is <physics><Type>. NEVER "transientStructural" or "steadyStateThermal".

Material Properties

model.MaterialProperties = materialProperties(Material="steel");           % Catalog (preferred)
model.MaterialProperties = materialProperties(YoungsModulus=70e9, PoissonsRatio=0.33);  % Explicit
model.MaterialProperties(cellID) = materialProperties(Material="copper");  % Per-cell

Catalog names: "copper", "Invar", "steel", "aluminum", "brass", "tungsten", "iron", "gold", "silver", "lead", "zinc", "glass", "concrete", "wood".

Orthotropic (structural): All three as 3-element vectors — ShearModulus has no scalar form:

model.MaterialProperties = materialProperties(YoungsModulus=[Ex Ey Ez], ...
    PoissonsRatio=[nu_xy nu_yz nu_xz], ShearModulus=[Gxy Gyz Gxz], MassDensity=rho);

Note: Structural analysis is linear-elastic only — no plasticity or hyperelastic models.

Structural Damping

Three forms: Hysteretic (HystereticDamping=η on materialProperties, frequency response only), Rayleigh (model.DampingAlpha, model.DampingBeta, transient/frequency), Modal (DampingZeta=ζ at solve time with ModalResults).

% Hysteretic — frequency response
model.MaterialProperties = materialProperties(..., HystereticDamping=0.05);
% Rayleigh — transient/frequency (C = alpha*M + beta*K)
model.DampingAlpha = 10; model.DampingBeta = 0.002;
% Modal — with modal superposition
R = solve(model, tlist, ModalResults=Rm, DampingZeta=0.02);

See references/damping-reference.md for frequency-dependent damping and when to use each form.

Nonconstant (Function) Parameters

Any material property, BC, or load can be a function handle. Signatures: @(location, state) for BCs/loads/materials, @(location) for ICs only.

location: .x, .y, .z (always), .nx, .ny, .nz (boundary only). state: .time, .frequency, .u, .ux/.uy/.uz, .NormFluxDensity (analysis-dependent).

NaN convention (critical): Solver probes with NaN in state fields. Function MUST return NaN of correct size when state fields are NaN. Pure arithmetic propagates NaN automatically; only conditional logic needs an explicit check:

pressure = @(location, state) ...
    ifelse(isnan(state.time), NaN(1,numel(location.x)), ...
    1e5*sin(2*pi*100*state.time)*ones(1,numel(location.x)));

Vectorization: Output as 1×Np row (scalars) or 3×Np matrix (vectors). See references/nonconstant-parameters.md for state fields and output size rules. See references/time-dependent-loads.md for full NaN patterns with conditionals.

Boundary Conditions vs Loads

Rule: Prescribed values (Dirichlet) → BC. Everything else → Load.

Boundary dimension rule:

• 3-D: boundary = faces, edges, vertices (region = cell); body load = cellLoad
• 2-D: boundary = edges, vertices (region = face); body load = faceLoad

model.FaceBC(topFace) = faceBC(Temperature=100);          % 3-D thermal
model.EdgeBC(edgeID) = edgeBC(Temperature=100);            % 2-D thermal
model.VertexBC(vtxID) = vertexBC(Constraint="fixed");      % Structural pin
model.FaceLoad(sideFace) = faceLoad(ConvectionCoefficient=30, AmbientTemperature=25);
model.VertexLoad(tipVtx) = vertexLoad(Force=[0; -1000; 0]);  % Full vector, no YForce
model.CellLoad(1) = cellLoad(Heat=1e6);                    % Body load (3-D)
model.CellLoad(1) = cellLoad(Gravity=[0; 0; -9.81]);       % Self-weight (structural 3-D)
model.CellLoad(1) = cellLoad(Temperature=thermalResults);  % Thermal stress coupling (3-D)
model.FaceLoad(1) = faceLoad(Gravity=[0; -9.81]);          % Body loads in 2-D use faceLoad

vertexBC is structural-only: Constraint="fixed", XDisplacement, YDisplacement, ZDisplacement.

Common mistakes: Convection on faceBC → error (it's a load). faceBC in 2-D → error (use edgeBC). cellLoad in 2-D → error (use faceLoad).

See references/bc-load-reference.md for complete property tables.

Radiation

When using Emissivity on faceLoad/edgeLoad: set model.StefanBoltzmannConstant = 5.670374419e-8 and use all temperatures in Kelvin.

Thermal Stress

Temperature load generates strain only when ReferenceTemperature is set:

model.ReferenceTemperature = 20;  % Strain = CTE × (T - Tref)

Initial Conditions

Required for all transient analyses. Structural transient needs BOTH Displacement and Velocity.

model.CellIC = cellIC(Temperature=25);                                    % Thermal 3-D
model.FaceIC = faceIC(Temperature=25);                                    % Thermal 2-D
model.CellIC = cellIC(Displacement=[0;0;0], Velocity=[0;0;0]);           % Structural 3-D
model.FaceIC = faceIC(Displacement=[0;0], Velocity=[0;0]);               % Structural 2-D
model.CellIC = cellIC(MagneticVectorPotential=Rlin);                     % Nonlinear mag 3-D

PlanarType (2-D)

Set AFTER construction: model.PlanarType = "planeStress"; (or "planeStrain", "axisymmetric")

Model Reusability — ALWAYS Reuse

When a task requires multiple analysis types on the same geometry, change AnalysisType on the existing model — NEVER create a second femodel:

model = femodel(AnalysisType="structuralStatic", Geometry=gm);
model.MaterialProperties = materialProperties(Material="brass");
model.FaceBC(bottomFace) = faceBC(Constraint="fixed");
model = generateMesh(model);
staticResult = solve(model);
model.AnalysisType = "structuralModal";  % Reuse — geometry, material, mesh carry over
modalResult = solve(model, FrequencyRange=[-Inf, 5000]);
model.AnalysisType = "structuralTransient";  % Modal transient
transientResult = solve(model, tlist, ModalResults=modalResult, DampingZeta=0.02);

Multi-physics coupling sequences:

Mesh regeneration rule: When switching TO magnetostatic, magneticHarmonic, or electricHarmonic (Nedelec elements), the mesh must be linear. If the prior analysis used quadratic (default), call generateMesh(model) again — the solver will auto-select linear order for these types.

Electromagnetic Setup

EM types require vacuum constants:

model.VacuumPermittivity = 8.854187817e-12;   % Required for electrostatic, harmonic
model.VacuumPermeability = 4*pi*1e-7;         % Required for magnetostatic, harmonic

Nonlinear permeability (magnetostatic): Express μr as @(~,s) muRfcn(s.NormFluxDensity). Two-step: solve linear first with constant μr, then use result as initial guess via cellIC(MagneticVectorPotential=Rlin).

See references/electromagnetic-setup.md for CurrentDensity format (scalar in 2-D, vector in 3-D), nonlinear B-H workflow, harmonic EM, far-field BCs.

Phase 3: Solve

Units convention (critical): FrequencyRange, flist, and omega are all in rad/s, not Hz. Convert: omega = 2*pi*freqHz.

R = solve(model);                                          % Steady/static
R = solve(model, linspace(0, 2, 100));                     % Transient
R = solve(model, FrequencyRange=[-Inf, 2*pi*500]);         % Modal
R = solve(model, 2*pi*linspace(10, 1000, 200));            % Frequency response
Rm = solve(model, FrequencyRange=[-Inf, 2*pi*500]);        % Modal superposition...
R = solve(model, tlist, ModalResults=Rm, DampingZeta=0.02);
coarseR = solve(model, linspace(0, 40, 10));               % POD thermal...
R = solve(model, linspace(0, 40, 200), Snapshots=coarseR.Temperature);

SolverOptions: model.SolverOptions.RelativeTolerance, model.SolverOptions.AbsoluteTolerance.

Incremental range strategy (modal/decay): Never guess a large upper bound. Start small, check mode count, widen only if needed. Same for DecayRange. See references/modal-and-reduced-order.md for full details and modal transient workflow.

Phase 4: Post-Processing

Result Objects

Critical Rules

1. Stress.Magnitude is EMPTY — use result.VonMisesStress (static) or evaluateVonMisesStress(result) (transient)
2. evaluateHeatFlux requires 3 outputs: [qx, qy, qz] = evaluateHeatFlux(result)
3. Transient stress/strain are methods, not properties: evaluateVonMisesStress(result), evaluateStress(result)
4. FEStruct field names — ALWAYS use domain-specific names:

Interpolation

Query points are 3×N matrices (each column = one point):

Tq = interpolateTemperature(result, queryPts);        % Steady: no iT needed
Tq = interpolateTemperature(result, queryPts, iT);    % Transient: iT REQUIRED
d = interpolateDisplacement(result, queryPts);        % Returns all times (no iT accepted)

iT asymmetry: interpolateTemperature (transient) REQUIRES time index iT. interpolateDisplacement does NOT accept iT — returns all times. For intermediate times, use interp1 over SolutionTimes.

filterByIndex — Extract Time Steps

Required for transient visualization (pdeplot3D needs N×1, not N×T):

resultF = filterByIndex(result, numel(result.SolutionTimes));
vm = evaluateVonMisesStress(resultF);
pdeplot3D(resultF.Mesh, ColorMapData=vm, Deformation=resultF.Displacement, DeformationScaleFactor=100);

Not available on FrequencyStructuralResults — index columns directly:

iFreq = 10;
vm = evaluateVonMisesStress(result);  % N×F
pdeplot3D(result.Mesh, ColorMapData=vm(:,iFreq), ...
    Deformation=struct("ux", result.Displacement.ux(:,iFreq), ...
                       "uy", result.Displacement.uy(:,iFreq), ...
                       "uz", result.Displacement.uz(:,iFreq)), ...
    DeformationScaleFactor=500);

Visualization

pdeplot3D(result.Mesh, ColorMapData=result.Temperature, Mesh="on");
colorbar; title("Temperature");
pdeplot3D(result.Mesh, ColorMapData=result.VonMisesStress, ...
    Deformation=result.Displacement, DeformationScaleFactor=100);
[qx, qy, qz] = evaluateHeatFlux(result);
pdeplot3D(result.Mesh, FlowData=[qx qy qz]);
% Partial mesh — hide regions (e.g., air in EM)
elemIDs = findElements(result.Mesh, "region", Cell=[1 3]);
pdeplot3D(result.Mesh.Nodes, result.Mesh.Elements(:,elemIDs), ColorMapData=data);

See references/visualization-reference.md for pdeplot (2-D), pdeviz, FlowData, cross-sections, streamlines. See references/em-postprocessing.md for EM post-processing (Maxwell stress tensor, energy methods, all interpolation methods).

Common Mistakes

Copyright 2026 The MathWorks, Inc.