geomap.m

非连续波导时间域有限差分法仿真（TD-FDM）

% add some comments here
% describe clearly why we need the matrix 'index'

% the disadvantange of this method here is to 
% waste storage on unvalid nodes of our geometry.
% if the geometry is large and complicated,
% it's not appropriate. :(
% Advantange is just for simplicity of programming.

nmbr = zeros(MAX_Y,MAX_X);         % for node numbering
type = zeros(MAX_Y,MAX_X);        % for node type (inner/boundary/invalid)

% type each node with increament number
cnt=1;
for j=1:Nx
    for i=1:Ny
        nmbr(i,j) = cnt;
        cnt = cnt+1;
    end
end
for j=Nx+1:Nx+Nix
    for i=Nh:Ny-(Nh-1)
        nmbr(i,j) = cnt;
        cnt = cnt+1;
    end
end
for j=Nx+Nix+1:MAX_X
    for i=1:Ny
        nmbr(i,j) = cnt;
        cnt = cnt+1;
    end
end

% invalid nodes, exclude them!!!
for i=1:Nh-1
    for j=Nx+1:Nx+Nix
        type(i,j) = inf;
        type(Ny-i+1,j) = inf;
    end
end

% boundary edge classification
for j=2:Nx-1
    type(1,j) = 1;                 % y=0, bottom boundary
    type(1,j+Nx+Nix) = 1;
    type(Ny,j) = 2;                % y=H, top boundary
    type(Ny,j+Nx+Nix) = 2;
end

for j=Nx+1:Nx+Nix
    type(Nh,j) = 1;                % y=h, bottom of discontinuity
    type(Ny-(Nh-1),j) = 2;         % y=H-h, top of discontinuity
end

for i=2:Ny-1
    type(i,1) = 3;                 % x=0
    type(i,MAX_X) = 6;             % x=2L+l
end

for i=2:Nh-1
    type(i,Nx) = 4;                % x=L
    type(i+Nh+Niy-2,Nx) = 4;
    type(i,Nx+1+Nix) = 5;          % x=L+l
    type(i+Nh+Niy-2,Nx+1+Nix) = 5;
end

% corner nodes%%%%%%%%%%%%%%%%%%%%%%
type(1,1) = 7;
type(1,Nx+1+Nix) = 7;
type(Ny,1) = 8;
type(Ny,Nx+1+Nix) = 8;
type(Ny,Nx) = 9;
type(Ny,MAX_X) = 9;
type(1,Nx) = 10;
type(1,MAX_X) = 10;


%%%%%%%%%%%%%%%%%%%%%%????????????????????????%%%%%%%%%%%%%%%%%%%%%%%%%%
%for j=1:MAX_X
%    type(:,j,3)=(j-1)*hx-(2*L+l)/2;
%end
%for i=1:Ny
%    type(i,:,4)=H/2-(i-1)*hy;
%end