% 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
% coefficient of friction is changed between two given values and the
% change in the slip-stick boundary is 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
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
%% Parameter
% Creating parameter structure to transfer variables
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;
%% Varying the coefficient of friction
Steps=25;                                            % Number of different values of coefficient of friction(Default is 25)
% Minimum permitted coefficient of friction(Note, small values here produce errors, default is 0.2)
fmin=0.2;       
% Maximum permitted coefficient of friction(Recall critical coefficient of friction and try not to go too close, default is 0.45)
fmax=0.45;                                           
f=linspace(fmin,fmax,Steps);                         % Vector of coefficient of friction values
for i=1:Steps
    c0d0 = (-(f(i)-g1)/(f(i)-g2))^(1/(lam2-lam1));              % Normalised slip zone length by Williams solution
    C0   = c0d0*d0;                                       % True slip zone size by Williams solution
    % Stroing the coefficient of frition in the parameter structure.
    parameter.f     = f(i);
    %% Find slip zone size
    Options.Display = 'iter';
    [c,fval,exitflag(i),output] = fsolve(@(d) FindTolSeq(parameter,d),2.5*c0d0,Options)
    cRat(i)  = c./c0d0
    %% Slip interface properties
    % If desired this section can calculate and store the interfacial
    % tractions and slip displacements for each different coefficient of
    % friction. To do so simply un-comment this section
%     % Calculate dislocation density
%     [Phi,om,W,t,s] = CalculatePhiSeq(parameter,c);
%     % Calculating stresses
%     parameter.t     = t;
%     parameter.s     = s;
%     parameter.W     = W;
%     parameter.om    = om;
%     parameter.Phi   = Phi;
%     [Stress]  =  FindStresses(parameter,c);
%     SlipInterface(:,i) = Stress.x./c;
%     SlipDisp(:,i)      = Stress.CSlip;
end
%% Plotting
figure
plot(f,cRat,'LineWidth',2)
xlabel(strcat('$f$'),'interpreter','latex','FontSize',16)
ylabel(strcat('$\frac{c}{c_0}$'),'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         );