%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ME599 PROJECT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%THIS PROGRAM CALCULATES THE VELOCITIES OF THE ELEMENTAL NODES AND
%%%%%CALCULATES THE NEXT POSITION USING EULER'S METHOD. THIS CODE DOESNOT 
%%%%%INCORPORATES THE STRESS EFFECTS ON THE NODES 

function [] = microproject()

time = 30;    %%%Total time for which the program is to be run
N = 400;      %%%Number of time steps
h = time/N;   %%%Length of each time step

nelements = 30;   %%%%Total number of elements

nperelement = 2;  %%%%Number of nodes per element, in this case, is two for
% % % % % % % % % % % each element

nnodes = (nperelement-1)*nelements; %%Total number of nodes are dependent
% % % % % % % % % % % % % % % % % % % % % % % % % on the number of elements

%%%%%%%%%%%%%%%%%%%%Execution of the Euler's Method to find the next
%%%%%%%%%%%%%%%%%%%%position of the nodes%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t(1) = 0;
for n = 1:N
    for A = 1:nnodes
        NNODES(A) = A;
    end
    t(n+1) = t(n) + h;   %%%%%Time increment

    load circle.dat;     %%%%%Circle.dat contains the initial positions of
    % % % % % % % % % % % % %the nodes% or the initial shape of the microcrack.

    QDOT = me599fem(nelements);  %Nodal velocities

    for A = 1:nnodes
        X(A) = circle(A,2) + h*QDOT(A,1);   %%%%Next x-position
        Y(A) = circle(A,3) + h*QDOT(A,2);   %%%%Next x-position
    end

    J = X;
    K = Y;
    J(nnodes+1) = X(1)
    K(nnodes+1) = Y(1)

    %%J and K are required to plot the shape as a closed curve

    %%%%Plot the shape%%%%%%
    plot(J,K);
    axis([-0.75, 0.75, -0.75, 0.75])
    axis off
    %%%%%%%%%%%%%%%%%%%%%%%

    Q = [NNODES;X;Y]
    fid = fopen('circle.dat','wt');
    fprintf(fid, '%8.0f %20.18f %20.18f\n',Q);
    fclose(fid)

    %%%The fprint command outputs the next/new nodal positions in the
    %%%format which are required by the ANSYS to calculate the
    %%%stress/strain values to calculate the strain energy densities

end

function [QDOT] = me599fem(nelements)

%This function calculates the nodal velocities with nelements as the
%input

%%%%%%%%%%%%%%%%%%%%%Step 1: Read nodes/coordinates%%%%%%%%%%%%%%%%%%%%%%%%

nperelement = 2;
nnodes = (nperelement-1)*nelements;

M = 0.00000001;
g0 = 1;
m = 0.01;
ki = 0;

load circle.dat;  %%Load the nodal positions

circle(:,:);

%%Vector X(A) and Y(A) specifies the x and y-coordinates of the nnodes%%%%%

for A = 1:nnodes
    X(A)= circle(A,2);
    Y(A)= circle(A,3);
end

J = X;
K = Y;
J(nnodes+1) = J(1);
K(nnodes+1) = K(1);

%%J and K are required to plot the shape as a closed curve

%%%%%%%%%%Step 2: Read boundary conditions and create ID array%%%%%%%%%%%%%

%%%%%%%%%%ID array gives the equation number related to each node%%%%%%%

p = 1;
for A=1:nnodes
    for i = 1:3
        G(i,A)=NaN;
        if isnan(G(i,A))
            ID(i,A) = p;
            p = p+1;
        else
            ID(i,A) = 0;  %%%%%%i = 1: 3 because each node has three degrees
            % %           of freedom, two related to positions and one related to flux
        end

    end
end

g = max(max(ID));       %%This will define the dimensions of K and F vector

%%%%%%%%%%%%%%%%%Step 4: Read IEN array and LM array%%%%%%%%%%%%%%%%%%%%%%%
%%%%%IEN array relates the local nodes of an element e to the global nodes%

for e=1:nelements
    for a=1:nperelement
        for i = 1:3
            if (e == nelements) && (a == nperelement)
                IEN(a,e) = 1;
                LM(i,a,e)=ID(i,IEN(a,e));
            else
                A=(e-1)*(nperelement-1)+a;
                IEN(a,e)=A;
                LM(i,a,e)=ID(i,IEN(a,e));
            end
        end

    end
end


%%Step 5:Define local and global(assembled) stiffness matrix & force vector

Hsize = nelements+3*nnodes;

%%%Size of global H matrix is nelements + 3*nnodes because each element has
%%%7 degrees of freedom, six associated with each of the two nodes and one
%%%related to flux (Jm) at the mid of the element

%Initialization of the global H matrix

H=zeros(nelements+3*nnodes,nelements+3*nnodes);

%Add matrix is set to ensure that the H matrix is not singular

add=[1 0 0 0 0 0 0; 0 1 0 0 0 0 0; 0 0 0 0 0 0 0; 0 0 0 1 0 0 0; 0 0 0 0 1 0 0; 0 0 0 0 0 0 0; 0 0 0 0 0 0 0];


for e=1:nelements
    
    %Find the node numbers
    
    node1= IEN(1,e);
    node2= IEN(2,e);

    %Nodal positions
    x1=X(node1);
    y1=Y(node1);
    x2=X(node2);
    y2=Y(node2);

    %The length of the element
    l(e)=sqrt((x2-x1)^2+(y2-y1)^2);

    %Angle of each element
    if and(x1<x2,y1==y2) 
        angle(e)=0;  %%%Horizontal element
        
    elseif and(x1<x2,y1<y2) %First quadrent
        angle(e)=asin((y2-y1)/l(e));
        
    elseif and(x1==x2,y1<y2) %+ vertical
        angle(e)=pi/2;
        
    elseif and(x1>x2,y1<y2) %Second quadrent
        angle(e)=pi-asin((y2-y1)/l(e));
        
    elseif and(x1>x2,y1==y2) %- horizontal
        angle(e)=pi;
        
    elseif and(x1>x2,y1>y2) %Third quadrent
        angle(e)=pi+asin((y1-y2)/l(e));
        
    elseif and(x1==x2,y1>y2) %- vertical
        angle(e)=1.5*pi
        
    elseif and(x1<x2,y1>y2) %Fourth quadrent
        angle(e)=2*pi-asin((y1-y2)/l(e));
        
    end

end

% In order to ensure that the two consequtive elements are not inline with 
% each other, a random deviation is assigned to the x and y coordinates 
% random number between 0.99 and 1.01 times the current position

if abs(angle(1)-angle(nelements))<0.01
    X(IEN(1,1)) = X(IEN(1,1))*(0.99+(1.01-0.99)*rand);
    Y(IEN(1,1)) = Y(IEN(1,1))*(0.99+(1.01-0.99)*rand);
end

for e = 2:nelements
    if abs(angle(e)-angle(e-1))<0.01
        X(IEN(1,e)) = X(IEN(1,e))*(0.99+(1.01-0.99)*rand);
        Y(IEN(1,e)) = Y(IEN(1,e))*(0.99+(1.01-0.99)*rand);
    end
end


%%Step 5:Assembly of global viscosity matrix

for e = 1:nelements;

    node1=IEN(1,e);
    node2=IEN(2,e);

    %The nodal posetions
    x1= X(node1);
    y1= Y(node1);
    x2= X(node2);
    y2= Y(node2);

    %Length of the element
    l=sqrt((y2-y1)^2+(x2-x1)^2);

    S=(y2-y1)/l;
    C=(x2-x1)/l;


    u=(m*(l^2))/(5*M);

    %Initialize the element viscosity matrix
    
    He=zeros(7);

    %Determine the elemental viscosity matrix
    
    He=(1/(6*m*l))*[2*l^2*S^2 -2*l^2*S*C  5*l*S  l^2*S^2 -l^2*S*C -l*S -4*l*S; -2*l^2*S*C 2*l^2*C^2 -5*l*C -l^2*S*C  l^2*C^2  l*C 4*l*C; 5*l*S -5*l*S 14+4*u l*S -l*C 2-u -16+2*u; l^2*S^2 -l^2*S*C l*S 2*l^2*S^2 -2*l^2*S*C -5*l*S 4*l*S; -l^2*S*C l^2*C^2 -l*C -2*l^2*S*C 2*l^2*C^2 5*l*C -4*l*C; -l*S l*C 2-u -5*l*C 5*l*C 14+4*u -16+2*u; -4*l*S 4*l*C -16+2*u 4*l*S -4*l*C -16+2*u 32+16*u];

    %Add the add matrix to He to supress singularity
    
    He=He+(1/10^6)*add;


    if e==1 %First element
        H(1:7,1:7)=H(1:7,1:7)+He(:,:);
    
    elseif e==2 %Second element
        H(4:6,4:6)=H(4:6,4:6)+He(1:3,1:3);
        H(4:6,8:11)=H(4:6,8:11)+He(1:3,4:7);
        H(8:11,4:6)=H(8:11,4:6)+He(4:7,1:3);
        H(8:11,8:11)=H(8:11,8:11)+He(4:7,4:7);
    
    elseif e==nelements %Last element
        H(1:3,1:3)=H(1:3,1:3)+He(4:6,4:6);
        H(Hsize-4:Hsize-2,Hsize-4:Hsize-2)=H(Hsize-4:Hsize-2,Hsize-4:Hsize-2)+He(1:3,1:3);
        H(Hsize,Hsize)=H(Hsize,Hsize)+He(7,7);
        H(1:3,Hsize-4:Hsize-2)=H(1:3,Hsize-4:Hsize-2)+He(4:6,1:3);
        H(1:3,Hsize)=H(1:3,Hsize)+He(4:6,7);
        H(Hsize-4:Hsize-2,1:3)=H(Hsize-4:Hsize-2,1:3)+He(1:3,4:6);
        H(Hsize,1:3)=H(Hsize,1:3)+He(7,4:6);
        H(Hsize-4:Hsize-2,Hsize)=H(Hsize-4:Hsize-2,Hsize)+He(1:3,7);
        H(Hsize,Hsize-4:Hsize-2)=H(Hsize,Hsize-4:Hsize-2)+He(7,1:3);
    
    else %Other elements
        H(4*(e-1):4*(e-1)+2,4*(e-1):4*(e-1)+2)=H(4*(e-1):4*(e-1)+2,4*(e-1):4*(e-1)+2)+He(1:3,1:3);
        H(4*(e-1):4*(e-1)+2,4*e:4*e+3)=H(4*(e-1):4*(e-1)+2,4*e:4*e+3)+He(1:3,4:7);
        H(4*e:4*e+3,4*(e-1):4*(e-1)+2)=H(4*e:4*e+3,4*(e-1):4*(e-1)+2)+He(4:7,1:3);
        H(4*e:4*e+3,4*e:4*e+3)=H(4*e:4*e+3,4*e:4*e+3)+He(4:7,4:7);
    
    end

end

%%Step 7:Assembly of global force matrix

F=zeros(nelements + 3*nnodes,1);

for e = 1:nelements;

    node1=IEN(1,e);
    node2=IEN(2,e);

    %The nodal positions
    
    x1= X(node1);
    y1= Y(node1);
    x2= X(node2);
    y2= Y(node2);

    %Length of the element
    
    l=sqrt((y2-y1)^2+(x2-x1)^2);

    S=(y2-y1)/l;
    C=(x2-x1)/l;
    
    gamma=g0*(1+ki*abs(S));

    if node1==1
        F(1)=F(1)+gamma*C;
        F(2)=F(2)+gamma*S;
    
    elseif node1==2
        F(4)=F(4)+gamma*C;
        F(5)=F(5)+gamma*S;
    
    else
        F(4*(node1-1))=F(4*(node1-1))+gamma*C;
        F(4*(node1-1)+1)=F(4*(node1-1)+1)+gamma*S;
    
    end

    if node2==1
        F(1)=F(1)-gamma*C;
        F(2)=F(2)-gamma*S;
    
    elseif node2==2
        F(4)=F(4)-gamma*C;
        F(5)=F(5)-gamma*S;
    
    else
        F(4*(node2-1))=F(4*(node2-1))-gamma*C;
        F(4*(node2-1)+1)=F(4*(node2-1)+1)-gamma*S;
    
    end

end

warning off MATLAB:nearlySingularMatrix;
qdot = H\F;


%Nodal velocities

QDOT=zeros(nnodes,2);

for i=1:nnodes
    if i==1
        QDOT(i,1)=qdot(1);
        QDOT(i,2)=qdot(2);
    elseif i==2
        QDOT(i,1)=qdot(4);
        QDOT(i,2)=qdot(5);
    else
        QDOT(i,1)=qdot(4*(i-1));
        QDOT(i,2)=qdot(4*(i-1)+1);
    end
end

