Hi @SCIUSCIA
I understand that you are facing this error while using MATLAB PDE Toolbox:
Error using pde.ElectromagneticMaterialAssignmentRecords/findElectromagneticProperties
Missing material properties for one or more subdomains.
This typically means not all faces (subdomains) in your geometry have assigned material properties.
Here is how you can fix this error:
1. Visualize Face IDs
Use this to see all face numbers:
pdegplot(model, 'FaceLabels', 'on');
2. Assign Default Material to All Faces
Assign air (ε=1, μ=1) to all faces first:
nFaces = model.Geometry.NumFaces;
electromagneticProperties(model, 'Face', 1:nFaces, ...
'RelativePermittivity', 1, ...
'RelativePermeability', 1, ...
'Conductivity', 0);
3. Override Specific Faces with Custom Materials
electromagneticProperties(model, 'Face', SiO2_face, ...
'RelativePermittivity', eps_sio2);
electromagneticProperties(model, 'Face', Gold_face, ...
'RelativePermittivity', eps_au);
4. Specify the vacuum permeability
The electromagnetic model needs you to explicitly set both vacuum permittivity (ε₀) and vacuum permeability (μ₀). After creating the model, include:
model.VacuumPermittivity = eps0;
model.VacuumPermeability = 1.2566370614e-6;
5. Field Magnitude Calculation in 2D
For a 2D geometry, MATLAB populates only two components of the electric field (E.Ex, E.Ey). If you try to calculate the magnitude using all three (Ex, Ey, Ez), you may get an array size mismatch error.
Use only the existing components for your magnitude calculation:
Emag = sqrt(abs(E.Ex).^2 + abs(E.Ey).^2);
In summary, always visualize and assign materials to every face, set both vacuum permittivity and permeability properties in your model, and use the correct components for field calculations. These steps will systematically address all the common errors encountered after geometry and material assignment in electromagnetic simulations with MATLAB PDE Toolbox.
Here is the complete MATLAB code for your reference:
clear; close all; clc;
% ===== PHYSICAL PARAMETERS =====
lambda = 785e-9; % Wavelength [m]
omega = 2*pi*3e8/lambda; % Angular frequency [rad/s]
eps0 = 8.8541878128e-12; % Vacuum permittivity [F/m]
% ===== MATERIAL PROPERTIES =====
eps_au = (-24.051 + 1.3i) + 25; % Gold permittivity (shifted to Re(ε)>0)
eps_sio2 = 2.25; % SiO₂ relative permittivity
% ===== CREATE MODEL =====
model = createpde('electromagnetic','harmonic');
model.VacuumPermittivity = eps0;
model.VacuumPermeability = 1.2566370614e-6;
% ===== GEOMETRY DEFINITION =====
% SiO₂ substrate (500nm x 100nm)
rect_sio2 = [3;4; 0;500e-9;500e-9;0; 0;0;100e-9;100e-9];
% Air gap (500nm x 5nm, starting at y=100nm)
gap_y_start = 100e-9;
gap_height = 5e-9;
rect_air = [3;4; 0;500e-9;500e-9;0; ...
gap_y_start+1e-12; gap_y_start+1e-12; ...
gap_y_start+gap_height-1e-12; gap_y_start+gap_height-1e-12];
% Gold nanoparticle (radius 50nm at x=250nm, y=149nm)
circ = [1;250e-9;149e-9;50e-9; zeros(6,1)];
% Combine geometry with constructive solid geometry
gd = [rect_sio2, rect_air, circ];
ns = char('sio2','airgap','gold')';
sf = '(sio2 + airgap) + gold'; % Proper CSG formula
dl = decsg(gd,sf,ns);
geometryFromEdges(model, dl);
% ===== VISUALIZE FACE LABELS (CRUCIAL STEP) =====
figure;
pdegplot(model,'FaceLabels','on','EdgeLabels','on');
title('Geometry Face Labels');
axis equal;
% ===== MATERIAL ASSIGNMENT (UPDATED) =====
% Assign air (ε=1) to all faces by default
nFaces = model.Geometry.NumFaces;
electromagneticProperties(model, 'Face', 1:nFaces, ...
'RelativePermittivity', 1, ...
'RelativePermeability', 1, ...
'Conductivity', 0);
% --- Replace these face numbers (below) with your actual SiO₂ and Au face numbers ---
% Check the face IDs in the geometry plot, then update:
SiO2_face = 1; % Example: replace with SiO₂ substrate face number found from plot
Au_face = 3; % Example: replace with gold nanoparticle face number found from plot
electromagneticProperties(model, 'Face', SiO2_face, ...
'RelativePermittivity', eps_sio2, ...
'RelativePermeability', 1, ...
'Conductivity', 0);
electromagneticProperties(model, 'Face', Au_face, ...
'RelativePermittivity', eps_au, ...
'RelativePermeability', 1, ...
'Conductivity', 0);
% ---------------------------------------------------------------------
% ===== BOUNDARY CONDITIONS =====
k_spp = 2*pi/lambda * sqrt(eps_au/(1 + eps_au));
jsFun = @(location,~) exp(1i * k_spp * location.x);
electromagneticBC(model, 'Edge', 1, 'SurfaceCurrentDensity', jsFun); % SPP excitation
electromagneticBC(model, 'Edge', [2,4], 'FarField', 'absorbing', ...
'Thickness', lambda/4); % Absorbing BCs
% ===== MESH GENERATION =====
generateMesh(model, ...
'Hmax', lambda/20, ...
'Hmin', gap_height/2, ...
'Hgrad', 1.3, ...
'GeometricOrder', 'linear');
figure;
pdemesh(model);
title('Mesh');
axis equal;
% ===== SOLVE & VISUALIZE =====
result = solve(model, 'Frequency', omega);
E = result.ElectricField;
Emag = sqrt(abs(E.Ex).^2 + abs(E.Ey).^2);
figure;
pdeplot(model,'XYData',Emag,'ColorMap','jet','Mesh','off');
title('Electric Field Magnitude |E|');
colorbar;
axis equal;
I hope this helps.
