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
! -------------------------------------------------------------------------
 IMPLICIT none
 INTEGER,   PARAMETER :: D_P = kind (1.d0) ! machine portable
 INTEGER,   PARAMETER :: n_incr = 5, m_iter = 20
 REAL(D_P), PARAMETER :: d = 2.5d3, x_0 = 2.5d3, area_0 = 1.d2,     &
                         E = 5.d5, f_incr = 0.3d7, c_norm = 1.d-20, &
                         tolerence = TINY (1.d0)

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

!  initialize geometry and stiffness

  x = x_0 ; area = area_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

  DO incrm = 1, n_incr
    force    = force + f_incr
    residual = residual - f_incr

!  Newton-Raphson iteration

    r_norm = c_norm * 2 ; n_iter = 0  ! initialize
    DO WHILE ( (r_norm > c_norm) .and. &
               (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)
      print "('increment = ', i3, ' iteration = ', i3,     &
      &       ' residual = ', g10.4, ' x = ', g10.4,       &
      &       ' force = ', g10.4)", incrm, n_iter, 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
        print *, 'Near zero stiffness'
        STOP     'Near zero stiffness'
      END IF ! unstable
    END DO ! Newton-Raphson iteration

  END DO ! incrm
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