version: "1.0.1" name: matlab-solve-pde description: > End-to-end finite element analysis in MATLAB PDE Toolbox — geometry creation, model setup, solve, and post-processing in one skill. Use when building geometry from primitives or file import, setting up femodel with BCs/loads/materials, solving thermal/structural/EM problems, and extracting or visualizing results. Covers fegeometry, multicuboid, multicylinder, multisphere, decsg, boolean ops, mesh generation, femodel, all AnalysisTypes (thermalSteady, thermalTransient, structuralStatic, structuralTransient, structuralModal, structuralFrequency, electrostatic, magnetostatic, dcConduction, harmonic EM), materialProperties, faceBC, faceLoad, cellLoad, vertexLoad, solve, interpolation, von Mises stress, principal stress, reaction forces, heat flux, pdeplot3D visualization. Triggers on: PDE Toolbox, finite element, FEA, thermal analysis, structural analysis, electromagnetic analysis, femodel, mesh, boundary conditions, stress, displacement, heat transfer, post-processing. license: https://www.mathworks.com/content/dam/mathworks/license/pmrl/license.md metadata: author: MathWorks version: "1.0"

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.