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

Function	Origin	Arguments
multicuboid(W, D, H)	x-y centered, base at z=0	Width, Depth, Height
multicylinder(R, H)	x-y centered, base at z=0	Radius, Height
multisphere(R)	Centered at origin	Radius

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:

Type code	Shape	Column format
1	Circle	[1; xc; yc; r; 0; ...]
2	Polygon	[2; N; x1;...;xN; y1;...;yN]
3	Rectangle	[3; 4; x1;x2;x3;x4; y1;y2;y3;y4] (CCW corners)
4	Ellipse	[4; xc; yc; a; b; angle; 0; ...]

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

AnalysisType	Use For
"thermalSteady" / "thermalTransient" / "thermalModal"	Heat transfer
"structuralStatic" / "structuralTransient" / "structuralModal" / "structuralFrequency"	Stress/displacement
"electrostatic" / "magnetostatic" / "dcConduction"	Static EM
"electricHarmonic" / "magneticHarmonic"	AC EM

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

BC (Dirichlet)	Load (Neumann/Robin)
faceBC, edgeBC, vertexBC	faceLoad, edgeLoad, cellLoad, vertexLoad
Temperature, Constraint, XDisplacement	Heat, ConvectionCoefficient+AmbientTemperature
Voltage, MagneticPotential	Pressure, SurfaceTraction, Force

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:

Sequence	Pattern	Caveat
static → modal → transient	Change AnalysisType, mesh reused	None
dcConduction → magnetostatic	cellLoad(CurrentDensity=conductionResults)	Must re-mesh: magnetostatic uses Nedelec (linear) elements; dcConduction defaults to quadratic
electrostatic → electricHarmonic	Change AnalysisType, add frequency	Add VacuumPermeability (harmonic needs both constants)
magnetostatic → magneticHarmonic	Change AnalysisType, add frequency	None if already linear mesh
thermalSteady → structuralStatic	cellLoad(Temperature=thermalResults)	Set ReferenceTemperature on structural model

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

Analysis type	Syntax	Notes
Steady/static	R = solve(model)	thermalSteady, structuralStatic, electrostatic, magnetostatic, dcConduction
Transient	R = solve(model, tlist)	tlist = monotonically increasing vector (seconds)
Modal (structural)	R = solve(model, FrequencyRange=[ω1, ω2])	Bounds in rad/s. Use -Inf lower bound
Modal (thermal)	R = solve(model, DecayRange=[-Inf, λmax])	Decay rates in s⁻¹. Use -Inf lower bound
Frequency response	R = solve(model, flist)	flist in rad/s (multiply Hz by 2π)
Modal superposition	R = solve(model, tlist, ModalResults=Rm)	Transient or frequency via prior modal solve
POD (thermal)	R = solve(model, tlist, Snapshots=Tmatrix)	Tmatrix = Temperature from coarse solve
Harmonic EM	R = solve(model, omega)	omega in rad/s

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

Analysis	Key Properties	Key Methods
thermalSteady	Temperature	evaluateHeatFlux, evaluateHeatRate, interpolateTemperature
thermalTransient	Temperature (N×T), SolutionTimes	Same + filterByIndex
structuralStatic	Displacement, VonMisesStress, Stress, Strain	evaluatePrincipalStress, evaluateReaction, interpolateDisplacement
structuralTransient	Displacement (N×T), SolutionTimes	evaluateVonMisesStress, evaluateStress, filterByIndex
structuralFrequency	Displacement (N×F), SolutionFrequencies	evaluateVonMisesStress, evaluateStress (no filterByIndex)
structuralModal	NaturalFrequencies, ModeShapes	Pass to solve(..., ModalResults=R)
electrostatic	ElectricPotential, ElectricField	interpolateElectricPotential, interpolateElectricField, interpolateElectricFlux
magnetostatic	MagneticFluxDensity, MagneticField	interpolateMagneticFlux, interpolateMagneticField
conduction	ElectricPotential, CurrentDensity	interpolateCurrentDensity

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:

Analysis	Fields	NEVER
Displacement	.ux, .uy, .uz	.x, .y, .z
Stress	.sxx, .syy, .szz, .syz, .sxz, .sxy	.xx, .yy
Strain	.exx, .eyy, .ezz, .eyz, .exz, .exy	.xx, .yy
Electric field	.Ex, .Ey, .Ez	.x, .y, .z
Magnetic flux density	.Bx, .By, .Bz	.x, .y, .z
Current density	.Jx, .Jy, .Jz	.x, .y, .z

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