
% Nettoyage
clear all
close all
clc

% Gypsilab path
run('../addpathGypsilab.m')

% Physical parameters
L  = [3 2];        % Room size
X0 = [0. 0.];    % Initial position
tf = 2;            % Final time
c  = 1;            % Sound celerity
N  = 10000;      % number of nodes

% 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

% Create mesh
% mesh = mshSquare(N,L)
mesh = mshDisk(N,1)


% Time discretization
dx = mesh.stp;          % Space step
dx = dx(2);
dt = 0.25*dx/c;      % Time step (CFL)

% Create quadrature
omega = dom(mesh,3)

% Create finite element
phi = fem(mesh,'P1')

% Mesh
figure
plot(mesh)
hold on
plot(omega)
plot(phi)

% Initial boundary condition
X  = mesh.vtx(:,1);
Y  = mesh.vtx(:,2);
U0 = ui(X,Y);
U1 = U0 + dt*vi(X,Y);

% Solution
figure
plot(mesh,U0)

% Finite element matrix
M = integral(omega,phi,phi);
K = integral(omega,grad(phi),grad(phi));

% boucle en temps
t = 0;
while t<tf
    % Schema
    U2 = M\(2*M*U1 - M*U0 - (c*dt)^2*K*U1);
    
    % Animation
    plot(mesh,U2)
    pause(0.001)
    
    % Incrementation
    tmp = U1;
    U1  = U2;
    U0  = tmp;
    t   = t + dt;
end

% % Solution
% figure
% plot(X,U2)





% % Physical parameters
% L  = 1%[3 2];        % Room size
% X0 = 0.7%[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)