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 1994 EIJI FUKUMORI REARRANGED
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=3, MXE=10, MXN=MXE*(ND-1)+1,NBW=ND,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 ( 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
DO I = 1 , NNODE
WRITE(*,*)' NODAL # =',I, ' DISPLACEMENT =', RHS(I)
END DO
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 SYSTEM ( MXN, MBAND, NUMNP, A, B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MBAND) , B(MXN)
C---------- ELIMINATION ------------------
DO 30 N = 1 , NUMNP
DO 20 L = 2 , MBAND
C = A(N,L) / A(N,1)
I = N + L - 1
IF ( I .LE. NUMNP ) THEN
J = 0
DO K = L , MBAND
J = J + 1
A(I,J) = A(I,J) - C * A(N,K)
END DO
A(N,L) = C
B(I) = B(I) - A(N,L) * B(N)
ENDIF
20 CONTINUE
B(N) = B(N) / A(N,1)
30 CONTINUE
C---------- BACKSUBSTITUTION -------------
DO WHILE ( N .GT. 0 )
DO K = 2 , MBAND
L = N + K - 1
IF ( L .LE. NUMNP ) B(N) = B(N) - A(N,K) * B(L)
END DO
N = N - 1
END DO
RETURN
END