% This function finds the `error' in the bounded-bounded quadrature. This
% is the amount by which the ratio of shear to normal stress differs from
% the coefficient of friction at the point in the slip zone for which the
% equation was deleted. This is a measure of how incorrect the estimated
% size of the slip zone is. This function will therefore be minised by
% calling scripts.
%
% This code requires WilliamsEigenSol and PolarKernel to run and was 
% written in Matlab R2017b.
%
% Daniel Riddoch, 26-01-23
% University of Oxford, Department of Engineering Science
function [Tol]=FindTolSim(parameter,d)
%% Calculate dislocation density
[Tol,~] = CalculatePhi(parameter,d);                        % Calculating the magnitude of the error to be set as an objective function
end
% Calculation of Dislocation distribution 
%
% Calculates Dislocation Distribution to clear free surfaces from tractions
% and to enforce slip conditions, and uses this to calculate stresses along
% the slip line
% -> slip at angle theta, length d -> normal normalisation 
% -> free surfaces at angles +-alpha , length inf -> mellin transormation 
%
%                             |
%                             | r,rho,alpha
%                             |   -> v,u,X,alpha
%                             |
%                             |
%                             |
%                             |
%      _______________________|_____________ r,rho,theta,d
%   r,rho,-alpha               
%     -> v,u,X,-alpha             
%
%    A*phi = F
function [Tol,parameter] = CalculatePhi(parameter,d)
%% Parameter
N     = parameter.N;                  % Number of Gauss points along crack
theta = parameter.theta;              % Angle of slip
n     = parameter.n;                  % Number of Gauss points along free surf
facl  = parameter.facl;               % multiplier for free surf transform
alpha = parameter.alpha;              % Angle of free surfaces
norm  = parameter.norm;               % Normalisation for free notch faces
lam1  = parameter.lam1;               % Mode I eigenvalue
lam2  = parameter.lam2;               % Mode II eigenvalue
g1    = parameter.g1;                 % Realigned eigenvector
g2    = parameter.g2;                 % Realigned eigenvector
fc    = parameter.fc;                 % Coefficient of friction
ID    = parameter.ID;                 % Logical vector of equation to be discarded
%% Quadrature for slip (bounded-bounded endpoint behaviour)
% Numerical quadrature for bounded-bounded
Idx = 1:N+1;                                                               % Index up to number of points
t   = cos(pi.*(2.*Idx'-1)/(2*(N+1)));                                      % Observation (collocation) points
s   = cos(pi.*Idx(1:end-1)./(N+1));                                        % Integration points
W   = (1-s.^2)./(N+1);                                                     % Weights
om  = (1-s.^2).^(1/2);                                                     % Fundamental function
[rt,rhos,~,drhods] = normalisation(t,s,d,'normal','reverse');              % Normalising coordinates
%% Quadrature for free surfaces (singular-bounded endpoint behaviour)
idx = 1:n;                                                                 % Indices
v   = cos(pi.*(2.*idx'-1)./(2*n+1));                                       % Collocation points
u   = cos((pi.*(2.*idx))./(2*n+1));                                        % Integration points
X   = 2.*(1-u)./(2*n+1);                                                   % Weights
psi = sqrt((1-u)./(1+u));                                                  % Fundamental function
[rv,rhou,~,drhodu] = normalisation(v,u,facl,norm,'reverse');               % Normalising coordinates
%% Assemble matrix A
A = -1.*ones(N+4*n,N+4*n);    % Preallocating matrix entries for speed

% Creating matrix entries for the integral equation observed on the slip interface
[~,Krtt,Krrt,~,~,~] = PolarKernel(rt(~ID),theta,rhos,theta);                     % G_{ijk}(0<r<d,theta,0<rho<d,theta)
A(1:N    ,1:N)                   = pi.*W.*drhods.*(fc.*Krtt-Krrt);
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rt(~ID),theta,rhou,alpha);               % G_{ijk}(0<r<d,theta,0<rho<infty,alpha)
A(1:N    ,N+1:4:end)             = pi.*X.*drhodu.*(fc.*Krtt-Krrt);
A(1:N    ,N+2:4:end)             = pi.*X.*drhodu.*(fc.*Kttt-Ktrt);
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rt(~ID),theta,rhou,-alpha);              % G_{ijk}(0<r<d,theta,0<rho<infty,-alpha)
A(1:N    ,N+3:4:end)             = pi.*X.*drhodu.*(fc.*Krtt-Krrt);
A(1:N    ,N+4:4:end)             = pi.*X.*drhodu.*(fc.*Kttt-Ktrt);

% Creating matrix entries for the integral equation observed on the free
% surface theta=alpha
[~,Krtt,Krrt,~,~,~] = PolarKernel(rv,alpha,rhos,theta);                     % G_{ijk}(0<r<infty,alpha,0<rho<d,theta)
A(N+1:N+n    ,1:N)               = pi.*W.*drhods.*Krtt;
A(N+n+1:N+2*n,1:N)               = pi.*W.*drhods.*Krrt;
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rv,alpha,rhou,alpha);               % G_{ijk}(0<r<infty,alpha,0<rho<infty,alpha)
A(N+1:N+n    ,N+1:4:end)         = pi.*X.*drhodu.*Krtt;
A(N+1:N+n    ,N+2:4:end)         = pi.*X.*drhodu.*Kttt;
A(N+n+1:N+2*n,N+1:4:end)         = pi.*X.*drhodu.*Krrt;
A(N+n+1:N+2*n,N+2:4:end)         = pi.*X.*drhodu.*Ktrt;
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rv,alpha,rhou,-alpha);              % G_{ijk}(0<r<infty,alpha,0<rho<infty,-alpha)
A(N+1:N+n    ,N+3:4:end)         = pi.*X.*drhodu.*Krtt;
A(N+1:N+n    ,N+4:4:end)         = pi.*X.*drhodu.*Kttt;
A(N+n+1:N+2*n,N+3:4:end)         = pi.*X.*drhodu.*Krrt;
A(N+n+1:N+2*n,N+4:4:end)         = pi.*X.*drhodu.*Ktrt;

% Creating matrix entries for the integral equation observed on the free
% surface theta=-alpha
[~,Krtt,Krrt,~,~,~] = PolarKernel(rv,-alpha,rhos,theta);                    % G_{ijk}(0<r<infty,-alpha,0<rho<d,theta)
A(N+2*n+1:N+3*n,1:N)             = pi.*W.*drhods.*Krtt;
A(N+3*n+1:end    ,1:N)           = pi.*W.*drhods.*Krrt;
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rv,-alpha,rhou,alpha);              % G_{ijk}(0<r<infty,-alpha,0<rho<infty,alpha)
A(N+2*n+1:N+3*n,N+1:4:end)       = pi.*X.*drhodu.*Krtt;
A(N+2*n+1:N+3*n,N+2:4:end)       = pi.*X.*drhodu.*Kttt;
A(N+3*n+1:end    ,N+1:4:end)     = pi.*X.*drhodu.*Krrt;
A(N+3*n+1:end    ,N+2:4:end)     = pi.*X.*drhodu.*Ktrt;
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rv,-alpha,rhou,-alpha);             % G_{ijk}(0<r<infty,-alpha,0<rho<infty,-alpha)
A(N+2*n+1:N+3*n,N+3:4:end)       = pi.*X.*drhodu.*Krtt;
A(N+2*n+1:N+3*n,N+4:4:end)       = pi.*X.*drhodu.*Kttt;
A(N+3*n+1:end    ,N+3:4:end)     = pi.*X.*drhodu.*Krrt;
A(N+3*n+1:end    ,N+4:4:end)     = pi.*X.*drhodu.*Ktrt;

%% Right hand side
F             = zeros(N+4*n,1);   % Preallocating for speed
F(1:N,1)      = -((fc-g1).*rt(~ID).^(lam1-1)+...
                  (fc-g2).*rt(~ID).^(lam2-1));              % Slip stress condition
%% Solve system of equations
phi = A\F;

%% Find error
AID = zeros(1,N+4*n);

% Creating kernels for observation on the slip line, only for the discarded
% equation
[~,Krtt,Krrt,~,~,~] = PolarKernel(rt(ID),theta,rhos,theta);                     % G_{ijk}(0<r<d,theta,0<rho<d,theta)
AID(1    ,1:N)                   = pi.*W.*drhods.*(fc.*Krtt-Krrt);
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rt(ID),theta,rhou,alpha);               % G_{ijk}(0<r<d,theta,0<rho<infty,alpha)
AID(1    ,N+1:4:end)             = pi.*X.*drhodu.*(fc.*Krtt-Krrt);
AID(1    ,N+2:4:end)             = pi.*X.*drhodu.*(fc.*Kttt-Ktrt);
[~,Krtt,Krrt,~,Kttt,Ktrt] = PolarKernel(rt(ID),theta,rhou,-alpha);              % G_{ijk}(0<r<d,theta,0<rho<infty,-alpha)
AID(1    ,N+3:4:end)             = pi.*X.*drhodu.*(fc.*Krtt-Krrt);
AID(1    ,N+4:4:end)             = pi.*X.*drhodu.*(fc.*Kttt-Ktrt);

FID      = -((fc-g1).*rt(ID).^(lam1-1)+...
                  (fc-g2).*rt(ID).^(lam2-1));              % Slip stress condition on for the discarded equation
              
Tol=AID*phi-FID;        % The error in the quation found by subtracting the slip stress condition from the result of the integral equation
%% Parameter
% This section stores all the coordinates, densities and weight functions
% for use in other functions        
parameter.phi   = phi;
parameter.om    = om;
parameter.psi   = psi;
parameter.N     = N;
parameter.rt    = rt;
parameter.s     = s;
parameter.rhos  = rhos;
parameter.rhou  = rhou;
parameter.drhods= drhods;
parameter.drhodu= drhodu;
parameter.W     = W;
parameter.X     = X;
end
% This functions normalises the coordinates using appropriate transforms to
% the integral range [-1,1]
function [coord,intcord,dcoorddt,dintcordds] = normalisation(t,s,l,norm,rev)

switch norm
    case 'normal'   % Linear transform
        switch rev  % Choosing whether or not to reverse the integral range
            case 'reverse'
                coord      = 1./2.*l.*(t+1);
                intcord    = 1./2.*l.*(s+1);
                dcoorddt   = 1./2.*l;
                dintcordds = 1./2.*l;
            case 'normalise'
                coord      = 2.*t./l-1;
                intcord    = 2.*s./l-1;
                dcoorddt   = 0;
                dintcordds = 0;
        end
    case 'mellin'   % Mellin transform - option 1
        switch rev  % Choosing whether or not to reverse the integral rangeption 2
            case 'reverse'
                coord      = l./2.*log(2./(1-t))./log(2);
                intcord    = l./2.*log(2./(1-s))./log(2);
                dcoorddt   = l./2./((1-t).*log(2));
                dintcordds = l./2./((1-s).*log(2));
            case 'normalise'
                coord      = 1-2.^(1-2.*t./l);
                intcord    = 1-2.^(1-2.*s./l);
                dcoorddt   = 0;
                dintcordds = 0;
        end
    case 'mellin2' % Mellin transform - option 2
        switch rev  % Choosing whether or not to reverse the integral range
            case 'reverse'
                coord      = l.*log(2./(1-t));
                intcord    = l.*log(2./(1-s));
                dcoorddt   = l./(1-t);
                dintcordds = l./(1-s);
            case 'normalise'
                coord      = 1-2./exp(t./l);
                intcord    = 1-2./exp(s./l);
                dcoorddt   = 0;
                dintcordds = 0;
        end
end
end
% Kernelfunctions 
%
% Kernels for dislocations near a notch
% -> origin at the tip of the notch
% -> dislocation at r = rho and theta = phi
% Checked against StressTrans_finalDan.m
%
% Date: 23/09/20   
% 
function [Krrr,Krtt,Krrt,Ktrr,Kttt,Ktrt] = PolarKernel(r,theta,rho,phi)

R = sqrt(rho.^2 + r.^2 -2.*rho.*r.*cos(theta-phi));

Krrr = sin(theta-phi)     ./R.^4.*(2.*rho.^2.*r - r.^3 - 2.*rho.^3.*cos(theta-phi) ...
     + rho.^2.*r.*cos(2.*(theta-phi)));
Krtt =-sin(theta-phi)     ./R.^4.*(4.*rho.^2.*r + r.^3 - 2.*rho.*(rho.^2+2.*r.^2).*cos(theta-phi) ...
     + rho.^2.*r.*cos(2.*(theta-phi)));
Krrt = (r.*cos(theta-phi)-rho)./R.^4.*(          r.^2 - 2.*rho.*r.*cos(theta-phi) ...
     + rho.^2.*cos(2.*(theta-phi)));
Ktrr = (r.*cos(theta-phi).*(7.*rho.^2+2.*r.^2) - rho.*(2.*rho.^2 + 6.*r.^2 ...
     + rho.*r.*cos(3.*(theta-phi))))./R.^4./2;
Kttt = (r.*cos(theta-phi).*(5.*rho.^2+2.*r.^2) - rho.*(2.*(rho.^2 + r.^2 + 2.*r.^2.*cos(2.*(theta-phi))) ...
     - rho.*r.*cos(3.*(theta-phi))))./R.^4./2;
Ktrt = r.*sin(theta-phi)./R.^4.*(r.^2 - 2.*rho.*r.*cos(theta-phi) ...
     + rho.^2.*cos(2.*(theta-phi)));
end