- trimesh (to plot a triangular mesh) - https://www.mathworks.com/help/matlab/ref/trimesh.html
- createpde (for creating an electromagnetic model) - https://www.mathworks.com/help/pde/ug/createpde.html
- electromagneticBC (define boundary conditions) - https://www.mathworks.com/help/pde/ug/pde.electromagneticmodel.electromagneticbc.html
- electromagneticSource (specify current density, charge density, or magnetization values) - https://www.mathworks.com/help/pde/ug/pde.electromagneticmodel.electromagneticsource.html
Importing geometry from cell data for 3D model
9 views (last 30 days)
Show older comments
Hi,
I wonder if anyone may be able to help with this problem. I'm trying to create a 3D Finite Difference Time Domain (FDTD) simulation of a wireless communications link. I have put together the following program (see below).
The first part of the program imports the geometry from Openstreetmap using the uavScenario tool. This seems to work fine and produces the following showing the buildings on a ground plane:
From this plot I would like to import the geometry into my 3D FDTD simulator.
I have noticed on exploration of the file uavScenario generates, that it produces a cell called under the structure 'scene.Meshes' which appears to contain the coordinates of the vertices and faces of all of the 27 buildings and the ground plane as seen in the plot above.
In order to carry out the 3D FDTD simulation, I wish to import the above buildings etc. into the simulation space and assign dielectric conditions in order to simulate the material properties e.g. 'concrete'. Then, in theory, I should be able to observe the effects of how a plane wave from a source, propagates through the map geometry.
The 3D FDTD code with PML boundaries appears to work fine, as I'm able to produce a sensible simulation with a test signal as below:
I hope someone may be able to advise how this problem may be resolved (if it can), as no doubt it would be useful in many applications, not just a 3D FDTD simulation. I have attached the Map file (I had to put it into a .zip format as attachments would not except .osm file format)
Thanking you in advance.
Best regards,
Andy
My code so far, is as follows. I have also attached the Openstreetmap file I generated if anyone wants to explore.
%Import Map from Openstreetmap using uavScenario
scene = uavScenario(ReferenceLocation=[52.76269 -1.23959 110],UpdateRate=5);
xTerrainLimits = [-200 200];
yTerrainLimits = [-200 200];
color = [0.6 0.6 0.6];
addMesh(scene,"terrain",{"gmted2010",xTerrainLimits,yTerrainLimits},color)
xBuildingLimits = [-150 150];
yBuildingLimits = [-150 150];
color = [0.6431 0.8706 0.6275];
addMesh(scene,"buildings",{"LU.osm",xBuildingLimits,yBuildingLimits,"auto"},color)
figure(1)
show3D(scene)
Xscale=max(xTerrainLimits);
Yscale=max(yTerrainLimits);
%% Define Simulation Based on Source and Wavelength
f0 = 2.45e9; % Frequency of Source [Hertz]
Lf = 1; % Divisions per Wavelength [unitless]
[Lx,Ly,Lz] = deal(Xscale,Yscale,3); % Wavelengths x,y,z [unitless]
nt = 500; % Number of time steps (simulation length) [unitless]
%% PML Parameters
pml_thickness = 32; % Thickness of PML layers (in grid points)
pml_alpha = 0.01; % PML absorption coefficient
%% Simulation Performance
single = 1; % Use Single Precision | 1 = yes , 0 = no |
usegpu = 1; % Use gpuArray | 1 = yes , 0 = no |
%% Spatial and Temporal System
e0 = 8.854*10^-12; % Permittivity of vacuum [farad/meter]
u0 = 4*pi*10^-7; % Permeability of vacuum [henry/meter]
c0 = 1/(e0*u0)^.5; % Speed of light [meter/second]
L0 = c0/f0; % Freespace Wavelength [meter]
t0 = 1/f0; % Source Period [second]
%% Geometry
% This is where I would like to import the geometry to
%% Add PML Boundaries
[Nx,Ny,Nz] = deal(Lx*Lf,Ly*Lf,Lz*Lf); % Points in x,y,z [unitless]
x = linspace(0,Lx,Nx+1)*L0; % x vector [meter]
y = linspace(0,Ly,Ny+1)*L0; % y vector [meter]
z = linspace(0,Lz,Nz+1)*L0; % z vector [meter]
% Extend domain with PML layers
Nx_pml = Nx + 2*pml_thickness;
Ny_pml = Ny + 2*pml_thickness;
Nz_pml = Nz + 2*pml_thickness;
% Update dimensions and spacing for the extended domain
dx = L0 / Nx;
dy = L0 / Ny;
dz = L0 / Nz;
dt = min([dx, dy, dz]) / (sqrt(3) * c0) * 0.99; % CFL-condition time step
% Create coordinates for the extended domain
x_pml = linspace(-pml_thickness*dx, Lx+dx+pml_thickness*dx, Nx_pml+1);
y_pml = linspace(-pml_thickness*dy, Ly+dy+pml_thickness*dy, Ny_pml+1);
z_pml = linspace(-pml_thickness*dz, Lz+dz+pml_thickness*dz, Nz_pml+1);
%% Initialize Magnetic and Electric Field Vectors and Coefficients
[Ex] = deal(zeros(Nx_pml,Ny_pml+1,Nz_pml+1)); % Ex cells
[Ey] = deal(zeros(Nx_pml+1,Ny_pml,Nz_pml+1)); % Ey cells
[Ez] = deal(zeros(Nx_pml+1,Ny_pml+1,Nz_pml)); % Ez cells
[Hx] = deal(zeros(Nx_pml+1,Ny_pml,Nz_pml)); % Hx cells
[Hy] = deal(zeros(Nx_pml,Ny_pml+1,Nz_pml)); % Hy cells
[Hz] = deal(zeros(Nx_pml,Ny_pml,Nz_pml+1)); % Hz cells
[udx,udy,udz] = deal(dt/(u0*dx),dt/(u0*dy),dt/(u0*dz)); % H Coeffcients
[edx,edy,edz] = deal(dt/(e0*dx),dt/(e0*dy),dt/(e0*dz)); % E Coeffcients
%% Performance Enhancement
if single == 1; vars2Single(); end
if usegpu == 1; gpu2Single(); end
%% Start loop
for t = 1:nt
%% Magnetic Field Update
Hx = Hx + udz .* diff(Ey, 1, 3) - udy .* diff(Ez, 1, 2);
Hy = Hy + udz .* diff(Ez, 1, 1) - udx .* diff(Ex, 1, 3);
Hz = Hz + udy .* diff(Ex, 1, 2) - udx .* diff(Ey, 1, 1);
%% Electric Field Update
Ex(:,2:end-1,2:end-1) = Ex(:,2:end-1,2:end-1) + ...
edy .* diff(Hz(:,:,2:end-1), 1, 2) - ...
edz .* diff(Hy(:,2:end-1,:), 1, 3);
Ey(2:end-1,:,2:end-1) = Ey(2:end-1,:,2:end-1) + ...
edz .* diff(Hx(2:end-1,:,:), 1, 3) - ...
edx .* diff(Hz(:,:,2:end-1), 1, 1);
Ez(2:end-1,2:end-1,:) = Ez(2:end-1,2:end-1,:) + ...
edx .* diff(Hy(:,2:end-1,:), 1, 1) - ...
edy .* diff(Hx(2:end-1,:,:), 1, 2);
%% PML Absorption
% Update fields in PML regions
Ex = pml_absorption(Ex, Nx_pml, Ny_pml, Nz_pml, pml_thickness, pml_alpha, dx, dy, dz);
Ey = pml_absorption(Ey, Nx_pml, Ny_pml, Nz_pml, pml_thickness, pml_alpha, dx, dy, dz);
Ez = pml_absorption(Ez, Nx_pml, Ny_pml, Nz_pml, pml_thickness, pml_alpha, dx, dy, dz);
%% Point Source
Ez(round(Nx_pml/2)+1,round(Ny_pml/2)+1,round(Nz_pml/2)) = ...
sin(2*pi*f0*dt*t)./(1+exp(-.2*(t-60)));
%% Plotting
if t == 1
[X,Y,Z] = meshgrid(y_pml,x_pml,z_pml(1:end-1)); % for plotting Ez
end
figure(2)
slice(X,Y,Z,Ez,y_pml(end)/2,x_pml(end)/2,z_pml(end)/2);
axis([0 y_pml(end) 0 x_pml(end) 0 z_pml(end)]);
view([1 1 1]); caxis([-.001 .001]); shading interp; drawnow;
end
%% PML Absorption Function
function F = pml_absorption(F, Nx_pml, Ny_pml, Nz_pml, pml_thickness, pml_alpha, dx, dy, dz)
% Apply PML absorption
pml_profile = @(n, thickness) pml_alpha * ((thickness - n + 0.5) / thickness)^2;
% Apply PML to x direction
for i = 1:pml_thickness
alpha_x = pml_profile(i, pml_thickness);
F(i,:,:) = F(i,:,:) * exp(-alpha_x * (pml_thickness - i) * dx);
F(end-i+1,:,:) = F(end-i+1,:,:) * exp(-alpha_x * (pml_thickness - i) * dx);
end
% Apply PML to y direction
for j = 1:pml_thickness
alpha_y = pml_profile(j, pml_thickness);
F(:,j,:) = F(:,j,:) * exp(-alpha_y * (pml_thickness - j) * dy);
F(:,end-j+1,:) = F(:,end-j+1,:) * exp(-alpha_y * (pml_thickness - j) * dy);
end
% Apply PML to z direction
for k = 1:pml_thickness
alpha_z = pml_profile(k, pml_thickness);
F(:,:,k) = F(:,:,k) * exp(-alpha_z * (pml_thickness - k) * dz);
F(:,:,end-k+1) = F(:,:,end-k+1) * exp(-alpha_z * (pml_thickness - k) * dz);
end
end
0 Comments
Accepted Answer
Alan
on 21 May 2024
Hi Andy,
I tried looking for inbuilt methods for FDTD simulations, but it seems like it is there are none as of now. Like you mentioned, you can extract the points in 3D from scenario.Meshes. I’ve created a script which concatenates points from all buildings:
% Create the UAV scene as you have provided
% Initialize containers for vertices and faces across all meshes
allVertices = [];
allFaces = [];
% Initialize counters for vertices and faces
vertexCount = 0;
% Loop through the building meshes in the scene
for meshIdx = 1:length(scene.Meshes)
% for meshIdx = 1:1
% Extract vertices and faces from the current mesh
vertices = scene.Meshes{meshIdx}.Vertices;
faces = scene.Meshes{meshIdx}.Faces;
% Check if adjustment is needed for the face indices
if vertexCount > 0
% Adjust face indices based on the current vertex count
faces = faces + vertexCount;
end
% Append to the allVertices and allFaces containers
allVertices = [allVertices; vertices];
allFaces = [allFaces; faces];
% Update the vertex count
vertexCount = size(allVertices, 1);
end
The above script also combines the face (triangle connectivity matrix) data with which you can use view the 3D meshes. You can use trimesh function to view the grid:
trimesh(allFaces, allVertices(:,1),allVertices(:,2),allVertices(:,3))
Which will give the following output:
If the points are your only requirement, you can use the allVertices variable.
But I’ve also tried an alternative method to FDTD by modelling the above mesh using createpde, which can be used for electromagnetic analyses, but it requires the mesh to have closed surfaces. Hence, it will require removal of the ground plane (the first mesh) by changing the for loop to:
for meshIdx = 2:length(scene.Meshes)
Then we get the following grid:
Now the vertices and faces can be imported as a mesh, and a model can be created and solved as follows:
% Create a PDE model for electromagnetic wave analysis
model = createpde('electromagnetic', 'electrostatic');
% Create geometry from mesh and append to the PDE model
geometryFromMesh(model, allVertices', allFaces');
% Display the geometry
figure(2);
pdegplot(model);
title('Geometry Imported into PDE Model');
% Define electromagnetic properties for the buildings
electromagneticProperties(model, 'RelativePermittivity', 4, 'RelativePermeability', 1, 'Conductivity', 0.01);
% Apply boundary condition
bc = electromagneticBC(model, 'Face', 1:model.Geometry.NumFaces, 'Voltage', 0);
model.VacuumPermittivity = 8.8541878128e-12;
% Define ChargeDensity source
sc1 = electromagneticSource(model,"Cell",1,"ChargeDensity",0.3)
% Generate mesh
generateMesh(model); % Adjust 'Hmax' based on the desired resolution
% Solve the PDE
result = solve(model);
% Plot the solution
figure(3);
pdeplot3D(model, 'ColorMapData', abs(result.NodalSolution));
title('Field Distribution');
The above is just a sample electromagnetic scenario which you can modify. To get results across time-domain, you can modify and solve the model for each timestep. This could get computationally expensive, but the solution can be parallelized using parfor (https://www.mathworks.com/help/parallel-computing/parfor.html), as solutions for each timestep are independent of each other.
Here are documentation links for some of the functions I used:
I hope this helps.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!