% Created by ILL RYU (04/29/2010)
% -------------------------------------------------------------------
% This input is prepared for checking the traction boundary condition
% To do that, dislocation loop which normal is perpendicular to
% the cylinder top surface is made using this input. 
% 
% Using this input, ParaDiS input will be generated. 
% 

clc

clear all;

global Bscrew Bedge Beclimb Bline


R=0.5; % Radius of dislocation loop
N=100; % Number of nodes 

% base vectors
e1=[1 0 0];
e2=[0 1 0];                     % Flat loop
e2=[0 1 1];                    % Inclined loop

% Coordinate vector X
X=zeros(N,3);

% x coordinate:  vector X(:,1)
% y coordinate:  vector X(:,2)
% z coordinate:  vector X(:,3)

for i=1:N
    X(i,:)=R*(e1*cos(((i-1)/N)*2*pi)+e2*sin(((i-1)/N)*2*pi));
end

% Plots nodes 
figure(1)
limit=R*2;

plot3(X(:,1),X(:,2),X(:,3),'o');
xlim([-limit limit]);ylim([-limit limit]);zlim([-limit limit]);
xlabel('X');ylabel('Y');zlabel('Z');
grid on
%view([30 -30 40]);


% Making dislocation nodes

rn = zeros(N,4);
for i=1:3
    rn(:,i)=X(:,i);
end
rn(:,4)=0;

% Burgers vector & slip plane
b1=[1 0 0];
n1=[0 0 1];

% Making dislocation segments
% 
links = zeros(N-1,8);
for i=1:N-1
    links(i,1)=i;
    links(i,2)=i+1;
    links(i,3:5)=b1;
    links(i,6:8)=n1;

end
links(N-1,2)=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MU = 1;
NU = 0.305;
maxconnections=8;
lmax = 1000;
lmin = 200;
areamin=lmin*lmin*sin(60/180*pi)*0.5; 
areamax=20*areamin;
a=lmin/sqrt(3)*0.5;
Ec = MU/(4*pi)*log(a/0.1);
totalsteps=200;
dt0=1e7;  
mobility='mobbcc0';
%Drag (Mobility) parameters
Bscrew=1e0;
Bedge=1e0;
Beclimb=1e2;
Bline=1.0e-4*min(Bscrew,Bedge);

integrator='int_trapezoid';
rann = 0.5*a;    
rntol = 0.5*rann;    
doremesh=1;
docollision=1;
doseparation=1;
elasticinteraction=1; %line tension model when this is zero
plotfreq=1;       
plim=10000;         
%appliedstress =1e-3.*[2 0 1; 0 2 -1; 1 -1 0];
appliedstress =1.*[1 0 0; 0 0 0; 0 0 0];
viewangle=[45 -45 ];
printfreq=1;      
printnode=3;
rmax=100;
b_mag=1;
% simulation box size (for paradis)
L = 40000;
minSideX = -L/2;
minSideY = -L/2;
minSideZ = -L/2;
maxSideX =  L/2;
maxSideY =  L/2;
maxSideZ =  L/2;

% boundary conditions (for paradis)
xBoundType = 1; % free
yBoundType = 1; % free
zBoundType = 1; % free

% output and communication settings
paradis_dir = './';
%paradis_input_dir = strcat(paradis_dir,'/Runs');
paradis_input_dir = strcat(paradis_dir,'/');
paradis_output_dir = strcat('Outputs/frank_read_results'); 
