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