% PROGRAM truss_box_1_2
% "Introduction to Nonlinear Continuum Mechanics for Finite Element Analysis"
% J. Bonet & R.D. Wood, Cambridge University Press, ISBN 0-521-57272-X, 1997
%
% Translated from f77 syntax to f90 syntax by J. E. Akin,
% Rice University, June 2003, akin@rice.edu. 
% -------------------------------------------------------------------------
% Newton-Raphson solver for 1 d.o.f. nonlinear truss example
% Input:
%  d        = horizontal span
%  x_0      = initial height
%  area_0   = initial area
%  E        = Young's modulus
%  n_incr   = number of load increments
%  f_incr   = force increment
%  c_norm   = residual force convergence norm
%  m_iter   = max number of Newton-Raphson iterations
% -------------------------------------------------------------------------
 %B D_P = kind (1.d0) % machine portable
n_incr = 5, m_iter = 20
d = 2.5d3, x_0 = 2.5d3, area_0 = 1.d2,     
%B E = 5.d5, f_incr = 0.3d7, c_norm = 1.d-20,
E = 5.d5, f_incr = 1.5d7, c_norm = 1.d-20,
% Running with parameters E = 5.d5, f_incr = 1.5d7, c_norm = 1.d-20:
tolerence = 1e-10

% incrm            % current increment
% n_iter           % current N-R iteration number
% force = 0.d0     % current force
% vertical         % current vertical truss force
% residual = 0.d0  % current residual 
% r_norm           % current norm of residual 
% length, length_0 % current and initial length
% stiff            % current stiffness
% stress           % current stress
% x, area          % current height and area
% volume           % current volume
% u                % current displacement increment

%  initialize geometry and stiffness

x = x_0 ; area = area_0; force = 0; residual = 0;
length_0 = sqrt ( x ^2 + d ^2 );
volume   = area_0 * length_0;
stiff    = (area / length_0) * E * (x / length_0 ) ^2;

%  start load increments

%B DO incrm = 1, n_incr
for  incrm = 1:2 % n_incr
  force    = force + f_incr ;
  residual = residual - f_incr ;

%  Newton-Raphson iteration

  r_norm = c_norm * 2 ; n_iter = 0 ;  % initialize
  %B DO WHILE ( (r_norm > c_norm) .and. &
  while ( (r_norm > c_norm) & (n_iter < m_iter) )
    n_iter = n_iter + 1 ;

%     Find geometry increment

    u      = -residual / stiff ;      % displacement increment
    x      = x + u ;                  % update height
    length = sqrt ( x ^2 + d ^2 ) ;   % update length
    area   = volume / length        ; % update area
      
%     Find stresses and residual forces

    stress   = E * log (length / length_0) ;  % axial stress
    vertical = stress * area * x / length ;   % vertical force
    residual = vertical - force ;
    r_norm   = abs (residual / force) ;
    %B print "('increment = ', i3, ' iteration = ', i3,     ...
    %B &       ' residual = ', g10.4, ' x = ', g10.4,       ...
    %B &       ' force = ', g10.4)", incrm, n_iter, r_norm, ...
    %B &       x, force
    fprintf ('increment = %i iteration = %i', incrm, n_iter)
    fprintf ('residual = %9.2e, x = %9.2e, force = %9.2e \n', r_norm, x, force)

%     Find stiffness and check stability

    stiff = (area/length) * (E - 2 * stress) * (x/length) ^2 ...
          + (stress * area / length) ;
    if ( abs ( stiff) < tolerence ) % THEN
      disp ('Near zero stiffness')
      error   ('Near zero stiffness')
    end %  % unstable
  end % DO % Newton-Raphson iteration

end % DO % incrm
%B END PROGRAM truss_box_1_2
%
% Running with parameters E = 5.d5, f_incr = 1.5d7, c_norm = 1.d-20:
% increment = 1 iteration = 1 residual = 0.2184     x = 4621. force = 0.1500E+08
% increment = 1 iteration = 2 residual = 0.4239E-01 x = 5540. force = 0.1500E+08
% increment = 1 iteration = 3 residual = 0.2973E-02 x = 5822. force = 0.1500E+08
% increment = 1 iteration = 4 residual = 0.1806E-04 x = 5844. force = 0.1500E+08
% increment = 1 iteration = 5 residual = 0.6780E-09 x = 5845. force = 0.1500E+08
% increment = 1 iteration = 6 residual = 0.1242E-15 x = 5845. force = 0.1500E+08
% increment = 1 iteration = 7 residual =  0.000     x = 5845. force = 0.1500E+08

%B increment = 1, iteration = 1 residual =  2.184e-01, x =  4.621e+03,  force =  1.500e+07
%B increment = 1, iteration = 2 residual =  4.239e-02, x =  5.540e+03,  force =  1.500e+07
%B increment = 1, iteration = 3 residual =  2.973e-03, x =  5.822e+03,  force =  1.500e+07
%B increment = 1, iteration = 4 residual =  1.806e-05, x =  5.844e+03,  force =  1.500e+07
%B increment = 1, iteration = 5 residual =  6.780e-10, x =  5.845e+03,  force =  1.500e+07
%B increment = 1, iteration = 6 residual =  1.242e-16, x =  5.845e+03,  force =  1.500e+07
%B increment = 1, iteration = 7 residual =  0.000e+00, x =  5.845e+03,  force =  1.500e+07
%B increment = 2, iteration = 1 residual =  4.287e-01, x =  1.364e+04,  force =  3.000e+07
%B increment = 2, iteration = 2 residual =  1.407e+00, x = -3.147e+04,  force =  3.000e+07
%B increment = 2, iteration = 3 residual =  1.105e+00, x = -2.353e+05,  force =  3.000e+07
%B increment = 2, iteration = 4 residual =  1.012e+00, x = -3.482e+06,  force =  3.000e+07
%B increment = 2, iteration = 5 residual =  1.000e+00, x = -3.568e+08,  force =  3.000e+07
%B increment = 2, iteration = 6 residual =  1.000e+00, x = -2.054e+12,  force =  3.000e+07
