% This fuction calculates the slip stick boundary for a complete contact.
% This is done by optimising to minimise the error found by reintroducing
% the discarded equation from the bounded-bounded quadrature. The tractions 
% along the contact interface are then plotted. 
%
% This code requires WilliamsEigenSol, CalculatePhiSeq, FindStresses and 
% FindTolSeq to run and was written in Matlab R2017b.
%
% Daniel Riddoch, 03/06/21
% Department of Engineering Science, University of Oxford
close all
clear
clc
%% Parameter specification
% Stress intensity factors
K1    = -1;
K2    = 1;
% Contact properties
alpha = 0.75*pi;                                 % Half angle of notch
f     = 0.3;                                     % Coefficient of friction
theta = pi-alpha;                                % Angle of crack
% Crack properties
N     = 2000;                                        % Number of Gauss points(default value is 2000)
ID    = round(N/2);                                  % The number of the equation to exclude
ID    = logical([zeros(ID-1,1);1;zeros(N-ID+1,1)]);  % Logical vector to determine which equation to discard
% Input for notch faces
n     = 10000;                                       % Number of Gauss points(default value is 10000)
facl  = 150;                                         % Multiplicative constant(default value is 150)
%% Williams solution
[EigenVal,EigenVec]=WilliamsEigenSol(alpha,theta);  % Williams eigensolution
lam1 = EigenVal.M1;                                 % Williams eigenvalue - mode I   
lam2 = EigenVal.M2;                                 % Williams eigenvalue - mode II
g1   = EigenVec.rtheta1/EigenVec.thetatheta1;       % Realigned eigenvector
g2   = EigenVec.rtheta2/EigenVec.thetatheta2;       % Realigned eigenvector
K1o  = EigenVec.thetatheta1*K1;                     % Realigned eigenvalue
K2o  = EigenVec.thetatheta2*K2;                     % Realigned eigenvalue
d0   = (K1o/K2o)^(1/(lam2-lam1));                   % Characteristic length dimension
c0d0 = (-(f-g1)/(f-g2))^(1/(lam2-lam1));            % Normalised slip zone length by Williams solution
C0   = c0d0*d0;                                     % True slip zone size by Williams solution
%% Parameter
% Creating parameter structure to transfer variables
parameter.f     = f;
parameter.N     = N;
parameter.ID    = ID;
parameter.theta = theta;
parameter.n     = n;
parameter.facl  = facl;
parameter.alpha = alpha;
parameter.lam1  = lam1;
parameter.lam2  = lam2;
parameter.g1    = g1;
parameter.g2    = g2;
%% Find slip zone size
% Calculate the optimal slip zone size to minimise the error
Options.Display = 'iter';
[c,fval,exitflag,output] = fsolve(@(d) FindTol(parameter,d),2*c0d0,Options)
cRat  = c./c0d0;
% Calculate dislocation density
[Phi,om,W,t,s] = CalculatePhi(parameter,c);
% Calculating stresses
parameter.t     = t;
parameter.s     = s;
parameter.W     = W;
parameter.om    = om;
parameter.Phi   = Phi;
[Stress]  =  FindStresses(parameter,c);
%% Plots
% This section creates the plot of stresses and displacements
figure
subplot(1,2,1)
cla
hold on
plot(Stress.x,f.*Stress.N,'b-','LineWidth',1.5)
plot(Stress.x,Stress.S,'r-','LineWidth',1.5)
plot(Stress.x,Stress.S./Stress.N/f,'k-','LineWidth',1.5)
plot(Stress.x,f.*Stress.Nb,'b--','LineWidth',1.5)
plot(Stress.x,Stress.Sb,'r--','LineWidth',1.5)
plot(Stress.x,Stress.Sb./Stress.Nb/f,'k--','LineWidth',1.5)
% Axis labels
xlabel(strcat('$x/d_{0}$'),'interpreter','latex','FontSize',16)
ylabel(strcat('$\frac{\sigma_{i\theta}}{K_{I}^{o}d_{0}^{\lambda_{I}-1}}$'),'interpreter','latex','FontSize',16)
ax = gca;
set(ax, ...
  'Box'         , 'on'     , ...
  'TickDir'     , 'in'     , ...
  'TickLength'  , [.02 .02] , ...
  'TickLabelInterpreter'  , 'latex' , ...
  'FontSize'    , 16 , ...
  'XMinorTick'  , 'off'      , ...
  'YMinorTick'  , 'off'      , ...
  'XGrid'       , 'on'      , ...
  'YGrid'       , 'on'      , ...
  'XColor'      , [0 0 0], ...
  'YColor'      , [0 0 0], ...
  'YLim'        , [0 5], ...
  'LineWidth'   , 1.0         );
legend({'$p(r)$-Slip','$q(r)$-Slip','$\frac{q(r)}{fp(r)}$-Slip',...
    '$p(r)$-Stuck','$q(r)$-Stuck','$\frac{q(r)}{fp(r)}$-Stuck'},'interpreter','latex','FontSize',16)
    
subplot(1,2,2)
cla
hold on
plot(Stress.x,Stress.CSlip,'k-','LineWidth',1.5)
% Axis labels
xlabel(strcat('$x/d_{0}$'),'interpreter','latex','FontSize',16)
ylabel(strcat('$\frac{2\mu}{\pi(\kappa+1)}\frac{h(x)}{K_{I}^{o}d_{0}^{\lambda_{I}-1}}$'),'interpreter','latex','FontSize',16)
ax = gca;
set(ax, ...
  'Box'         , 'on'     , ...
  'TickDir'     , 'in'     , ...
  'TickLength'  , [.02 .02] , ...
  'TickLabelInterpreter'  , 'latex' , ...
  'FontSize'    , 16 , ...
  'XMinorTick'  , 'off'      , ...
  'YMinorTick'  , 'off'      , ...
  'XGrid'       , 'on'      , ...
  'YGrid'       , 'on'      , ...
  'XColor'      , [0 0 0], ...
  'YColor'      , [0 0 0], ...
  'LineWidth'   , 1.0         );