%
% P5.13 Uniaxial bar with contact--total Lagrangian formulation
%
tol = 1.0e-5;  iter = 0;  E = 2000;  omega=1E6;
u = 0;         uold = u;  del=0.1;
strain = u + 0.5*u^2;
stress = E*strain;
P = stress*(1+u) ;
Fc = omega*max(del+u,0);
R = -P - Fc;
conv= R^2;
fprintf('\n iter      u1      E11      S11        Fc     conv');
fprintf('\n %3d  %8.5f  %8.5f %8.3f %8.3f %12.3e',iter,u,strain,stress,Fc,conv);
while conv > tol && iter < 20
 Kt = E*(1+u)^2 + stress + omega;
 delu = R/Kt;
 u = uold + delu;
 strain = u + 0.5*u^2;
 stress = E*strain;
 P = stress*(1+u);
 Fc = omega*max(del+u,0);
 R = -P - Fc;
 conv= R^2;
 uold = u;
 iter = iter + 1;
 fprintf('\n %3d  %8.5f  %8.5f %8.3f %8.3f %12.3e',iter,u,strain,stress,Fc,conv);
end
%
% P5.14 Uniaxial bar with contact--updated Lagrangian formulation
%
tol = 1.0e-5;  iter = 0;  E = 2000; omega=1E6;
u = 0;         uold = u;  del=0.1;
strain = u/(1+u);
stress = E*(u+.5*u^2)*(1+u);
P = stress;
Fc = omega*max(del+u,0);
R = -P - Fc;
conv= R^2;
fprintf('\n iter      u1      E11      S11        Fc     conv');
fprintf('\n %3d  %8.5f  %8.5f %8.3f %8.3f %12.3e',iter,u,strain,stress,Fc,conv);
while conv > tol && iter < 20
 Kt = E*(1+u)^2 + stress/(1+u) + omega;
 delu = R/Kt;
 u = uold + delu;
 strain = u/(1+u);
 stress = E*(u+.5*u^2)*(1+u);
 P = stress;
 Fc = omega*max(del+u,0); 
 R = -P - Fc;
 conv= R^2;
 uold = u;
 iter = iter + 1;
 fprintf('\n %3d  %8.5f  %8.5f %8.3f %8.3f %12.3e',iter,u,strain,stress,Fc,conv);
end
%
% Problem 5.16
%******************************************************************************
% Two-element contact example 2D plane strain
%******************************************************************************
XYZ=[0 0; 1 0; 0 1; 1 1; 0 2; 1 2];                         % Nodal coordinates
LE=[1 2 4 3;3 4 6 5];                                    % Element connectivity
SDISPT=[1 1 0;1 2 0;2 2 0;3 1 0;5 1 0];  % Prescribed displ. [Node, DOF, Value]
TIMS=[0.0 1.0 0.2 0.0 1.0]';      % Load increments [Start End Increment IF FF]
MID=0; PROP=[1.1538E6, 7.6923E5];    % Linear elastic material PROP=[LAMBDA MU]
XYZC=[-0.1 2; 1.1 2];                                  % Rigid body coordinates
CNODE=[5 6];                                                    % Contact nodes
LEC=[8 7];                                                    % Target segments
CDISP=[7 2 -0.2;8 2 -0.2];                   % Prescribed displ. of rigid nodes
CONTPARM=[1E8 1E7 0.1];     % Contact parameters [wN, wT, friction coefficient]
CONTACT=struct('XYZC',XYZC,'XYZP',XYZC,'CDISP',CDISP,'CNODE',CNODE,...
               'LEC',LEC,'PARM',CONTPARM);                  % Contact structure
% Set simulation parameters
PARAM=struct('ITRA',70,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);
clc; NOUT=fopen('output.txt','w');                                % Output file
NLFEA(PARAM,NOUT,MID,PROP,[],SDISPT,XYZ,LE,CONTACT);            % Calling NLFEA
%******************************************************************************
%
% Problem 5.17
%******************************************************************************
% Two-element contact example 2D plane strain
%******************************************************************************
XYZ=[0 0; 1 0; 0 1; 1 1; 0 2; 1 2];                         % Nodal coordinates
LE=[1 2 4 3;3 4 6 5];                                    % Element connectivity
SDISPT=[1 1 0;1 2 0;2 2 0;3 1 0;5 1 0];  % Prescribed displ. [Node, DOF, Value]
TIMS=[0.0 1.0 0.2 0.0 1.0]';      % Load increments [Start End Increment IF FF]
MID=0; PROP=[1.1538E6, 7.6923E5];    % Linear elastic material PROP=[LAMBDA MU]
XYZC=[-0.1 2; 1.1 2];                                  % Rigid body coordinates
CNODE=[5 6];                                                    % Contact nodes
LEC=[8 7];                                                    % Target segments
CDISP=[7 2 -0.2;8 2 -0.2];                   % Prescribed displ. of rigid nodes
CONTPARM=[1E8 1E7 0.0];     % Contact parameters [wN, wT, friction coefficient]
CONTACT=struct('XYZC',XYZC,'XYZP',XYZC,'CDISP',CDISP,'CNODE',CNODE,...
               'LEC',LEC,'PARM',CONTPARM);                  % Contact structure
% Set simulation parameters
PARAM=struct('ITRA',70,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);
NOUT=fopen('output.txt','w');                                     % Output file
NLFEA(PARAM,NOUT,MID,PROP,[],SDISPT,XYZ,LE,CONTACT);            % Calling NLFEA
%******************************************************************************
%
% Problem 5.18
%******************************************************************************
% Two-element contact example 3D
%******************************************************************************
XYZ=[0 0 0; 1 0 0; 0 1 0; 1 1 0; 0 2 0; 1 2 0;              % Nodal coordinates
     0 0 1; 1 0 1; 0 1 1; 1 1 1; 0 2 1; 1 2 1];
LE=[1 2 4 3 7 8 10 9;                                    % Element connectivity
    3 4 6 5 9 10 12 11];
SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 1 0;3 3 0;   % Prescribed displacements
        4 3 0;5 1 0;5 3 0;6 3 0;7 1 0;7 2 0;7 3 0;
        8 2 0;8 3 0;9 1 0;9 3 0;10 3 0;11 1 0;11 3 0;12 3 0];
TIMS=[0.0 1.0 0.2 0.0 1.0]';      % Load increments [Start End Increment IF FF]
MID=0; PROP=[1.1538E6, 7.6923E5];    % Linear elastic material PROP=[LAMBDA MU]
XYZC=[-0.1 2 -0.1; 1.1 2 -0.1; -0.1 2 1.1; 1.1 2 1.1]; % Rigid body coordinates
CNODE=[5 6 11 12];                                              % Contact nodes
CDISP=[13 2 -0.2;14 2 -0.2;15 2 -0.2;16 2 -0.2]; % Displacementd of rigid nodes
LEC=[13 14 16 15];                                             % Target element
CONTPARM=[0.5E8 0.5E7 0.1];      % Contact parameters [wN, wT, friction coeff.]
CONTACT=struct('XYZC',XYZC,'XYZP',XYZC,'CDISP',CDISP,'CNODE',CNODE,...
               'LEC',LEC,'PARM',CONTPARM);                  % Contact structure
% Set simulation parameters
PARAM=struct('ITRA',70,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);
NOUT=fopen('output.txt','w');                                     % Output file
NLFEA(PARAM,NOUT,MID,PROP,[],SDISPT,XYZ,LE,CONTACT);            % Calling NLFEA
%******************************************************************************
%
% Problem 5.19
%******************************************************************************
% Two-element contact example 3D
%******************************************************************************
XYZ=[0 0 0; 1 0 0; 0 1 0; 1 1 0; 0 2 0; 1 2 0;              % Nodal coordinates
     0 0 1; 1 0 1; 0 1 1; 1 1 1; 0 2 1; 1 2 1];
LE=[1 2 4 3 7 8 10 9;                                    % Element connectivity
    3 4 6 5 9 10 12 11];
SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 1 0;3 3 0;   % Prescribed displacements
        4 3 0;5 1 0;5 3 0;6 3 0;7 1 0;7 2 0;7 3 0;
        8 2 0;8 3 0;9 1 0;9 3 0;10 3 0;11 1 0;11 3 0;12 3 0];
TIMS=[0.0 1.0 0.2 0.0 1.0]';      % Load increments [Start End Increment IF FF]
MID=0; PROP=[1.1538E6, 7.6923E5];    % Linear elastic material PROP=[LAMBDA MU]
XYZC=[-0.1 2 -0.1; 1.1 2 -0.1; -0.1 2 1.1; 1.1 2 1.1]; % Rigid body coordinates
CNODE=[5 6 11 12];                                              % Contact nodes
CDISP=[13 2 -0.2;14 2 -0.2;15 2 -0.2;16 2 -0.2]; % Displacementd of rigid nodes
LEC=[13 14 16 15];                                             % Target element
CONTPARM=[0.5E8 0.5E7 0.0];      % Contact parameters [wN, wT, friction coeff.]
CONTACT=struct('XYZC',XYZC,'XYZP',XYZC,'CDISP',CDISP,'CNODE',CNODE,...
               'LEC',LEC,'PARM',CONTPARM);                  % Contact structure
% Set simulation parameters
PARAM=struct('ITRA',70,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);
NOUT=fopen('output.txt','w');                                     % Output file
NLFEA(PARAM,NOUT,MID,PROP,[],SDISPT,XYZ,LE,CONTACT);            % Calling NLFEA
%******************************************************************************
%
% Problem 5.20
%******************************************************************************
% Two-element contact example 3D
%******************************************************************************
XYZ=[0 0 0; 1 0 0; 0 1 0; 1 1 0; 0 2 0; 1 2 0;              % Nodal coordinates
     0 0 1; 1 0 1; 0 1 1; 1 1 1; 0 2 1; 1 2 1];
LE=[1 2 4 3 7 8 10 9;                                    % Element connectivity
    3 4 6 5 9 10 12 11];
SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 1 0;3 3 0;   % Prescribed displacements
        4 3 0;5 1 0;5 3 0;6 3 0;7 1 0;7 2 0;
        8 2 0;9 1 0;11 1 0];
TIMS=[0.0 1.0 0.2 0.0 1.0]';      % Load increments [Start End Increment IF FF]
MID=0; PROP=[1.1538E6, 7.6923E5];    % Linear elastic material PROP=[LAMBDA MU]
XYZC=[-0.1 2 -0.1; 1.1 2 -0.1; -0.1 2 1.1; 1.1 2 1.1]; % Rigid body coordinates
CNODE=[5 6 11 12];                                              % Contact nodes
CDISP=[13 2 -0.2;14 2 -0.2;15 2 -0.2;16 2 -0.2]; % Displacementd of rigid nodes
LEC=[13 14 16 15];                                             % Target element
CONTPARM=[0.5E8 0.5E7 0.4];      % Contact parameters [wN, wT, friction coeff.]
CONTACT=struct('XYZC',XYZC,'XYZP',XYZC,'CDISP',CDISP,'CNODE',CNODE,...
               'LEC',LEC,'PARM',CONTPARM);                  % Contact structure
% Set simulation parameters
PARAM=struct('ITRA',70,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);
NOUT=fopen('output.txt','w');                                     % Output file
NLFEA(PARAM,NOUT,MID,PROP,[],SDISPT,XYZ,LE,CONTACT);            % Calling NLFEA
%******************************************************************************