% <FEM for DIV EQUATION>
% governing equation: du/dt + k d2u/dx2 = 0
% solution: u(t,x)
% boundary condition: free

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

% initial condition: u(0,x)=?
u = [zeros(100,1); 200*ones(1,1); zeros(99,1)];

% divergence coef: k=?
k = 0.1*(0.01)^2/2/0.0001;

% time increment: dt=?
dt = 0.0001;

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

% fem: initialization
n = length(u); %-- num of elem
t = 0:dt:1;
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*1/3;
        M(ii,ii+1)= dx*1/6;
    elseif ii==n
        M(ii,ii-1)= dx*1/6;
        M(ii,ii)  = dx*1/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: stifness matrix
K = zeros(n,n);
for ii=1:n
    if ii==1
%         K(ii,n)   = -1/dx;
        K(ii,ii)  = 1/dx;
        K(ii,ii+1)= -1/dx;
    elseif ii==n
        K(ii,ii-1)= -1/dx;
        K(ii,ii)  = 1/dx;
%         K(ii,1)   = -1/dx;
    else
        K(ii,ii-1)= -1/dx;
        K(ii,ii)  = 2/dx;
        K(ii,ii+1)= -1/dx;
    end
end

% fem: forward/backward differential
A1 = eye(n)-k*dt/dx*K;
A2 = (M+k*dt*K)\M;
u1 = [u zeros(n,T-1)];
u2 = [u zeros(n,T-1)];
for tt=1:T-1
    u1(:,tt+1)=A1*u1(:,tt);
    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_04_' 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=200;
for tt=1:round(T/Fnum):T
    set(gcf, 'position', [100   100   480   480 ])
    set(gcf, 'color',[1 1 1])
    plot((0:dx:1)',[u1(:,tt);u1(1,tt)],'b-');hold on
    plot((0:dx:1)',[u2(:,tt);u2(1,tt)],'r-');hold off
    set(gca, 'position', [0.10  0.08  0.87  0.87])
    set(gca, 'fontname', 'Arial', 'fontsize', 12)
    xlim([0 1])
    ylim([0 10])
    frame = getframe(figure(1)); %-- get frame
    writeVideo(writerObj,frame); %-- add frame to video obj
end

close(writerObj);













