%
% Example 3.4 Polar decomposition of simple shear deformation
F=[1 2/sqrt(3);0 1];
C=F'*F;
[P. L]=eig(C);
U=P*sqrt(L)*P';
Q=F*inv(U);
V=F*Q';

X1=[0 0]';
X11 = U*X1;
X12 = Q*X11;
X2=[1 0]';
X21 = U*X2;
X22 = Q*X21;
X3=[1 1]';
X31 = U*X3;
X32 = Q*X31;
X4=[0 1]';
X41 = U*X4;
X42 = Q*X41;

px1=[X1(1) X2(1) X3(1) X4(1) X1(1)];
py1=[X1(2) X2(2) X3(2) X4(2) X1(2)];
line(px1,py1);
hold on;
px2=[X11(1) X21(1) X31(1) X41(1) X11(1)];
py2=[X11(2) X21(2) X31(2) X41(2) X11(2)];
line(px2,py2);
px3=[X12(1) X22(1) X32(1) X42(1) X12(1)];
py3=[X12(2) X22(2) X32(2) X42(2) X12(2)];
line(px3,py3);

figure;
X1=[0 0]';
X11 = Q*X1;
X12 = V*X11;
X2=[1 0]';
X21 = Q*X2;
X22 = V*X21;
X3=[1 1]';
X31 = Q*X3;
X32 = V*X31;
X4=[0 1]';
X41 = Q*X4;
X42 = V*X41;

px1=[X1(1) X2(1) X3(1) X4(1) X1(1)];
py1=[X1(2) X2(2) X3(2) X4(2) X1(2)];
line(px1,py1);
hold on;
px2=[X11(1) X21(1) X31(1) X41(1) X11(1)];
py2=[X11(2) X21(2) X31(2) X41(2) X11(2)];
line(px2,py2);
px3=[X12(1) X22(1) X32(1) X42(1) X12(1)];
py3=[X12(2) X22(2) X32(2) X42(2) X12(2)];
line(px3,py3);
%
% Example 3.8 Uniaxial bar--total Lagrangian formulation
tol=1.0e-5; iter=0; E=200; f = 100; u=0; uold=u;  % Parameters & initialization  
strain=u+0.5*u^2; stress=E*strain; % Calculate Lagrange strain & 2nd P-K stress
P = stress*(1+u); R = f - P;                        % Internal force & residual
conv= R^2/(1+f^2);                                      % Convergence criterion
fprintf('\n iter      u1      E11      S11         conv');         % Print head
fprintf('\n %3d  %7.5f  %7.5f %8.3f %12.3e %7.5f',iter,u,strain,stress,conv); %
while conv > tol && iter < 20                      % Convergence iteration loop
 Kt = E*(1+u)^2 + stress; delu = R/Kt;  % Tangent stiffness & displacement inc.
 u = uold + delu;                                         % Update displacement
 strain=u+0.5*u^2; stress=E*strain;% Calculate Lagrange strain & 2nd P-K stress
 P = stress*(1+u); R = f - P;                       % Internal force & residual
 conv= R^2/(1+f^2);                                     % Convergence criterion
 uold = u; iter = iter + 1;     % Store displacement and move to next iteration
 fprintf('\n %3d  %7.5f  %7.5f %8.3f %12.3e %7.5f',iter,u,strain,stress,conv);%
end                                                                           %
%
% Example 3.9 Uniaxial bar--updated Lagrangian formulation
tol=1.0e-5; iter=0; E=200; f = 100; u=0; uold=u;  % Parameters & initialization  
strain=u/(1+u); stress=E*(u+.5*u^2)*(1+u);   % Calculate strain & Cauchy stress
P = stress; R = f - P;                              % Internal force & residual
conv= R^2/(1+f^2);                                      % Convergence criterion
fprintf('\n iter      u1      E11      S11         conv');         % Print head
fprintf('\n %3d  %7.5f  %7.5f %8.3f %12.3e %7.5f',iter,u,strain,stress,conv); %
while conv > tol && iter < 20                      % Convergence iteration loop
 Kt = E*(1+u)^2 + stress/(1+u); delu = R/Kt;  % Tangent stiffness & displ. inc.
 u = uold + delu;                                         % Update displacement
 strain=u/(1+u); stress=E*(u+.5*u^2)*(1+u);  % Calculate strain & Cauchy stress
 P = stress; R = f - P;                             % Internal force & residual
 conv= R^2/(1+f^2);                                     % Convergence criterion
 uold = u; iter = iter + 1;     % Store displacement and move to next iteration
 fprintf('\n %3d  %7.5f  %7.5f %8.3f %12.3e %7.5f',iter,u,strain,stress,conv);%
end                                                                           %
%
% Example 3.16 Hyperelastic tension example
XYZ=[0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1];      % Nodal coordinates
LE=[1 2 3 4 5 6 7 8];                                    % Element connectivity
EXTFORCE=[];                                                % No external force
SDISPT=[1 1 0;4 1 0;5 1 0;8 1 0;                              % u1=0 for Face 6
        1 2 0;2 2 0;5 2 0;6 2 0;                              % u2=0 for Face 3
        1 3 0;2 3 0;3 3 0;4 3 0;                              % u3=0 for Face 1
        2 1 5;3 1 5;6 1 5;7 1 5];                             % u1=5 for Face 4
TIMS=[0.0 1.0 0.05 0.0 1.0]';                 % Load increments [t0 T Dt f0 fL]
MID=-1; PROP=[80   20   1E7];               % Material properties [A10, A01, K]
PARAM=struct('ITRA',30,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);         %
clc; NOUT=fopen('output.txt','w');                           % Open output file
NLFEA(PARAM,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE,[],[]);        % Calling NLFEA
fclose(NOUT);                                               % Close output file
[DMAX,DMIN]=PLOTELEM(XYZ,LE,'output.txt',1.0,0,1)      % Plot deformed geometry
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% P3.20 Two bar elements--displacement controlled procedure
tol = 1.0e-5;  conv = 0;  u2 = 0;  u2old = u2;
E = 1E8; A1 = 1E-4; A2 = .5E-4;
fprintf('\n step      u2      u3        F');
% Displacement increment loop
for i=1:10
 u3 = 0.1*i;
 S1 = E*(u2+.5*u2^2);
 S2 = E*(u3-u2+.5*(u3-u2)^2);
 P = S1*A1*(1+u2)+S2*A2*(u2-u3-1);
 F = S2*A2*(1+u3-u2);
 R = -P;
 conv = R^2;
 % Convergence loop
 iter = 0;
 while conv > tol && iter < 50
  Kt = A1*(E*(1+u2)^2+S1) + A2*(E*(u2-u3-1)^2+S2);
  delu2 = R/Kt;
  u2 = u2old + delu2;
  S1 = E*(u2+.5*u2^2);
  S2 = E*(u3-u2+.5*(u3-u2)^2);
  P = S1*A1*(1+u2)+S2*A2*(u2-u3-1);
  R = -P;
  conv= R^2;
  u2old = u2;
  iter = iter + 1;
 end
 F = S2*A2*(1+u3-u2);
 fprintf('\n %3d  %7.5f %7.5f %8.3f',i,u2,u3,F);
end
%
% P3.21 Uniaxial bar--updated Lagrangian formulation
%
tol = 1.0e-5;  iter = 0;  E = 200;
u = 0;         uold = u;  f = 100;
strain = u/(1+u);
stress = E*strain;
P = stress;
R = f - P;
conv= R^2/(1+f^2);
fprintf('\n iter      u1    Strain   Stress         conv');
fprintf('\n %3d  %7.5f  %7.5f %8.3f %12.3e %7.5f',iter,u,strain,stress,conv);
while conv > tol && iter < 20
 Kt = E/(1+u)^2;
 delu = R/Kt;
 u = uold + delu;
 strain = u/(1+u);
 stress = E*strain;
 P = stress;
 R = f - P;
 conv= R^2/(1+f^2);
 uold = u;
 iter = iter + 1;
 fprintf('\n %3d  %7.5f  %7.5f %8.3f %12.3e %7.5f',iter,u,strain,stress,conv);
end