% 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
% wedge angle is changed between two specificed values
%
% 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
f     = 0.2;                                    % Coefficient of friction
% 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)
%% Varying the wedge angle
Steps=25;                                            % Number of different values of wedge angle
amin=0.725;                                            % Minimum permitted wedge angle(/pi)
amax=0.775;                                           % Maximum permitted wedge angle(/pi)
Alpha=pi.*linspace(amin,amax,Steps);                         % Vector of wedge angles
for i=1:Steps
    alpha=Alpha(i);         % Wedge angle
    theta=pi-Alpha(i);      % Interface angle
    %% 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(i) = (K1o/K2o)^(1/(lam2-lam1));                   % Characteristic length dimension
    c0d0(i)= (-(f-g1)/(f-g2))^(1/(lam2-lam1));           % Normalised slip zone length by Williams solution
    C0(i) = c0d0(i)*d0(i);                               % 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
    % An initial guess is required. This script will use 2.5 times the
    % Williams solution the first time, but to speed things up, the
    % solution for the previous alpha value is used from then on.
    if i==1
        IG=2.5*c0d0;
    else
        IG=c(i-1);
    end
    Options.Display = 'iter';
    [c(i),fval(i),exitflag(i),output] = fsolve(@(d) FindTolSeq(parameter,d),IG,Options);          % c is normalised by d0
    cRat(i)  = c(i)./c0d0(i)
    %% 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
hold on
plot(Alpha./pi,cRat,'-b','LineWidth',2)
plot(Alpha./pi,cRat.*c0d0,'-r','LineWidth',2)
xlabel(strcat('$\frac{\alpha}{\pi}$'),'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], ...
  'YLim'        , [0 10], ...
  'LineWidth'   , 1.0         );
legend({'$\frac{c}{0_0}$','$\frac{c}{d_0}$'},'interpreter','latex','FontSize',16)