clc;

%% define constant properties
% discretization
dz = 0.1; % length of each element
Lpile = 20; %m
nel = Lpile/dz; %m
ndof = nel + 1; % number of degrees of freedom

% pile properties
E = 33e9; %Pa
b = 0.275; %m
% >>> Need to convert square pile to circular pile
D = sqrt((4*b^2)/pi);
Ap = (pi*D^2)/4;
Ep = (E*b^2)/Ap;

% soil properties
gamma_cu = 1.4; % partial safety factor (EC DA3)
z_clay1_top = 15; %m
z_top = 20; %m
z_bottom = z_top + Lpile; %m
quake = 0.01*D;

% load properties
npiles = 108; %108
P = (671.4e3)/npiles; %N, applied load (tension)
nsteps = 150; % number of load steps
dP = P/nsteps;
loadsteps = linspace(0, P, nsteps);

%% define varying known properties
z_bottom = z_top + Lpile; %m
Su_top = (25 + (2/gamma_cu)*(z_top - z_clay1_top)); %kPa
Su_bottom = (25 + (2/gamma_cu)*(z_bottom - z_clay1_top)); %kPa
Su = linspace(Su_top, Su_bottom, nel)*1000; %Pa

%% initialize unknowns
w = zeros(1, ndof); % vertical displacements
v = zeros(1, nel); % avg. vertical displacements per element
p = zeros(1, nel); % switch parameter to determine elastic or plastic 
                   % behaviour for each element
cp = 0;
s = zeros(1, nel);
t = zeros(1, nel);
vbar = zeros(1, nel); % initialize accumulated plastic deformations per element

%% compute solution using load stepping
for ll = 2:nsteps
    % set load for current loadstep, ll
    load = loadsteps(ll);
        
    % update si and ti for all elements
    % if p == 0
    s(p == 0) = (100*Su(p == 0))/D;
    t(p == 0) = 0;
    % if p == 1
    s(p == 1) = Su(p == 1);
    t(p == 1) = Su(p == 1)*pi*D^2*dz;
    % if p == -1
    s(p == -1) = Su(p == -1);
    t(p == -1) = Su(p == -1)*pi*D^2*dz;

    % compute system of equations
    [A, b] = matrix(nel, Ep, Ap, dz, s, t, cp, vbar, load);

    % solve system of equations
    w = A\b;

    % update avg. displacements
    v = mean([w(1:end-1)';w(2:end)']);

    % accumulate plastic deformations
    vbar((v - vbar) > quake) = v((v - vbar) > quake) - quake;
    vbar((v - vbar) < -quake) = v((v - vbar) < -quake) + quake;

    % update switch parameter
    p(abs(v - vbar) <= quake) = 0;
    p((v - vbar) > quake) = 1;
    p((v - vbar) < -quake) = -1;
end

%% plotting
z = linspace(z_top, z_bottom, nel+1)';
plot(w*1000, -z);
hold on; grid on;
xlabel('Axial deformation [mm]');
ylabel('Depth below the ground surface [m]');
title('Axial deformation of a pile');