% This script determines the calibration matrix entries for a single crack
% placed in a wedge of internal angle 2*alpha, at an angle theta. This
% script allows user to vary the angle of the crack theta between two
% predetermined points. The calibration matrix entries for each calculated
% point are generated, and a figure plotting the stress intensity factor
% ratios is produced.
%
% A numerical error can occur depending on the choice of N, n and facl.
% To mitigate this a solution is patched in, to correct the affected values
% of theta.
% 
% This script was written in MATLAB R2017b and requires SingleCrackCalib
% and WilliamsEigenSol
%
% The results produced are displayed in figure 3. Numerical values may 
% differ from those published if different parameter values are used.
% 
% Daniel Riddoch, 2021-03-01
% University of Oxford, Department of Engineering Science
clc
clear
close all
%% Input vectors
% This section contains all the inputs to the problem, both physical and
% numerical
Pts=100;                                    % This is the number of points that will be calculated
alpha=0.75*pi;                              % alpha, the wedge half-angle(radians)
thetaLower=0;                               % Lower limit of theta(radians)
thetaUpper=0.4*pi;                          % Upper limit of theta(radians)
theta=linspace(thetaLower,thetaUpper,Pts);  % theta range vector(radians)
N=100;                                      % Number of points on the crack face(default value is 100)
n=500;                                      % Number of Gauss points for notch faces(default value is 500)
facl=80;                                    % Length of free surfaces for notch faces(default value is 80)
%% Calibration loop
% This section determines the calibration matrix entries and stores them in
% the Mat matrix

% Initialising matrices
Mat=zeros(Pts,4);
EigenVectors=zeros(Pts,2);
% Looping for theta values
for i=1:Pts
    [a,b,c,d,EigenVec]=SingleCrackCalib(alpha,theta(i),1,n,facl,N);      % Generating calibration entries
    % Storing values
    Mat(i,1)=a;
    Mat(i,2)=b;
    Mat(i,3)=c;
    Mat(i,4)=d;
    EigenVectors(i,1)=EigenVec.thetatheta1;
    EigenVectors(i,2)=EigenVec.thetatheta2;
end
%% Patch
% This section corrects a numerical error which occurs due to the small
% number of points and low free surface length. The constants N, n and facl
% are increased in value for s small region to correct this defect without
% slowing the code overall. Changing the default values of N,n or facl cmay
% change the range over which this defect is visible, requiring a change of
% PthetaLower and/or PthetaUpper.
Pts=10;                                    % This is the number of points that will be calculated
PthetaLower=0.24*pi;                         % Lower limit of alpha(radians)
PthetaUpper=0.28*pi;                          % Upper limit of alpha(radians)
Ptheta=linspace(PthetaLower,PthetaUpper,Pts);  % alpha range vector(radians)
N=500;                                      % Number of points on the crack face(default value is 500)
n=2500;                                      % Number of Gauss points for notch faces(default value is 2500)
facl=250;                                    % Length of free surfaces for notch faces(default value is 250)
%% Calibration loop
% This section determines the calibration matrix entries and stores them in
% the Mat matrix for the patch

% Initialising matrices
PMat=zeros(Pts,4);
EigenVectors=zeros(Pts,2);
% Looping for alpha values
for i=1:Pts
    [a,b,c,d,EigenVec]=SingleCrackCalib(alpha,Ptheta(i),1,n,facl,N);      % Generating calibration entries
    % Storing values
    PMat(i,1)=a;
    PMat(i,2)=b;
    PMat(i,3)=c;
    PMat(i,4)=d;
    EigenVectors(i,1)=EigenVec.thetatheta1;
    EigenVectors(i,2)=EigenVec.thetatheta2;
end
%% Inserting patch
% This section splices the corrected patch results into the original set of
% results
CutBelow=find(theta<PthetaLower);
CutAbove=find(theta>PthetaUpper);
PatchT=[theta(CutBelow),Ptheta,theta(CutAbove)];
PatchM1=[Mat(CutBelow,1);PMat(:,1);Mat(CutAbove,1)];
PatchM2=[Mat(CutBelow,2);PMat(:,2);Mat(CutAbove,2)];
PatchM3=[Mat(CutBelow,3);PMat(:,3);Mat(CutAbove,3)];
PatchM4=[Mat(CutBelow,4);PMat(:,4);Mat(CutAbove,4)];
%% Plotting results
% This section plots the calibration matrix entries
Font=36;
figure
hold on
plot(PatchT./pi,PatchM1,'--r','LineWidth',2)
plot(PatchT./pi,PatchM2,'--b','LineWidth',2)
plot(PatchT./pi,PatchM3,':r','LineWidth',2)
plot(PatchT./pi,PatchM4,':b','LineWidth',2)
axis([thetaLower/pi thetaUpper/pi -2.5 2.5])
xlabel('$\frac{\theta}{\pi}$','interpreter','latex','FontSize',Font)
ylabel('Matrix coefficient','interpreter','latex','FontSize',Font)
legend({'a','b','c','d'},'interpreter','latex','FontSize',Font)
set(gca,'FontSize',24)