
% Nettoyage
clear all
close all
clc

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

% Physical parameters
L  = 1;        % String size
x0 = 0.0;      % Initial position
tf = 1;        % Final time
c  = 1;        % Sound celerity
N  = 100;      % number of nodes

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

% Create mesh
mesh = mshSegment(N,L)

% Time discretization
dx = mesh.stp;          % Space step
dx = dx(2)
dt = 0.5*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);
U0 = ui(X);
U1 = U0 + dt*vi(X);

% Solution
figure
h = plot(X,U0)
grid on
axis([-L/2 L/2 -1 1])

% 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
    set(h,'YData',U2)
    pause(0.1)
    
    % 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)