%ThinFilm is a numerical simulation of a grain boundary in a thin film.
%Chris Welch 3/30/10

clear all;

%Surface boundary
x=[0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40];
y=[3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3];
%Grain boundary
xb = [10 10 10 10 10 10 10 10 10 10];
yb=[2.75 2.5 2.25 2 1.75 1.5 1 .75 .25 0];
%Define the point in the set that is where the grain boundary meets the
%surface boundary
TriPt = 6;

NumberOfIterations = 3000;
lowerElementLength = 0.4;
upperElementLength = 1.5;

deltaT = .0000002;
L = 1;
M = .1;
gammaS=.75;
gammaB=1;
gammaI=.999;
dgammaS=0;
dgammaB=0;
g1=1;g2=1;


for v = 1:NumberOfIterations
v;
close;
figure;
%xb and xy do not contain connecting point.  This command includes it in
%the vector for plotting.
xbp = [x(TriPt), xb(1,:)];
ybp = [y(TriPt), yb(1,:)];

plot(x,y,'ko-',xbp,ybp,'ko-');
axis([0 35 0 8]);

%Record the plot
F(v)=getframe;


NumSurfE = length(x)-1;
NumBoundE = length(xb);
%Create the empty H matrix and f vector
H = zeros((length(x)+length(xb))*3+NumSurfE+NumBoundE);
f = zeros((length(x)+length(xb))*3+NumSurfE+NumBoundE,1);

%add a small mass to the nodes to avoid zero/zero division
%(eye is a command to create an identity matrix)
H = H + .000001*eye(size(H));

%*************************************************************
%Create the H and f matrix for the surface boundary
%*************************************************************

%Create boundary matricies and add to global matrices
for e = 1:(length(x)-1)
    
    %Calculate theta and l for each element in the surface
    theta = atan((y(e+1)-y(e))/(x(e+1)-x(e)));
    if x(e+1)<x(e)
        theta = pi + theta;
    end
    l = ((y(e+1)-y(e))^2+(x(e+1)-x(e))^2)^(1/2);
    
    %LocalHmatrix calculates the 4x4 local H matrix based on theta, l and L
    LocalH = LocalHmatrixDiffusion(theta,l,L,M);
   
    %Insert local H matrix into global matrix
    H(4*e-3:4*e+3,4*e-3:4*e+3) = ...
        H(4*e-3:4*e+3,4*e-3:4*e+3) + LocalH(:,:);
    %for every iteration of this loop, the next local matrix gets added
    %into the pre-existing global matrix, allowing for summation
    
    %Create local f vector for individual element
    Localf = gammaS*[cos(theta);sin(theta);-cos(theta);-sin(theta)]...
        + dgammaS*[-sin(theta);cos(theta);sin(theta);-cos(theta)]...
        + l*(g1 - g2)/2*[-sin(theta);cos(theta);-sin(theta);cos(theta)];
    
    %Insert the local vector into the large global vector
    f(4*e-3:4*e-2,1) = f(4*e-3:4*e-2,1)+Localf(1:2,1);
    f(4*e+1:4*e+2,1) = f(4*e+1:4*e+2,1)+Localf(3:4,1);
    
end

%Create boundary conditions by adding huge masses to endpoints  
M = 1000000;
H(1,1) = H(1,1) + M;
H(2,2) = H(2,2) + M;
%Un-Comment this next line to fix the Triple point in the vert. direction
%H(TriPt*2,TriPt*2) = H(TriPt*2,TriPt*2) + M;
H(4*length(x)-3,4*length(x)-3) = H(4*length(x)-3,4*length(x)-3) + M;
H(4*length(x)-2,4*length(x)-2) = H(4*length(x)-2,4*length(x)-2) + M; 

%*************************************************************
%Create the H and f matrix for the grain boundary
%*************************************************************

%Create grain boundary matricies and add to global matrices.
%This is exactly as the loop above, except it reads the boundary vector of
%x and y positions and adds the local matricies to the end of the global
%matricies that were started above.
%Note: this does not include the element connecting the grain boundary to the
%surface
for e = 1:(length(xb)-1)
    theta = atan((yb(e+1)-yb(e))/(xb(e+1)-xb(e)));
    if xb(e+1)<xb(e)
        theta = pi + theta;
    end
    l = ((yb(e+1)-yb(e))^2+(xb(e+1)-xb(e))^2)^(1/2);
    
    %LocalHmatrix calculates the 4x4 local H matrix based on theta, l and L
    LocalH = LocalHmatrixDiffusion(theta,l,L,M);
  
    H(4*(length(x)+e)-3:4*(length(x)+e)+3,4*(length(x)+e)-3:4*(length(x)+e)+3) =...
        H(4*(length(x)+e)-3:4*(length(x)+e)+3,4*(length(x)+e)-3:4*(length(x)+e)+3)...
        + LocalH(:,:);
    
    Localf = gammaB*[cos(theta);sin(theta);-cos(theta);-sin(theta)]...
        + dgammaB*[-sin(theta);cos(theta);sin(theta);-cos(theta)]...
        + l*(g2 - g2)/2*[-sin(theta);cos(theta);-sin(theta);cos(theta)];
    
    f(4*(length(x)+e)-3:4*(length(x)+e)-2,1) = f(4*(length(x)+e)-3:4*(length(x)+e)-2,1)+Localf(1:2,1);
    f(4*(length(x)+e)+1:4*(length(x)+e)+2,1) = f(4*(length(x)+e)+1:4*(length(x)+e)+2,1)+Localf(3:4,1);
end

%Fix the bottom point to the ground by adding a huge mass to y
H(4*(length(x)+length(xb))-2,4*(length(x)+length(xb))-2) = ...
    H(4*(length(x)+length(xb))-2,4*(length(x)+length(xb))-2) + M;

%Calculate the element connecting the surface to the grain boundary
theta = atan((yb(1)-y(TriPt))/(xb(1)-x(TriPt)));
if xb(1)<x(TriPt)
    theta = pi + theta;
end
l = ((yb(1)-y(TriPt))^2+(xb(1)-x(TriPt))^2)^(1/2);

%Calculate the local matrix and f vector for the connecting element
LocalH = LocalHmatrixDiffusion(theta,l,L,M);

Localf = gammaB*[cos(theta);sin(theta);-cos(theta);-sin(theta)]...
        + dgammaB*[-sin(theta);cos(theta);sin(theta);-cos(theta)]...
        + l*(g2 - g2)/2*[-sin(theta);cos(theta);-sin(theta);cos(theta)];

%Four corners of locations in global matrix for local matrix
lower = 4*TriPt-3;
upper = 4*(length(x)+1)-4;

H(lower:(lower+2),lower:(lower+2)) = ...
    H(lower:(lower+2),lower:(lower+2)) + LocalH(1:3,1:3);
H(lower:(lower+2),upper:(upper+3)) = ...
    H(lower:(lower+2),upper:(upper+3)) + LocalH(1:3,4:7);
H(upper:(upper+3),lower:(lower+2)) = ...
    H(upper:(upper+3),lower:(lower+2)) + LocalH(4:7,1:3);
H(upper:(upper+3),upper:(upper+3)) = ...
    H(upper:(upper+3),upper:(upper+3)) + LocalH(4:7,4:7); 

f(lower:(lower+1),1) = f(lower:(lower+1),1) + Localf(1:2,1);
f(upper+1:(upper+2),1) = f(upper+1:(upper+2),1) + Localf(3:4,1);

%Add gamma from the floor
f(length(f)-2) = f(length(f)-2) + gammaI;

%solve for velocity vector
vel = H\f;

%update positions
for i = 1:length(x)
    x(i) = x(i) + deltaT*vel(4*i-3);
    y(i) = y(i) + deltaT*vel(4*i-2);
end

for i = 1:length(xb)
    xb(i) = xb(i) + deltaT*vel(4*(length(x)+i)-3);
    yb(i) = yb(i) + deltaT*vel(4*(length(x)+i)-2);
end


%


%**************************************************************
%CHECK FOR ELEMENTS TOO LARGE OR SMALL
%**************************************************************

%******Checking in surface boundary*********
 
%too small: condense two nodes into the bisecting point on their
%element
for e = 1:(length(x)-1)
    l = ((y(e+1)-y(e))^2+(x(e+1)-x(e))^2)^(1/2);
    
    if l < lowerElementLength && e < (length(x)-1)
        x(e) = (x(e)+x(e+1))/2;
        y(e) = (y(e)+y(e+1))/2;
        x(e+1:length(x)-1) = x(e+2:length(x));
        y(e+1:length(y)-1) = y(e+2:length(y));
        x = x(1:length(x)-1);
        y = y(1:length(y)-1);
        %check if triple point changes
        if e < TriPt
            TriPt = TriPt - 1;
        end
        
    end
    %when an element is removed, the loop will try to access an element
    %that doesn't exist at the end, so break as soon as the end is reached
    if e == length(x)-1
        break
    end
end

%End element too small: delete second to last node
sln = length(x) - 1;
l = ((y(sln+1)-y(sln))^2+(x(sln+1)-x(sln))^2)^(1/2);
if l < lowerElementLength
    x(sln) = x(sln + 1);
    y(sln) = y(sln + 1);
    x = x(1:length(x)-1);
    y = y(1:length(y)-1);
    if TriPt == sln
        TriPt = length(x);
    end
end


%too large: add a node at the bisecting point of the element

%for every element added, the e index needs to be increased in order to hit
%all points in the position vectors
add = 0;
for e = 1:(length(x)-1)
    ee = e + add;
    l = ((y(ee+1)-y(ee))^2+(x(ee+1)-x(ee))^2)^(1/2);
    
    
    if l > upperElementLength
        newx = (x(ee)+x(ee+1))/2;
        newy = (y(ee)+y(ee+1))/2;
        x(ee+2:length(x)+1) = x(ee+1:length(x));
        x(ee+1) = newx;
        y(ee+2:length(y)+1) = y(ee+1:length(y));
        y(ee+1) = newy;
        add = add + 1;
        
        if ee < TriPt
            TriPt = TriPt + 1;
        end
    end
end


%********Checking in grain boundary**********

xb = [x(TriPt), xb(1,:)];
yb = [y(TriPt), yb(1,:)];
%too small: condense two nodes into the bisecting point on their
%element
for e = 1:(length(xb)-1)
    l = ((yb(e+1)-yb(e))^2+(xb(e+1)-xb(e))^2)^(1/2);
    
    if l < lowerElementLength && e < (length(xb)-1)
        xb(e) = (xb(e)+xb(e+1))/2;
        yb(e) = (yb(e)+yb(e+1))/2;
        xb(e+1:length(xb)-1) = xb(e+2:length(xb));
        yb(e+1:length(yb)-1) = yb(e+2:length(yb));
        xb = xb(1:length(xb)-1);
        yb = yb(1:length(yb)-1);
        
    end
    %when an element is removed, the loop will try to access an element
    %that doesn't exist at the end, so break as soon as the end is reached
    if e == length(xb)-1
        break
    end
end

%End element too small: delete second to last node
sln = length(xb) - 1;
l = ((yb(sln+1)-yb(sln))^2+(xb(sln+1)-xb(sln))^2)^(1/2);
if l < lowerElementLength
    xb(sln) = xb(sln + 1);
    yb(sln) = yb(sln + 1);
    xb = xb(1:length(xb)-1);
    yb = yb(1:length(yb)-1);
end

%too large: add a node at the bisecting point of the element

%for every element added, the e index needs to be increased in order to hit
%all points in the position vectors
add = 0;
for e = 1:(length(xb)-1)
    ee = e + add;
    l = ((yb(ee+1)-yb(ee))^2+(xb(ee+1)-xb(ee))^2)^(1/2);
    

    if l > upperElementLength
        newx = (xb(ee)+xb(ee+1))/2;
        newy = (yb(ee)+yb(ee+1))/2;
        xb(ee+2:length(xb)+1) = xb(ee+1:length(xb));
        xb(ee+1) = newx;
        yb(ee+2:length(yb)+1) = yb(ee+1:length(yb));
        yb(ee+1) = newy;
        add = add + 1;
    end
end

xb = xb(2:length(xb));
yb = yb(2:length(yb));

%}

end

movie2avi(F,'test34','fps',60,'compression','Cinepak','quality',15);
close;