% <FEM for ADVECTION EQUATION>
% governing equation: du/dt + c du/dx = 0
% solution: u(t,x)
% boundary condition: u(t,0) = u(t,1)

% ---------------------------------------------------------------------- %

% initial condition: u(0,x)=?
u = [zeros(100,1); ones(300,1); zeros(600,1)];
% u = sin(2*pi*(0:1/100:(1-1/100)))';

% advection speed: c=?
c = 1.0;

% time increment: dt=?
dt = 0.0001;

% ---------------------------------------------------------------------- %

% fem: initialization
n = length(u);
t = 0:dt:(2*1/c);
T = length(t);
dx= 1/n;
x = 0:dx:(n-1)*dx;

% fem: geometry matrix
M = zeros(n,n);
for ii=1:n
    if ii==1
        M(ii,n)   = dx*1/6;
        M(ii,ii)  = dx*2/3;
        M(ii,ii+1)= dx*1/6;
    elseif ii==n
        M(ii,ii-1)= dx*1/6;
        M(ii,ii)  = dx*2/3;
        M(ii,1)   = dx*1/6;
    else
        M(ii,ii-1)= dx*1/6;
        M(ii,ii)  = dx*2/3;
        M(ii,ii+1)= dx*1/6;
    end
end

% fem: convection matrix
C = zeros(n,n);
for ii=1:n
    if ii==1
        C(ii,n)   = -c/2;
        C(ii,ii+1)= c/2;
    elseif ii==n
        C(ii,ii-1)= -c/2;
        C(ii,1)   = c/2;
    else
        C(ii,ii-1)= -c/2;
        C(ii,ii+1)= c/2;
    end
end

% fem: forward differential
% A1 = eye(n)-dt*M\C;
% u1 = [u zeros(n,T-1)];
% for tt=1:T-1
%     u1(:,tt+1)=A1*u1(:,tt);
% end

% fem: backward differential
A2 = (M+dt*C)\M;
u2 = [u zeros(n,T-1)];
for tt=1:T-1
    u2(:,tt+1)=A2*u2(:,tt);
end

% ---------------------------------------------------------------------- %

% open file
[y,m,d]   = ymd(datetime);
ymdstr    = [num2str(y) num2str(m,'%0.2d') num2str(d,'%0.2d')];
filename  = ['ex_01_' ymdstr '.avi'];
writerObj = VideoWriter(filename,'Motion JPEG AVI');
writerObj.FrameRate = 20; %-- flame rate (slide/sec)
open(writerObj); %-- initialization of movie object

% avi visualization
figure(1)
Fnum=100;
for tt=1:round(T/Fnum):T
    set(gcf, 'position', [100   100   480   480 ])
    set(gcf, 'color',[1 1 1])
    plot((0:dx:1)',[u2(:,tt);u2(1,tt)],'k-')
    set(gca, 'position', [0.10  0.08  0.87  0.87])
    set(gca, 'fontname', 'Arial', 'fontsize', 12)
    xlim([0 1])
    ylim([-0.1 1.1])
    set(gca, 'YTick', -1:0.25:1)
    frame = getframe(figure(1)); %-- get frame
    writeVideo(writerObj,frame); %-- add frame to video obj
end

close(writerObj);













