---

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

| 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, `ConvectionCoe