% fdtdWave2D.m - Matthieu Aussal - Copyright 2018 - Licence CC BY-NC-SA 4.0
% CMAP - Ecole Polytechnique, Route de Saclay, 91128 Palaiseau CEDEX France
% For any informations, please contact  : matthieu.aussal@polytechnique.edu

clear all
close all
clc

% Physical parameters
L  = [3 2];        % Room size
X0 = [2.2 1.2];    % Initial position
tf = 1;            % Final time
c  = 1;            % Sound celerity

% Initial condition for magnitude and speed
ui = @(x,y) exp(-(x-X0(1)).^2/1e-2) .* exp(-(y-X0(2)).^2/1e-2); % gaussian
vi = @(x,y) zeros(size(x));                                     % null

% Numerical discretization
dx = 0.01;          % Space step (both for x and y)
dt = 0.5*dx/c;      % Time step (CFL)

% Regular grid (x moving first)
x       = 0:dx:L(1);
y       = 0:dx:L(2);
[xx,yy] = ndgrid(x,y);
X       = [xx(:) yy(:)];

% Dimensions
Nx = length(x);
Ny = length(y);
N  = Nx * Ny;

% 2D Laplacian operator 
e  = ones(N,1);
Mx = spdiags(e*[1 1 -4 1 1]/(dx^2),[-Nx,-1,0,1,Nx],N,N);

% Delete values non connected vertices (end of x line)
I       = Nx:Nx:N-1;  % row indices
J       = Nx+1:Nx:N;  % column indices 
ind     = sub2ind([N,N],[I;J],[J;I]); % convert to vector indices
Mx(ind) = 0;
% figure
% spy(Mx)

% Vertex multiplicity (number of neighbours)
mlt = full(sum(abs(Mx)>0,2)) - 1;

% Add Neumman condition using multiplicity : 4 -> multiplicity
Mx(sub2ind([N,N],1:N,1:N)) = -mlt/(dx^2);  

% Initial magnitude and speed (u0 and v0)
u0 = ui(X(:,1),X(:,2));
v0 = vi(X(:,1),X(:,2));

% Force initial condition
u0(mlt<4) = 0;

% First step scheme (Forward Euler) 
u1 = u0 + dt * v0;

% Graphical representation
figure(10)
uu   = reshape(u1,Nx,Ny);
Hfig = surf(xx,yy,uu);
axis([0 L(1) 0 L(2) 0 1])
axis equal
caxis([0 1])
shading interp
view(-10,60)
camlight
material dull

% Video output
vidObj = VideoWriter('2Dneumann.avi');
open(vidObj);

% Order two sheme in times
t = dt;
while t < tf
    % Leap-frog sheme
    u2 = 2*u1 - u0  + (c*dt)^2 * (Mx * u1);
    t  = t + dt;
    
    % Graphical representation
    uu = reshape(u2,Nx,Ny);
    set(Hfig,'ZData',uu,'CData',uu)
    title(['Wave propagation t = ',num2str(t,'%6.3f')])
    
    % Write each frame to the file.
    currFrame = getframe(gcf);
    writeVideo(vidObj,currFrame);
    pause(dt)
    
    % Incrementation
    u0 = u1;
    u1 = u2;
end

 % Close the file.
 close(vidObj);





