      PROGRAM BUCKLE2
C======================================================================
C AN FEM SOLVER FOR BUCKLING PROBLEM WITH MOMENT AT BOTH ENDS OF BEAM
C        EQUATION: D2U/DXDX + ALPHA*U=0; ALPHA = P/(EI)
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    ************* PENTA-DIAGONAL MATRIX SOLVER******************
C    **************** 3-NODED PARABOLIC ELEMENT USED ***************
C MAY 16, 1994   EIJI FUKUMORI
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER (ND=3,MXE=100,MXN=MXE*(ND-1)+1,NBW=2*(ND-1)+1,INTEPT=2)
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN),
     * IBTYPE(2), BV(2), STIFF(ND,ND),SAI(INTEPT),W(INTEPT),
     * F0(ND), F1(ND), SF(ND,INTEPT), BP(ND,INTEPT), B(ND),SX(ND)
C======================================================================
      DATA SAI / -0.5773502691896D0, 0.5773502691896D0 /
      DATA W   /  1.D0             , 1.D0              /
C======================================================================
      CALL DERIV ( ND, INTEPT, F0, F1, SAI, BP )
      CALL SHAPEF( ND, INTEPT, F0, SAI, SF )
C======================================================================
C (1) READING OF DATA
      CALL INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV )
C======================================================================
C (2) CONSTRUCTION OF FEM-MATRIX EQUATION
      CALL MATRIX ( MXE,MXN,INTEPT,ND,NBW,P,NE,NNODE,STIFF,
     *                    NODEX,EI,X,A,RHS, BP,W,SX,B,SF )
C======================================================================
C (3) IMPLEMENTATION OF BOUNDARY CONDITION
      CALL FORM ( ND, MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
C======================================================================
C (4) SOLVE FOR UNKNOWN VARIABLES
      IPD = ND - 1
      CALL SYSTEMA ( MXN , NBW , NNODE , IPD , A , RHS )
C======================================================================
C (5) PRINTING RESULTS
      OPEN ( 1,FILE='SOLUTION.FEM', STATUS='UNKNOWN')
      WRITE (*,*) 'SOLUTION IS IN A FILE OF SOLUTION.FEM.'
      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 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='BUCKLE2.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*(ND-1) + 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,INTEPT,ND,NBW,P,NE,NNODE,STIFF,
     *                    NODEX,EI,X,A,RHS, BP,W,SX,B,SF )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN),
     * STIFF(ND,ND),BP(ND,INTEPT),W(INTEPT),SX(ND),B(ND),SF(ND,INTEPT)
      DO I = 1 , NNODE
      DO J = 1 , NBW
      A(I,J) = 0.
      END DO
      RHS(I) = 0.
      END DO
      DO IEL = 1 , NE
      SX(1) = X(NODEX(IEL,1))
      SX(2) = X(NODEX(IEL,2))
      SX(3) = X(NODEX(IEL,3))
      ALPHA = P / EI(IEL)
      CALL SGSM ( INTEPT,ND,BP,W,SX,B,SF,ALPHA, STIFF )
      DO I = 1 , ND
      DO J = 1 , ND
      II = NODEX(IEL,I)
      JJ = ND - I + J
      A(II,JJ) = A(II,JJ) + STIFF(I,J)
      END DO
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE FORM ( ND, 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
      DO J = ND+1, NBW
      A(1,J) = 0.
      END DO
      A(1,ND) = 1.
      RHS(1) = BV(1)
      ELSE
      RHS(1) = RHS(1) - BV(1)
      END IF
      IF ( IBTYPE(2) .EQ. 1 ) THEN
      DO J = 1 , ND-1
      A(NNODE,J) = 0.
      END DO
      A(NNODE,ND) = 1.
      RHS(NNODE) = BV(2)
      ELSE
      RHS(NNODE) = RHS(NNODE) - BV(2)
      END IF
      RETURN
      END
C
C
      SUBROUTINE DERIV ( ND, INTEPT, F0, F1, SAI, BPP )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),BPP(ND,INTEPT), F0(ND),F1(ND)
C------- COMPUTATION OF BP(J) = D N(J) / D ETA
      DO K = 1 , INTEPT
      CALL ISOPARA ( ND , SAI(K)+0.5 , F1 )
      CALL ISOPARA ( ND , SAI(K)-0.5 , F0 )
      DO I = 1 , ND
      BPP(I,K) = F1(I) - F0(I)
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE SHAPEF ( ND , INTEPT , F , SAI , SF )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION F(ND) , SAI(INTEPT) , SF(ND,INTEPT)
      DO K = 1 , INTEPT
      CALL ISOPARA ( ND , SAI(K), F )
      DO I = 1 , ND
      SF(I,K) = F(I)
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE ISOPARA ( ND , SAI , F )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION F(ND)
      F(1) =   0.5 * SAI * ( SAI - 1.)
      F(2) = ( 1.- SAI ) * ( 1.+ SAI )
      F(3) =   0.5 * SAI * ( SAI + 1.)
      RETURN
      END
C
C
      SUBROUTINE SGSM ( INTEPT,ND,BP,W,X,B,SF,ALPHA, STIFF ) 
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION X(ND),W(INTEPT),STIFF(ND,ND),B(ND),BP(ND,INTEPT),
     *          SF(ND,INTEPT)
C-------- RESET
      DO I = 1 , ND
      DO J = 1 , ND
      STIFF(I,J) = 0.
      END DO
      END DO
C------- GAUSS INTEGRATION
      DO K = 1 , INTEPT
      YACOB = 0.
      DO I = 1 , ND
      YACOB = YACOB + BP(I,K)*X(I)
      END DO
      DO J = 1 , ND
      B(J) = BP(J,K) / YACOB
      END DO
      DO I = 1 , ND
      DO J = 1 , ND
      STIFF(I,J) = STIFF(I,J) + W(K)*YACOB *
     *            ( -B(I)*B(J) + ALPHA*SF(I,K)*SF(J,K) )
      END DO
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE SYSTEMA ( MXN , MXW , N , IPD , A , C )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,MXW), C(MXN)
C------- FULL BANDWIDTH = 2 * IPD + 1
      N1 = N - 1
      IPU = IPD + 1
      NIP = N - IPD
      DO K = 1 , N1
      IP = IPD
      IF ( K .GT. NIP ) IP = N - K
      IS = K + 1
      IE = K + IP
      DO I = IS , IE
      JS = IPD - ( I - IS )
      JE = JS + IP
      P = A(I,JS) / A(K,IPU)
      IF ( P .NE. 0. ) THEN
      DO J = JS , JE
      L = J - JS + IPU
      A(I,J) = A(I,J) - P * A(K,L)
      END DO
      C(I) = C(I) - P * C(K)
      END IF
      END DO
      END DO
C------- SOLUTION OF X(N)
      C(N) = C(N) / A(N,IPU)
C------- BACK SUBSTITUTION
      JS = IPU + 1
      DO K = 1 , N1
      I = N1 - K + 1
      NI = N - I
      IF ( NI .GT. IPD ) NI = IPD
      JE = IPU + NI
      T = 0.
      DO J = JS , JE
      L = J - JS + I + 1
      T = T + A(I,J) * C(L)
      END DO
      C(I) = ( C(I) - T ) / A(I,IPU)
      END DO
      RETURN
      END
  
