clc;
clear all; 
close all
global N_step
n_iteration=500; %no.of iterations
%declaring constants
deltaT=0.0030;%time step
g0=3;
gI=2;
gII=.5;
m=1.5; %mobility
T = deltaT;
min_length=0.5;
max_length=3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% MAIN CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial nodes, elements;
[node,elem]=mesh_generator;


% start of iteration
for i=1:n_iteration
%     N_step = v;
    % update
    n_node=size(node,1);
    n_elem=size(elem,1);
%     n_DOF= 4*n_node;    
    %update LM
    [ID, IEN, LM]=noderelation(node, elem);
    % plot and save configuration
    t=deltaT*i*g0*m/10^2; % time
%     Particle_plot(node,t);
    draw(node,t);
    F(i)=getframe;
        
    %Frame(i)=getframe;
    filename=['Evolution',num2str(i)];
    hgsave(eval(['filename']));
    % compute global Viscosity Matrix H
    clear H; % for n_node change
    Hglobal = H_global(ID, LM, node, elem, m);
   
    clear f; % for n_node change
    fglobal = f_global(elem, node, LM, g0, gI, gII, i);
    % get q
    q = q_assemble(node, ID);
    % compute qdot
    clear qdot;
    qdot=Hglobal\fglobal;
    % update general coordinate position or DOF, q
    q=q+qdot*deltaT;
    % update nodes
    for e=1:n_elem;
        node(e,2)=q(LM(1,e));
        node(e,3)=q(LM(2,e));
    end
    % check element size & add or remove nodes(elements)
    [node, elem] = nodes_check(node,elem, max_length, min_length);
end

% make movie file
movie(F);
movie2avi(F,'meramovie2.avi');
% save filename.mat F
% mpgwrite(F,jet,'movie.mpg');
% mpeg_play movie.mpg

%movie2avi(movie,Evolution,'colormap',236,'compresion', 'Indeo5','fps',15,'quality',75);
