%
% Example 4.4 Elastoplastic deformation of two bars in parallel
A1=.75; A2=1.25;
E1=10000; E2=5000;
Et1=1000; Et2=500;
sY1=5; sY2=7.5;
P=15;L1=100;L2=100;
H1=E1*Et1/(E1-Et1);
H2=E2*Et2/(E2-Et2);

iter = 1

du=P/(E1*A1/L1+E2*A2/L2)
de=du/L1
u = du;

ds1=E1*de
str1=ds1;
R=1-abs(str1-sY1)/ds1
s1=R*ds1+Et1*(1-R)*de
ep1=(1-R)/(1+H1/E1)*abs(de)

ds2=E2*de
str2=ds2;
s2 = str2

Residual = P - (s1*A1+s2*A2)

iter = 2

Eep1 = R*E1 + (1-R)*Et1
du=Residual/(Et1*A1/L1+E2*A2/L2)
de=du/L1
u = u + du

ds1=E1*de
str1=s1 + ds1
sYmax = sY1 + H1*ep1
sYmin = sYmax - 2*sY1
R=0
s1=s1+Et1*de
ep1=ep1+1/(1+H1/E1)*abs(de)

ds2=E2*de
str2=s2 + ds2
R2 = 1-(abs(str2)-sY2)/ds2
s2=s2+R2*ds2+Et2*(1-R2)*de
ep2=(1-R2)/(1+H2/E2)*abs(de)

Residual = P - (s1*A1+s2*A2)

iter = 3

du=Residual/(Et1*A1/L1+Et2*A2/L2)
de=du/L1
u = u + du

ds1=E1*de
str1=s1 + ds1
sYmax = sY1 + H1*ep1
sYmin = sYmax - 2*sY1
R=0
s1=s1+Et1*de
ep1=ep1+1/(1+H1/E1)*abs(de)

ds2=E2*de
str2=s2 + ds2
sY = sY2 + H2*ep2
R2 = 0;
s2=s2+Et2*de
ep2=ep2+1/(1+H2/E2)*abs(de)

Residual = P - (s1*A1+s2*A2)
%
% Example 4.5 Elastoplastic deformation of two bars in parallel
E1=10000; Et1=1000; sY1=5;E2=5000; Et2=500; sY2=7.5;      % Material properties
A1=0.75; L1=100; A2=1.25; L2=100;                        % Geometric parameters
tol=1.0E-5; u=0; P=15; iter=0;                           % Initialize variables
nS1=0; nAlp1=0; nep1=0; nS2=0; nAlp2=0; nep2=0;          % Initialize variables
mp1=[E1, 1, E1*Et1/(E1-Et1), sY1];                   % Bar1 material properties
mp2=[E2, 0, E2*Et2/(E2-Et2), sY2];                   % Bar2 material properties
Res=P-nS1*A1-nS2*A2;                                         % Initial residual
Dep1=E1; Dep2=E2;                                   % Initial tangent stiffness
conv=Res^2/(1+P^2);                                     % Convergence criterion
fprintf('\niter        u      S1      S2  Alpha1  Alpha2');        % Print head
fprintf('      ep1      ep2   Residual');                          % Print head
fprintf('\n %3d  %7.4f %7.3f %7.3f %7.3f %7.3f %8.6f %8.6f %10.3e',...        %
      iter,u,nS1,nS2,nAlp1,nAlp2,nep1,nep2,Res);                % Print results
while conv>tol && iter<20                       % Newton-Raphson iteration loop
  delu=Res/(Dep1*A1/L1 + Dep2*A2/L2);        % Calculate displacement increment
  delE=delu/L1; u=u+delu;            % Strain increment and update displacement
  [S1,Alp1,ep1]=COMBHARD1D(mp1,delE,nS1,nAlp1,nep1);      % State determination
  [S2,Alp2,ep2]=COMBHARD1D(mp2,delE,nS2,nAlp2,nep2);      % State determination
  Res=P-S1*A1-S2*A2;                                                 % Residual
  conv=Res^2/(1+P^2);                                   % Convergence criterion
  iter=iter+1;                                         % Move to next iteration
  Dep1=E1; if ep1>nep1; Dep1=Et1; end                       % Tangent stiffness
  Dep2=E2; if ep2>nep2; Dep2=Et2; end                       % Tangent stiffness
  nS1=S1; nAlp1=Alp1; nep1=ep1;                        % Update state variables
  nS2=S2; nAlp2=Alp2; nep2=ep2;                        % Update state variables
  fprintf('\n %3d  %7.4f %7.3f %7.3f %7.3f %7.3f %8.6f %8.6f %10.3e',...      %
      iter,u,nS1,nS2,nAlp1,nAlp2,nep1,nep2,Res);                % Print results
end                                                                           %
%
% Example 4.16 Shear deformation of elastoplastic square
E=24000; NU=0.2; MU=E/2/(1+NU); LAMBDA=NU*E/((1+NU)*(1-2*NU)); % Material prop.
MP = [LAMBDA MU 0 1000 200*sqrt(3)];        % Elastoplastic material properties
IDEN=[1 1 1 0 0 0]';                                % 2nd-order identity matrix
D=MU*diag([2 2 2 1 1 1]) + LAMBDA*(IDEN*IDEN');             % Elastivity matrix
STRESS=zeros(6,1); ALPHA=zeros(6,1); EP=0;       % Initialize history variables
DEPS= [0 0 0 0.004 0 0]';                                    % Strain increment
X=zeros(1,15); Y=zeros(1,16);                              % Strain-stress data
for i=2:16                                                % Load increment loop
 [STRESS,ALPHA,EP]=COMBHARD(MP,D,DEPS,STRESS,ALPHA,EP,false);% Elastoplasticity
 X(i)=(i-1)*DEPS(4); Y(i)=STRESS(4);            % Store strain & stress history
end                                                % End of load increment loop
plot(X,Y,'-o');                                  % Plot the strain-stress graph
%
% Example 4.20 - Shear deformation of a square (finite rotation)
E=24000; NU=0.2; MU=E/2/(1+NU); LAMBDA=NU*E/((1+NU)*(1-2*NU)); % Material prop.
MP = [LAMBDA MU 0 1000 200*sqrt(3)];        % Elastoplastic material properties
IDEN=[1 1 1 0 0 0]';                                % 2nd-order identity matrix
D=MU*diag([2 2 2 1 1 1]) + LAMBDA*(IDEN*IDEN');             % Elasticity matrix
STRESS=zeros(6,1); ALPHA=zeros(6,1); EP=0;  % Initialize small strain variables
STRESSR=STRESS; ALPHAR=ALPHA; EPR=EP;     % Initialize large rotation variables
L=[0 0.024 0;-0.02 0 0;0 0 0];                              % Velocity gradient
DEPS= [0 0 0 0.004 0 0]';                                    % Strain increment
X=zeros(1,15); Y1=zeros(1,16); Y2=zeros(1,16);             % Strain-stress data
for i=2:16                                                % Load increment loop
 [STRESSR,ALPHAR]=ROTATESTRESS(L,STRESSR,ALPHAR);               % Rotate stress
 [STRESSR,ALPHAR,EPR]=COMBHARD(MP,D,DEPS,STRESSR,ALPHAR,EPR,false);           %
 [STRESS,ALPHA,EP]=COMBHARD(MP,D,DEPS,STRESS,ALPHA,EP,false);% Elastoplasticity
 X(i)=(i-1)*DEPS(4); Y1(i)=STRESS(4); Y2(i)=STRESSR(4); % Store strain & stress
end                                                % End of load increment loop
plot(X,Y1,'k-o',X,Y2,'k--');                     % Plot the strain-stress graph
%
% Example 4.23 - Shear deformation of a square
E=24000; NU=0.2; MU=E/2/(1+NU); LAMBDA=NU*E/((1+NU)*(1-2*NU)); % Material prop.
MP = [LAMBDA MU 0 1000 200*sqrt(3)];        % Elastoplastic material properties
IDEN=[1 1 1 0 0 0]';                                % 2nd-order identity matrix
D=MU*diag([2 2 2 1 1 1]) + LAMBDA*(IDEN*IDEN');             % Elasticity matrix
DM=2*MU*eye(3) + LAMBDA*[1 1 1]'*[1 1 1]; % Elasticity matrix (principal space)
L=[0 0.024 0;-0.02 0 0;0 0 0];                              % Velocity gradient
DEPS= [0 0 0 0.004 0 0]';                                    % Strain increment
STRN=zeros(6,1); ALPN=zeros(6,1); EPN=0;    % Initialize small strain variables
STRR=STRN; ALPR=ALPN; EPR=EPN;            % Initialize large rotation variables
STRM=STRN; ALPM=[0 0 0]';EPM=EPN;     % Initialize finite deformation variables
BM=[1 1 1 0 0 0]';                    % Elastic Cauchy-Green deformation tensor
X=zeros(1,16);Y1=zeros(1,16);Y2=zeros(1,16);Y3=zeros(1,16);% Strain-stress data
for i=2:16                                                % Load increment loop
 [STRR,ALPR]=ROTATESTRESS(L,STRR,ALPR);                         % Rotate stress
 [STRR,ALPR,EPR]=COMBHARD(MP,D,DEPS,STRR,ALPR,EPR,false);    % Elastoplasticity
 [STRN,ALPN,EPN]=COMBHARD(MP,D,DEPS,STRN,ALPN,EPN,false);    % Elastoplasticity
 [STRM,BM,ALPM,EPM]=MULPLAST(MP,DM,L,BM,ALPM,EPM,false);     % Elastoplasticity
 X(i)=(i-1)*DEPS(4);Y1(i)=STRN(4);Y2(i)=STRR(4);Y3(i)=STRM(4);     % Store data
end                                                % End of load increment loop
plot(X,Y1,'k-',X,Y2,'k-o',X,Y3,'k--');           % Plot the strain-stress graph
%
% Example 5.4 Two bar in parallel--Ramberg-Osgood model
%
A1=0.75; A2=1.25; E1=10000; E2=5000; n1=10; n2=7; sY1=5; sY2=7.5;
tol = 1.0e-10; iter = 0; u = 0; uold = u; f = 15; L1 = 100; L2 = 100;
s1=0;s2=0;ep=0;
P = s1*A1+s2*A2;
R = f - P;
conv= R^2/(1+f^2);
fprintf('\n iter      u     s1       s2         conv');
 fprintf('\n %3d  %7.5f %7.5f %7.5f %12.3e',iter,u,s1,s2,conv);
while conv > tol && iter < 20
 Eep1 = E1/(1+3/7*n1*(s1/sY1)^(n1-1));
 Eep2 = E2/(1+3/7*n2*(s2/sY2)^(n2-1));
 Kt = Eep1*A1/L1 + Eep2*A2/L2;
 delu = R/Kt;
 delep = delu/L1;
 ep = ep + delep;
 u = uold + delu;
%calculate stress using local N-R iteration
 residual = s1/E1 + 3/7*sY1/E1*(s1/sY1)^n1 - ep;
 k=0;
 while residual^2 > 1E-13
  dels1 = -Eep1*residual;
  s1 = s1 + dels1;
  residual = s1/E1 + 3/7*sY1/E1*(s1/sY1)^n1 - ep;
  Eep1 = E1/(1+3/7*n1*(s1/sY1)^(n1-1));
  k = k+1;
 end
 residual = s2/E2 + 3/7*sY2/E2*(s2/sY2)^n2 - ep;
 k=0;
 while residual^2 > 1E-13
  dels2 = -Eep2*residual;
  s2 = s2 + dels2;
  residual = s2/E2 + 3/7*sY2/E2*(s2/sY2)^n2 - ep;
  Eep2 = E2/(1+3/7*n2*(s2/sY2)^(n2-1));
  k = k + 1;
 end
 P = s1*A1+s2*A2;
 R = f - P;
 conv= R^2/(1+f^2);
 uold = u;
 iter = iter + 1;
 fprintf('\n %3d  %7.5f %7.5f %7.5f %12.3e',iter,u,s1,s2,conv);
end
%
% Figure 5.7  Stress-strain curve for the Ramberg-Osgood model
%
E = 10000;
sy = 40;
n = 11;
s = -50:1:50;
sr = s/sy;
eps = s/E + 3*sy/(7*E)*sr.^n;
plot(eps,s,'k-');
e0 = eps(101);
s0 = s(101);

sr=(s-s0)/(sy+abs(s0));
eps = e0 + (s-s0)/E + 3*(sy+abs(s0))/(7*E)*sr.^n;
hold on;
plot(eps,s,'b--');

s = -s;
eps = -eps;
plot(eps,s,'b--');

axis([-0.03 0.03 -55 55]);
line([-0.03 0.03],[0 0]);
line([0 0],[-55 55]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% P4.1 elastoplastic bar (MPa, mm)
%
delE=0.003; A=100;
mp = [1E5, 0, 1E4, 100];
nS=0; nA=0; nep=0;
[Snew, Anew, epnew]=combHard1D(mp,delE,nS,nA,nep)
eelast=delE - epnew
Force = Snew*A
%
% P4.4 Elastoplastic bar (isotropic hardening)
%
delE=-0.003; nS=200; nA=0; nep=1E-4;
mp = [2E5, 0, 2.5E4, 250];
[Snew, Anew, epnew]=combHard1D(mp,delE,nS,nA,nep)
%
% P4.5 Elastoplastic bar (kinematic hardening)
%
delE=-0.003; nS=200; nA=2.5; nep=1E-4;
mp = [2E5, 1, 2.5E4, 250];
[Snew, Anew, epnew]=combHard1D(mp,delE,nS,nA,nep)
%
% P4.6 Elastoplastic bar (combined hardening)
%
delE=-0.003; nS=200; nA=2.5; nep=1E-4;
mp = [2E5, 0.5, 2.5E4, 250];
[Snew, Anew, epnew]=combHard1D(mp,delE,nS,nA,nep)
%
% P4.8 Variable loadings
%
E=70000; H=10000;Et=1000; sYield=250;
mp = [E, 0, H, sYield];
Et=E*H/(E+H);
nS = 0; nep = 0; nA=0; epnew=0;
A = 100; L = 1000;
tol = 1.0E-5; u = 0; iter=0; Res=0;
Force = 1000*[5:5:30 25 20 25 30 35 30 25 20];
N = size(Force',1);
X=zeros(N,1);Y=zeros(N,1);
fprintf('\nstep  iter        u       S       ep   Residual');
fprintf('\n %3d %3d  %7.4f %7.3f %8.6f %10.3e',...
      i,iter,u,nS,nep,Res);
for i=1:N
 P = Force(i); iter = 0;
 Res = P - nS*A;
 conv = Res^2/(1+P^2);
 du=0;
 while conv > tol && iter < 20
   Eep = E; if epnew > nep; Eep = Et; end
   delu = Res / (Eep*A/L);
   du = du + delu;
   delE = du / L;
   [Snew, Anew, epnew]=combHard1D(mp,delE,nS,nA,nep);
   Res = P - Snew*A;
   conv = Res^2/(1+P^2);
   iter = iter + 1;
 end
 u=u+du;
 nS=Snew; nep=epnew;
 X(i) = u; Y(i) = nS;
 fprintf('\n %3d %3d  %7.4f %7.3f %8.6f %10.3e',...
      i,iter,u,nS,nep,Res);
end
X=[0;X];Y=[0;Y];plot(X,Y);
%
% P4.9 Variable displacement
%
E=70000; H=10000; sYield=250;
mp = [E, 0, H, sYield];
nS = 0; nA=0; nep = 0;
A = 100; L = 1000;
u = 0; iter=0; Res=0;
disp=[0 1 2 3 4 5 4 3 4 5 6 7 6];
N = size(disp',1);
X=zeros(N,1);Y=zeros(N,1);
fprintf('\nstep         u       S       ep');
fprintf('\n %3d %7.4f %7.3f %8.6f',i,u,nS,nep);
for i=2:N
 delu = disp(i) - disp(i-1);
 delE = delu / L;
 [Snew, Anew, epnew]=combHard1D(mp,delE,nS,nA,nep);
 nS = Snew; nep = epnew;
 X(i) = disp(i); Y(i) = nS;
 fprintf('\n %3d %7.4f %7.3f %8.6f',i,u,nS,nep);
end
plot(X,Y);
%
% P4.10 Two-bar in parallel - unloading
%
E1=10000; Et1=1000; sYield1=5;
E2=5000;  Et2=500;  sYield2=7.5;
mp1 = [E1, 0, E1*Et1/(E1-Et1), sYield1];
mp2 = [E2, 1, E2*Et2/(E2-Et2), sYield2];
nS1 = 0; nA1=0; nep1 = 0; epnew1=0;
nS2 = 0; nA2=0; nep2 = 0; epnew2=0;
A1 = 0.75; L1 = 100;
A2 = 1.25; L2 = 100;
tol = 1.0E-5; u = 0;
Force = [1:15 14:-1:0];
N = size(Force',1);
X=zeros(N,1);Y1=zeros(N,1);Y2=zeros(N,1);
fprintf('\nstep iter        u      S1      S2      ep1      ep2   Residual');
fprintf('\n %3d %3d  %7.4f %7.3f %7.3f %8.6f %8.6f %10.3e',...
      0,0,u,nS1,nS2,nep1,nep2,0);
for i=1:N
 P = Force(i); iter = 0;
 Res = P - nS1*A1 - nS2*A2;
 conv = Res^2/(1+P^2);
 while conv > tol && iter < 20
   Eep1 = E1; if epnew1>nep1; Eep1 = Et1; end
   Eep2 = E2; if epnew2>nep2; Eep2 = Et2; end
   delu = Res / (Eep1*A1/L1 + Eep2*A2/L2);
   u = u + delu;
   delE = delu / L1;
   [Snew1, Anew1, epnew1]=combHard1D(mp1,delE,nS1,nA1,nep1);
   [Snew2, Anew2, epnew2]=combHard1D(mp2,delE,nS2,nA2,nep2);
   Res = P - Snew1*A1 - Snew2*A2;
   conv = Res^2/(1+P^2);
   iter = iter + 1;
   nS1 = Snew1; nep1 = epnew1; nA1=Anew1;
   nS2 = Snew2; nep2 = epnew2; nA2=Anew2;
 end
 X(i) = u; Y1(i) = nS1; Y2(i) = nS2;
 fprintf('\n %3d %3d  %7.4f %7.3f %7.3f %8.6f %8.6f %10.3e',...
      i,iter,u,nS1,nS2,nep1,nep2,Res);
end
X=[0;X];Y1=[0;Y1];Y2=[0;Y2];plot(X,Y1,X,Y2);
%
% P4.11 Two-bar in serial
%
clear;
E=100000; H=10000; sYield=100;Et=E*H/(E+H);
mp = [E, 0, H, sYield];
nS1 = 0; nA1=0; nep1 = 0; epnew1=0; Eep1=E;
nS2 = 0; nA2=0; nep2 = 0; epnew2=0; Eep2=E;
A1 = 100; L1 = 100; A2 = 50; L2 = 100;
tol = 1.0E-8; iter = 0; u = [0 0]'; F = [0 12000]';
Res = F - [nS1*A1-nS2*A2;nS2*A2];
conv = norm(Res)^2/(1+norm(F)^2);
fprintf('\n iter      u1      u2      S1      S2      ep1      ep2   Residual');
fprintf('\n %3d  %7.4f %7.4f %7.3f %7.3f %8.6f %8.6f %10.2e %10.2e',...
      0,u(1),u(2),nS1,nS2,nep1,nep2,0);
while conv > tol && iter < 20
  Kt = [Eep1*A1/L1+Eep2*A2/L2,-Eep2*A2/L2;-Eep2*A2/L2,Eep2*A2/L2];
  delu = Kt\Res;
  u = u + delu;
  delE1 = delu(1) / L1;
  delE2 = (delu(2)-delu(1)) / L2;
  [Snew1, Anew1, epnew1]=combHard1D(mp,delE1,nS1,nA1,nep1);
  [Snew2, Anew2, epnew2]=combHard1D(mp,delE2,nS2,nA2,nep2);
  Eep1 = E; if epnew1>nep1; Eep1 = Et; end
  Eep2 = E; if epnew2>nep2; Eep2 = Et; end
  nS1 = Snew1; nep1 = epnew1; nA1=Anew1;
  nS2 = Snew2; nep2 = epnew2; nA2=Anew2;
  Res = F - [nS1*A1-nS2*A2;nS2*A2];
  conv = norm(Res)^2/(1+norm(F)^2);
  iter = iter + 1;
  fprintf('\n %3d  %7.4f %7.4f %7.3f %7.3f %8.6f %8.6f %10.2e %10.2e',...
      iter,u(1),u(2),nS1,nS2,nep1,nep2,Res);
end
%
% P4.17 shear deformation of a square
%
lambda=1000; mu=1000;
mp = [lambda, mu, 0, 100, 100];
Iden=[1 1 1 0 0 0]';
deps=[0 0 0 0.01 0 0]';
stressN=[0 0 0 50 0 0]';
alphaN=[0 0 0 0 0 0]';
epN=0;
D=2*mu*eye(6) + lambda*Iden*Iden';
D(4,4) = mu; D(5,5) = mu; D(6,6) = mu;
[stress, alpha, ep]=combHard(mp,D,deps,stressN,alphaN,epN)
%
% P4.18 shear deformation of a square
%
lambda=100; mu=100;
mp = [lambda, mu, 0, 10, sqrt(12)];
Iden=[1 1 1 0 0 0]';
stressN=[0 0 0 0 0 0]';
alphaN=[0 0 0 0 0 0]';
epN=0;
D=2*mu*eye(6) + lambda*Iden*Iden';
D(4,4) = mu; D(5,5) = mu; D(6,6) = mu;
deps=[0 0 0 0.016 0 0]';
[stress, alpha, ep]=combHard(mp,D,deps,stressN,alphaN,epN);
stressN=stress; alphaN=alpha;epN=ep;
deps=[0 0 0 0.008 0 0]';
[stress, alpha, ep]=combHard(mp,D,deps,stressN,alphaN,epN);
%
% P4.19 tension and shear of a square
%
lambda=100000; mu=100000;
mp = [lambda, mu, 0.5, 10000, 100];
Iden=[1 1 1 0 0 0]';
stressN=[100 0 0 0 0 0]';
alphaN=[0 0 0 0 0 0]';
epN=0;
D=2*mu*eye(6) + lambda*Iden*Iden';
D(4,4) = mu; D(5,5) = mu; D(6,6) = mu;
deps=[0 0 0 0.002 0 0]';
[stress, alpha, ep]=combHard(mp,D,deps,stressN,alphaN,epN)
%
% P4.22 shear deformation of a square (saturated isotropic hardening)
%
lambda=1000; mu=1000;epinf=0.05;Y0=100;Yinf=200;
mp = [lambda, mu, epinf, Y0, Yinf];
Iden=[1 1 1 0 0 0]';
deps=[0 0 0 0.01 0 0]';
stressN=[0 0 0 50 0 0]';
alphaN=[0 0 0 0 0 0]';
epN=0;
D=2*mu*eye(6) + lambda*Iden*Iden';
D(4,4) = mu; D(5,5) = mu; D(6,6) = mu;
[stress,ep,D]=SATHARD(mp,D,deps,stressN,epN,false)