      PROGRAM BUCKLE1B
C======================================================================
C          --------- NON-LINEAR DIFFERENTIAL EQUATION ---------
C AN FEM SOLVER FOR BUCKLING PROBLEM WITH MOMENT AT BOTH ENDS OF BEAM
C                  EQUATION: D2U/DXDX + ALPHA*U=0; 
C                ALPHA = P/(EI)*(1+(DU/DX)**2)**1.5
C P=APPLIED FORCE AT ENDS OF BEAM, E=YOUNG MODULUS, I=2ND MOMENT OF
C INERTIA. RL=LENGTH OF ELEMENT, IBTYPE(1)=BOUNDARY CONDITION AT THE
C LEFT END OF BEAM, IBTYPE(2)=BC AT RIGHT END. IBTYPE(I)=1 FOR FIXED 
C DISPLACEMENT, IBTYPE(I)=2 FOR PRESCRIBED SLOPE AT THE END.
C    ****** SYMMETRIC TRI-DIAGONAL MATRIX SOLVER******************
C                  DECEMBER 25, 1998   EIJI FUKUMORI
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=2, MXE=100, MXN=MXE+1,NBW=ND )
      PARAMETER ( MXITERA=80,ERRMAX=1.D-6 )
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN),
     * IBTYPE(2), BV(2), U(MXN)
C======================================================================
      CALL INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV )
      CALL INITIAL (MXN, NNODE, U,X,IBTYPE,BV )
C======================================================================
      WRITE(*,*) 'ITERATION MAXIMUM % ERROR'
      ITERA = 0
      UMAXERR = 2.*ERRMAX
      DO WHILE ( ITERA .LE. MXITERA .AND. UMAXERR .GT. ERRMAX )
      ITERA = ITERA + 1
      CALL MATRIX ( MXE,MXN,ND,NBW,P,NE,NNODE,NODEX,EI,X,A,RHS,U )
      CALL FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
      CALL SYSTEM ( MXN, NBW, NNODE, A, RHS )
      CALL NEWU ( MXN,NNODE,RHS ,U, UMAXERR )
      write (*,*) ITERA , UMAXERR
      END DO
C======================================================================
C (5) PRINTING RESULTS
      OPEN ( 1,FILE='SOLUTION.FEM', STATUS='UNKNOWN')
      DO I = 1 , NNODE
      WRITE(1,100) I, X(I), RHS(I)
  100 FORMAT ( ' NODAL #=', I3, '  X=',F13.6, ' DISPLACEMENT=',G20.11 )
      END DO
      CLOSE (1)
      STOP' NORMAL TERMINATION'
      END
C
C
      SUBROUTINE NEWU ( MXN,NNODE,RHS ,U, UMAXERR )
C------- COMPUTES RADIUS AT TX INFINITE STRENGTH.-------
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION U(MXN), RHS(MXN)
C------- MIN AND MAX OF RHS
      RHSMAX = RHS(1)
      RHSMIN = RHS(1)
      DO I = 2 , NNODE
      RHSMAX = DMAX1 ( RHSMAX , RHS(I) )
      RHSMIN = DMIN1 ( RHSMIN , RHS(I) )
      END DO
C------- MAXIMUM DIFFERENCE IN RHS
      BASE = RHSMAX - RHSMIN
      IF ( BASE .EQ. 0.) BASE = 1.
C------- FIND MAX DIFFRENCE BETWEEN RHS(I) AND U(I)
      UMAXERR = 0.
      DO I = 1 , NNODE
      UMAXERR = DMAX1 ( UMAXERR , DABS ( RHS(I) - U(I) )  )
      U(I) = RHS(I)
      END DO
C-------- PERCENT ERROR IN U(I)
      UMAXERR = UMAXERR / BASE * 100.
      RETURN
      END
C
C
      SUBROUTINE INITIAL (MXN, NNODE, U,X,IBTYPE,BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION U(MXN),X(MXN),IBTYPE(2), BV(2)
C--------FIND MAX X AND MIN X
      XMAX = X(1)
      XMIN = X(1)
      DO I = 2 , NNODE
      XMAX = DMAX1 ( XMAX , X(I) )
      XMIN = DMIN1 ( XMIN , X(I) )
      END DO
C--------LENGTH OF DOMAIN
      DL = XMAX - XMIN 
C--------CASE OF DIRICHILET - DIRICHILET
      IF ( IBTYPE(1) .EQ. 1 .AND. IBTYPE(2).EQ.1 ) THEN
      SLOPE = (BV(2)-BV(1))/DL
      DO I = 1 , NNODE
      U(I) = SLOPE*(X(I)-XMIN) + BV(1)
      END DO
      END IF
C--------CASE OF DIRICHILET - NEUMANN
      IF ( IBTYPE(1) .EQ. 1 .AND. IBTYPE(2).EQ.2 ) THEN
      DO I = 1 , NNODE
      U(I) =  BV(1)
      END DO
      END IF
C--------CASE OF NEUMANN - DIRICHILET
      IF ( IBTYPE(1) .EQ. 2 .AND. IBTYPE(2).EQ.1 ) THEN
      DO I = 1 , NNODE
      U(I) =  BV(2)
      END DO
      END IF
      RETURN
      END
C
C
      SUBROUTINE INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),IBTYPE(2), BV(2)
      OPEN ( 1,FILE='BUCKLE.DAT', STATUS='OLD')
      READ(1,*) P
      READ(1,*) NE
      DO I = 1 , NE
      READ(1,*) IEL, (NODEX(IEL,J),J=1,ND),EI(IEL)
      END DO
      NNODE = NE + 1
      DO I = 1 , NNODE
      READ(1,*) NODE, X(NODE)
      END DO
      READ(1,*) IBTYPE(1), BV(1)
      READ(1,*) IBTYPE(2), BV(2)
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE MATRIX (MXE,MXN,ND,NBW,P,NE,NNODE,NODEX,EI,X,A,RHS,U)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN),U(MXN)
      DO I = 1 , NNODE
      A(I,1) = 0.
      A(I,2) = 0.
      RHS(I) = 0.
      END DO
      DO IEL = 1 , NE
      I = NODEX(IEL,1)
      J = NODEX(IEL,2)
      DX = X(J) - X(I)
      DUDX = (U(J)-U(I))/DX
      ALPHA = P/EI(IEL)*(1.+DUDX*DUDX)**1.5
      A(I,1) = A(I,1) - 1./DX + ALPHA*DX/3.
      A(I,2) = A(I,2) + 1./DX + ALPHA*DX/6.
      A(J,1) = A(J,1) - 1./DX + ALPHA*DX/3.
      END DO
      RETURN
      END
C
C
      SUBROUTINE FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION IBTYPE(2),BV(2),A(MXN,NBW), RHS(MXN)
      IF ( IBTYPE(1) .EQ. 1 ) THEN
      A(1,1) = 1.
      RHS(1) = BV(1)
      RHS(2) = RHS(2) - BV(1)*A(1,2)
      A(1,2) = 0.
      ELSE
      RHS(1) = RHS(1) - BV(1)
      END IF
      IF ( IBTYPE(2) .EQ. 1 ) THEN
      A(NNODE,1) = 1.
      RHS(NNODE) = BV(2)
      RHS(NNODE-1) = RHS(NNODE-1) - BV(2)*A(NNODE-1,2)
      A(NNODE-1,2) = 0.
      ELSE
      RHS(NNODE) = RHS(NNODE) - BV(2)
      END IF
      RETURN
      END
C
C
      SUBROUTINE SYSTEM ( MXN, NBW, NNODE, A , B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,NBW) , B(MXN)
      B(1) = B(1) / A(1,1)
      A(1,1) = A(1,2) / A(1,1)
      DO I = 2 , NNODE
      P = A(I,1) - A(I-1,2) * A(I-1,1)
      A(I,1) = A(I,2) / P
      B(I) = ( B(I) - A(I-1,2)*B(I-1) ) / P
      END DO
C------ BACK SUBSTITUTION ----
      DO I = NNODE-1, 1,-1
      B(I) = B(I) - A(I,1) * B(I+1)
      END DO
      RETURN
      END
