% 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 wedge alpha 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 alpha.
% 
% This script was written in MATLAB R2017b and requires SingleCrackCalib
% and WilliamsEigenSol
%
% The results produced are displayed in figure 2. 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
alphaLower=0.75*pi;                         % Lower limit of alpha(radians)
alphaUpper=0.9*pi;                          % Upper limit of alpha(radians)
alpha=linspace(alphaLower,alphaUpper,Pts);  % alpha range vector(radians)
theta=pi-alpha;                             % theta value, note that this can depend on alpha, or be independent(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 alpha values
for i=1:Pts
    [a,b,c,d,EigenVec]=SingleCrackCalib(alpha(i),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
% PalphaLower and/or PalphaUpper.
Pts=10;                                    % This is the number of points that will be calculated
PalphaLower=0.78*pi;                         % Lower limit of alpha(radians)
PalphaUpper=0.79*pi;                          % Upper limit of alpha(radians)
Palpha=linspace(PalphaLower,PalphaUpper,Pts);  % alpha range vector(radians)
Ptheta=pi-Palpha;                             % theta value, note that this can depend on alpha, or be independent(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(Palpha(i),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(alpha<PalphaLower);
CutAbove=find(alpha>PalphaUpper);
PatchAP=[alpha(CutBelow),Palpha,alpha(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(PatchAP./pi,PatchM1,'--r','LineWidth',2)
plot(PatchAP./pi,PatchM2,'--b','LineWidth',2)
plot(PatchAP./pi,PatchM3,':r','LineWidth',2)
plot(PatchAP./pi,PatchM4,':b','LineWidth',2)
axis([alphaLower/pi alphaUpper/pi -2.5 2.5])
xlabel('$\frac{\alpha}{\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)