      PROGRAM BUCKLE2A
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    ********SYMM 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=ND,INTEPT=2,MXW=ND)
      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 ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
C======================================================================
C (4) SOLVE FOR UNKNOWN VARIABLES
      CALL SYSTEMQ ( MXN, MXW, NNODE, NBW, 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 = I , ND
      II = NODEX(IEL,I)
      JJ = J - I + 1
      A(II,JJ) = A(II,JJ) + STIFF(I,J)
      END DO
      END DO
      END DO
      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 33 I = 1 , ND
      DO 33 J = 1 , ND
      STIFF(I,J) = 0.
   33 CONTINUE
C------- GAUSS INTEGRATION
      DO 400 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
  400 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE FORM ( MXN, NBW, NNODE, A, RHS, IBC, BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN), A(MXN,NBW), IBC(2), BV(2), IBND(2)
C--------IBC(I) = 1 ---> DIRICHLET,  IBC(I) = 2 ---> NEUMANN
      NB = 2
      IBND(1) = 1
      IBND(2) = NNODE
      DO I = 1 , NNODE
      RHS (I) = 0.
      END DO
C------- 1ST KIND BC IMPLEMENTATION
      DO 50 K = 1 , NB
      IF ( IBC(K) .EQ. 1 ) THEN
      I = IBND(K)
      DO J = 2 , NBW 
      I = I - 1
      IF ( I.GT. 0 ) RHS(I) = RHS(I) - BV(K) * A(I,J)
      END DO
      I = IBND(K)
      DO J = 2 , NBW 
      I = I + 1 
      IF ( I .LE. NNODE ) RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
      END DO
      END IF
   50 CONTINUE
C-----REFORMATION OF MATRIX A
      DO 70 K = 1 , NB
      I = IBND (K)
      IF ( IBC(K) .EQ. 1 ) THEN
      RHS(I) = BV(K) 
      A(I,1) = 1. 
      DO J = 2 , NBW
      L = I - J + 1 
      A(I,J) = 0. 
      IF ( L  .GT. 0 ) A(L,J) = 0. 
      END DO
      END IF
      IF ( IBC(K) .EQ. 2 )  RHS (I) = RHS(I) - BV(K)
   70 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE SYSTEMQ ( MXN, MXW, NNODE, NBW, A, B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,MXW) , B(MXN)
C---------- GAUSS ELIMINATION --------------
      DO N = 1 , NNODE-1
      DO L = 2 , MIN((NNODE-N+1),NBW)
      C = A(N,L) / A(N,1)
      IF ( C .NE. 0. ) THEN
      I = N + L - 1 
      DO K = L , NBW
      J = K-L+1
      A(I,J) = A(I,J) - C * A(N,K)
      END DO
      A(N,L) = C
      B(I) = B(I) - C * B(N)
      END IF
      END DO
      B(N) = B(N) / A(N,1)
      END DO
C------ SOLUTION OF LAST UNKNOWN ---------
      B(NNODE) = B(NNODE) / A(NNODE,1)
C---------- BACK SUBSTITUTION -------------
      DO I = NNODE-1 , 1 , -1
      DO J = 2 , MIN((NNODE-I+1),NBW)
      B(I) = B(I)-A(I,J)*B(I+J-1) 
      END DO
      END DO
      RETURN
      END 
