PROGRAM BUCKLE1C
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 ****** SYMMETRIC TRI-DIAGONAL MATRIX SOLVER******************
C JUNE 16, 1999 EIJI FUKUMORI
C DOWN AND UPPER STREAM SIDE ELEMENTS: NEW-TWO-NODE-PARABOLIC.
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=2, MXE=100, MXN=MXE+1,NBW=ND )
DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN),
* IBTYPE(2), BV(2)
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,ND,NBW,P,NE,NNODE,NODEX,EI,X,A,RHS,IBTYPE )
C======================================================================
C (3) IMPLEMENTATION OF BOUNDARY CONDITION
CALL FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
C======================================================================
C (4) SOLVE FOR UNKNOWN VARIABLES
CALL SYSTEM ( MXN, NBW, NNODE, 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=',G16.7, '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='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,
* IBTYPE)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN),
* IBTYPE(2)
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)
RL = X(J) - X(I)
ALPHA = P / EI(IEL)
C-------- TOP ELEMENT
IF ( IEL .EQ. 1 ) THEN
IF ( IBTYPE(1) .EQ. 1 ) THEN
A(I,1) = A(I,1) - 1./RL + ALPHA*RL/3.
A(I,2) = A(I,2) + 1./RL + ALPHA*RL/6.
A(J,1) = A(J,1) - 1./RL + ALPHA*RL/3.
END IF
IF ( IBTYPE(1) .EQ. 2 ) THEN
A(I,1) = A(I,1) - 4./(3.*RL) + ALPHA*RL*8./15.
A(I,2) = A(I,2) + 4./(3.*RL) + ALPHA*RL*2./15.
A(J,1) = A(J,1) - 4./(3.*RL) + ALPHA*RL*3./15.
END IF
END IF
C------- LAST ELEMENT
IF ( IEL .EQ. NE ) THEN
IF ( IBTYPE(2) .EQ. 1 ) THEN
A(I,1) = A(I,1) - 1./RL + ALPHA*RL/3.
A(I,2) = A(I,2) + 1./RL + ALPHA*RL/6.
A(J,1) = A(J,1) - 1./RL + ALPHA*RL/3.
END IF
IF ( IBTYPE(2) .EQ. 2 ) THEN
A(I,1) = A(I,1) - 4./(3.*RL) + ALPHA*RL*3./15.
A(I,2) = A(I,2) + 4./(3.*RL) + ALPHA*RL*2./15.
A(J,1) = A(J,1) - 4./(3.*RL) + ALPHA*RL*8./15.
END IF
END IF
C-------- IN BETWEEN
IF ( (IEL .GT. 1) .AND. (IEL .LT. NE) ) THEN
A(I,1) = A(I,1) - 1./RL + ALPHA*RL/3.
A(I,2) = A(I,2) + 1./RL + ALPHA*RL/6.
A(J,1) = A(J,1) - 1./RL + ALPHA*RL/3.
END IF
C
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