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 MAY1994 EIJI FUKUMORI RETOUCHED
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER (ND=3,MXE=10,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
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 = 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