function FilmEvolution_1NucleationSite()
%FilmEvolution_1NucleationSite
%   This function runs a 2D finite element simulation of thin film growth
%   by vapor deposition. The simulation assumes that the film grows as a
%   product of condenstation on a substrate. Film growth is initiated at
%   one user defined nucleation site. The boundary conditions set on the 
%   system, are that the end nodes can only travel in the x direction, thus
%   confining them to the surface of the substrate. Upon completion the 
%   simulation can output an AVI file of the film growth, as well as plots
%   of the film's surface area and volume as a functin of time. 
%
%   ME574 final project
%   Written by: Michael Aldridge and Ibrahim Oraiqat
%   April 2009


% clear workspace
clear; close all; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                   User Defined Simulation Parameters                    %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%  Input Parameters  %%%%%%%%

% define starting configuration
C     = 0;      % particle1 center
I     = 2;      % particle1 Intensity
Xstep = 0.25;   % particle1 node spacing
FWHM  = 0.25;   % particle1 FWHM



% Set Simulation parameters
gamma     = 0.2;      % set coeficient of isotropic surface energy
dgamma    = 0;
g1        = 1;        % set coeficients of volume energy
g2        = 0.2;
m         = 1;        % set element mass
timesteps = 2000;     % set number of timesteps to calculate
dt        = 0.001;    % set timestep length

% Set node cleaning parameters
MinLength   = 0.2; % set minimum element length
MaxLength   = 0.8; % set maximum element length
CleanCycles = 2;   % set number of cleaning cycles to apply
InterpStyle = 2;   % Set 1 for Linear interpolation 2 for cubic Spline


%%%%%%%%  Output Parameters  %%%%%%%%

% Toggle outputs on/off
PlotVol    = 1;         % 0 for no volume plot 1 to turn on plot
PlotArea   = 1;         % 0 for no area plot 1 to turn on plot
recordAVI  = 1;         % 0 for no recording 1 to set recording

% Animation Settings
filename    = '1SiteEvolution'; % filename to save as
fps         = 40;               % set animation frames per second
compression = 'Cinepak';        % set 'none' on OsX, Linux & Unix, 'Indeo3' for WinXP, 'Cinepak' for Vista 
AVIquality  = 70;               % set animation quality

%Animation Window Dimentions
xmin = -4;               % set Animation xmin
xmax =  4;              % set Animation xmax
ymin = -2;              % set Animation ymin
ymax = 10;              % set Animation ymax

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                        Primary Simulation Code                          %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%  Initialize Simulation Conditions  %%%%%%%%

%create initial nucleation sites
x = C-1:Xstep:C+1;
y = I*exp(-(x-C).^2/(2*FWHM^2));


%%%%%%%%  Begin Structural Evolution Loop  %%%%%%%%

for t = 1:timesteps;
 
        % store current frame for animation
        Frames(t) = grabplot(x,y,xmin,xmax,ymin,ymax);

        %determine number of elements
        n = length(x)-1;

        %evolve nodes to next timestep
        [ x, y] = evolve(x,y,n,gamma,dgamma,g1,g2,m,dt);

        % Cleaning removes nodes that are too close together and inserts nodes
        % where they are too sparse. This allows for a larger dt value
        for i = 1:CleanCycles
            [x, y] = CleanNodes(x, y, MinLength, MaxLength, InterpStyle);

        end

        %calculate Surface areas and Volume
        if PlotVol == 1
        vol(t) = Volume(x,y);
        end
        if PlotArea == 1
        SArea(t) = SurfaceArea(x,y);
        end
        if PlotVol == 1 || PlotArea == 1
        time(t) = t*dt;
        end
        
end % end of the time evolution loop

% store final frame for animation
Frames(t) = grabplot(x,y,xmin,xmax,ymin,ymax);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                     Generate Simulation Outputs                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Record Animation of the simulation
if recordAVI == 1
    disp('Writing AVI File')
    movie2avi(Frames,filename,'fps',fps,'compression',compression,'quality',AVIquality);
    disp('Animation Recorded as : ')
    disp( filename )
else
    disp('No Animation Recorded')
end

% Plot Volume vs Time
if PlotVol == 1
figure(2)
plot(time,vol);
xlabel('Elapsed Time');
ylabel('Film Volume');
disp('Volume vs Time plotted in Figure 2') 
end

% Plot Surface Area vs Time
if PlotArea == 1
figure(3)
plot(time,SArea);
xlabel('Elapsed Time');
ylabel('Film Surface Area');
disp('Surface Area vs Time plotted in Figure 3')
end

disp('Simulation Completed Sucessfully')

end % end of simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                  Function to Save Animation Frames                      %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Frame] = grabplot(X1val, Y1val, xmin, xmax, ymin, ymax)
% Function grabplot saves frame of current state for animation
x0 = xmin:xmax;
y0 = zeros(1,length(x0));
figure(1);
plot(X1val,Y1val,'ko-',x0,y0);
axis([xmin, xmax, ymin, ymax]);
Frame = getframe;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                Function to Calculate Time Evolution                     %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ x y ] = evolve( x, y, n, gamma, dgamma, g1, g2, m, dt)
%Function evolve creates H&F matricies and updates x and y values to next
%timestep in structural evolution


% Create the global H and F matrices

%initialize H, F, theta, & L
H = zeros(2*n+2,2*n+2);
Hl = zeros(4);
f = zeros(2*n+2,1);
fl = zeros(4,1);
theta = zeros(n,1);
l = zeros(n,1);

%define elements
%calculate lenght and angle of nodes
for i = 1:n
    %calculate length
    l(i) = ((x(i+1)-x(i))^2+(y(i+1)-y(i))^2)^(1/2);
    %calculate angle
    theta(i) = asin((y(i+1)-y(i))/l(i));
    if x(i+1) < x(i)
        theta(i)=pi-theta(i);
    end
end

for j=1:n;

    %Calculate local H matrix elements
    Hl(1,1)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(1,2)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(1,3)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(1,4)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,1)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(2,2)=l(j)/(6*m)*2*(cos(theta(j)))^2;
    Hl(2,3)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,4)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(3,1)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(3,2)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(3,3)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(3,4)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,1)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(4,2)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(4,3)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,4)=l(j)/(6*m)*2*(cos(theta(j)))^2;

    %add local h matrix to global h matrix
    H(2*j-1:2*j+2,2*j-1:2*j+2) = H(2*j-1:2*j+2,2*j-1:2*j+2) + Hl;


    % calculate Force Vector
    fl(1)=gamma*cos(theta(j))+dgamma*(-sin(theta(j)))+l(j)*(g1-g2)/2*(-sin(theta(j)));
    fl(2)=gamma*sin(theta(j))+dgamma*(cos(theta(j)))+l(j)*(g1-g2)/2*(cos(theta(j)));
    fl(3)=gamma*-cos(theta(j))+dgamma*(sin(theta(j)))+l(j)*(g1-g2)/2*(-sin(theta(j)));
    fl(4)=gamma*-sin(theta(j))+dgamma*(-cos(theta(j)))+l(j)*(g1-g2)/2*(cos(theta(j)));

    %add local force vector to global force vector
    f(2*j-1:2*j+2)=f(2*j-1:2*j+2)+fl;


end

%Add Artificial Finite Mobility (small number to diagonal points in H matrix)
for m=1:(2*n+2);
    H(m,m)=0.000001+H(m,m);
end

%boundary conditions
%fix end points comment and uncomment lines to turn off/on boundary
%conditions

%H(end-1,end-1) = (1+H(end-1,end-1))*1E10;  %fix last point in x direction
H(end,end) = (1+H(2*n+2,2*n+2))*1E10;      %fix last point in y direction
%H(1,1) = (1+H(1,1))*1E10;                  %fix first point in x direction
H(2,2) = (1+H(2,2))*1E10;                  %fix first point in y direction


%Sometimes this error will print. Use this command to turn it off
warning off MATLAB:nearlySingularMatrix;

%solve for velocities
vel=H\f;

%update positions using Euler method
for i=1:(n+1);
    x(i)=x(i)+vel(2*i-1)*dt;
    y(i)=y(i)+vel(2*i)*dt;
end
end % end of evolve function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                    Function to Clean Nodal Points                       %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to clean up nodes that are either too short or too long
function [x y] = CleanNodes(x,y,MinLength,MaxLength,InterpStyle)

% Find number of nodal pooints
n = length(x);

% Find elements that are too small
flag = 0;
for i = 1:n-2
    if ((x(i+1-flag)-x(i-flag))^2+(y(i+1-flag)-y(i-flag))^2)^(1/2) < MinLength

        % Delete extra nodal point
        x(i+1-flag) = [];
        y(i+1-flag) = [];

        % fix indices
        i = i+1;
        flag = flag +1;
    end
end

% check check to see if last element is too small
if ((x(end)-x(end-1))^2+(y(end)-y(end-1))^2)^(1/2) < MinLength
    %delete extra nodal point
    x(end-1) = [];
    y(end-1) = [];
end

% Find new number of nodal points
n = length(x);

% find elements that are too big
for i = 2:n
    if ((x(i)-x(i-1))^2+(y(i)-y(i-1))^2)^(1/2) > MaxLength

        
        % Interpolate new nodal point
        % pick a new x location
        newx = (x(i)+x(i-1))/2;
        
        if InterpStyle == 1
        % Calculate linear interpolation for y value
        newy = (y(i)+y(i-1))/2;
        else
            if  InterpStyle == 2 
        % Calculate cubic spline interpolation for y value
        newy = spline(x,y,newx);
            end
        end


        % Build new list of nodal points
        xtemp = zeros(1,length(x)+1);
        xtemp(1,1:i-1)= x(1,1:i-1);
        xtemp(i) = newx;
        xtemp(1,i+1:end) = x(1,i:end);

        ytemp = zeros(1,length(y)+1);
        ytemp(1,1:i-1)= y(1,1:i-1);
        ytemp(i) = newy;
        ytemp(1,i+1:end) = y(1,i:end);

        % Update nodal points
        x = xtemp;
        y = ytemp;

        % Adjust loop index
        i=i-1;
    end
end

end %end of CleanNodes function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%          Functions to Calculate Film Surface Area and Volume            %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate surface area
function [sa] = SurfaceArea(x,y)
sa = 0;
for i = 1: length(x) -1
    sa = sa + ((x(i+1)-x(i))^2+(y(i+1)-y(i))^2)^(1/2);
end
end

%calculate volume
function [v] = Volume(x,y)
v = 0;
for i = 1: length(x)-1
    v = v + (x(i+1)-x(i))*(y(i+1)+y(i))/2;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
