Hello,
I have a piece of FEA code that came with an old text book I bought that I cannot get to compile. I wrote codes in Fortran 90/95 in college but I have never used Fortran 77. If someone can please help me get this to compile. I have attached the code and the input file that it reads from. The error that it is end-of-file during read. I tried a few things that I thought might work but I cannot get to work. It's probably something stupid that I am over looking.
Thanks for everybodys help
C ----------------------- FEM1D2 ------------------------
DIMENSION X(50),NOC(25,2),F(50),AREA(25),MAT(25),DT(25),
$ PM(10,2),NU(20),U(20),MPC(15,2),BT(15,3)
DIMENSION S(100,55)
C IMAX = FIRST DIMENSION OF THE S-MATRIX
CHARACTER*16 FILE1,FILE2
CHARACTER*81 DUMMY,TITLE
PRINT *, '***************************************'
PRINT *, '* PROGRAM FEM1D2 *'
PRINT *, '* WITH MULTI-POINT CONSTRAINTS *'
PRINT *, '* T.R.Chandrupatla and A.D.Belegundu *'
PRINT *, '***************************************'
IMAX=100
PRINT *,'Input Data File Name'
READ '(A)',FILE1
LINP=10
OPEN(UNIT=10,FILE=FILE1,STATUS='UNKNOWN',ACTION='READ')
IF(IOSTAT /=0)THEN
WRITE(*,*)'INPUT FILE DID NOT OPEN'
ELSE
WRITE(*,*)'IOSTAT', IOSTAT
END IF
PRINT *,'Output Data File Name'
READ(*,*),FILE2
LOUT=11
OPEN(UNIT=11,FILE=FILE2,STATUS='UNKNOWN',ACTION='WRITE')
if(IOSTAT /= 0)then
write(*,*)"OUTPUT file did not open"
else
write(*,*)"iostat",IOSTAT
end if
READ(LINP,*)TITLE
READ(LINP,'(A)')DUMMY
READ(LINP,*) NN, NE, NM, NDIM, NEN, NDN
READ(LINP,'(A)')DUMMY
READ(LINP,*) ND, NL, NCH, NPR, NMPC
write(*,*)'running'
C ----- Coordinates -----
READ(LINP,'(A)'),DUMMY
READ(LINP,*), (N,X(N),I=1,NN)
C ----- Connectivity -----
READ(LINP,'(A)'),DUMMY
READ(LINP,*), (N,NOC(N,1),NOC(N,2),MAT(N),AREA(N),DT(N),I=1,NE)
c ----- Specified Displacements -----
READ(LINP,'(A)')DUMMY
READ(LINP,*)(NU(I),U(I),I=1,ND)
C ----- Component Loads -----
READ(LINP,'(A)')DUMMY
DO 101 I=1,NN
101 F(I)=0.
READ(LINP,*) (N,F(N),I=1,NL)
C ----- Material Properties -----
READ(LINP,'(A)')DUMMY
READ(LINP,*) (N,(PM(N,J),J=1,NPR),I=1,NM)
C ----- Multi-point Constraints B1*Qi+B2*Qj=B0
IF (NMPC .GT. 0) THEN
READ(LINP,'(A)')DUMMY
DO 9 I=1,NMPC
READ(LINP,*)BT(I, 1), MPC(I, 1), BT(I, 2), MPC(I, 2), BT(I, 3)
9 CONTINUE
ENDIF
C print*,'file read'
CLOSE (LINP)
C ----- Bandwidth Evaluation -----
NBW = 0
DO 11 N=1,NE
NABS = IABS(NOC(N, 1) - NOC(N, 2)) + 1
IF (NBW .LT. NABS) NBW = NABS
11 CONTINUE
DO 13 I=1,NMPC
NABS = IABS(MPC(I, 1) - MPC(I, 2)) + 1
IF (NBW .LT. NABS) NBW = NABS
13 CONTINUE
DO 102 I=1,NN
DO 102 J=1,NBW
102 S(I,J)=0.
C ----- Stiffness Matrix -----
DO 25 N=1,NE
N1 = NOC(N, 1)
N2 = NOC(N, 2)
N3 = MAT(N)
X21 = X(N2) - X(N1)
EL = ABS(X21)
EAL = PM(N3, 1) * AREA(N) / EL
IF (NPR .GT. 1) C = PM(N3, 2)
TL = PM(N3, 1) * C * DT(N) * AREA(N) * EL / X21
C ----- Temperature Loads -----
F(N1) = F(N1) - TL
F(N2) = F(N2) + TL
C ----- Element Stiffness in Global Locations -----
S(N1, 1) = S(N1, 1) + EAL
S(N2, 1) = S(N2, 1) + EAL
IR = N1
IF (IR .GT. N2) IR = N2
IC = IABS(N2 - N1) + 1
S(IR, IC) = S(IR, IC) - EAL
25 CONTINUE
C ----- Decide Penalty Parameter CNST -----
CNST = 0
DO 27 I=1,NN
IF (CNST .LT. S(I, 1)) CNST = S(I, 1)
27 CONTINUE
CNST = CNST * 10000
C ----- Modify for Boundary Conditions -----
C --- Displacement BC ---
DO 29 I=1,ND
N = NU(I)
S(N, 1) = S(N, 1) + CNST
F(N) = F(N) + CNST * U(I)
29 CONTINUE
C --- Multi-point Constraints ---
DO 31 I=1,NMPC
I1 = MPC(I, 1)
I2 = MPC(I, 2)
S(I1, 1) = S(I1, 1) + CNST * BT(I, 1) * BT(I, 1)
S(I2, 1) = S(I2, 1) + CNST * BT(I, 2) * BT(I, 2)
IR = I1
IF (IR .GT. I2) IR = I2
IC = IABS(I2 - I1) + 1
S(IR, IC) = S(IR, IC) + CNST * BT(I, 1) * BT(I, 2)
F(I1) = F(I1) + CNST * BT(I, 1) * BT(I, 3)
F(I2) = F(I2) + CNST * BT(I, 2) * BT(I, 3)
31 CONTINUE
C ----- Equation Solving using Band Solver -----
CALL BANSOL(NN,NBW,IMAX,S,F)
WRITE(*,'(A)') TITLE
WRITE(LOUT,'(A)') TITLE
PRINT *,'NODE# DISPLACEMENT'
WRITE(LOUT,*)'NODE# DISPLACEMENT'
DO 33 I=1,NN
WRITE(*,'(I5,E15.5)') I, F(I)
WRITE(LOUT,'(I5,E15.5)') I, F(I)
33 CONTINUE
C ----- Stress Calculation -----
PRINT *, 'ELEM# STRESS'
WRITE(LOUT,*) 'ELEM# STRESS'
DO 35 N=1,NE
N1 = NOC(N, 1)
N2 = NOC(N, 2)
N3 = MAT(N)
EPS = (F(N2) - F(N1)) / (X(N2) - X(N1))
IF (NPR .GT. 1) C = PM(N3, 2)
STRESS = PM(N3, 1) * (EPS - C * DT(N))
WRITE(*,'(I5,E15.5)') N, STRESS
WRITE(LOUT,'(I5,E15.5)') N, STRESS
35 CONTINUE
C ----- Reaction Calculation -----
PRINT *, 'NODE# REACTION'
WRITE(LOUT,*) 'NODE# REACTION'
DO 37 I=1,ND
N = NU(I)
R = CNST * (U(I) - F(N))
WRITE(*,'(I5,E15.5)') N, R
WRITE(LOUT,'(I5,E15.5)') N, R
37 CONTINUE
CLOSE(LOUT)
PRINT *, 'RESULTS ARE IN FILE ', FILE2
END
SUBROUTINE BANSOL(NN,NBW,IMAX,S,F)
DIMENSION S(IMAX,1),F(1)
N = NN
C ----- Forward Elimination -----
DO 39 K=1,N-1
NBK = N - K + 1
IF ((N - K + 1) .GT. NBW) NBK = NBW
DO 43 I=K+1, NBK+K-1
I1 = I - K + 1
C = S(K, I1) / S(K, 1)
DO 41 J=I, NBK+K-1
J1 = J - I + 1
J2 = J - K + 1
S(I, J1) = S(I, J1) - C * S(K, J2)
41 CONTINUE
F(I) = F(I) - C * F(K)
43 CONTINUE
39 CONTINUE
C ----- Back Substitution -----
F(N) = F(N) / S(N, 1)
DO 47 II=1,N-1
I = N - II
NBI = N - I + 1
IF ((N - I + 1) .GT. NBW) NBI = NBW
SUM = 0.
DO 45 J=2,NBI
SUM = SUM + S(I, J) * F(I + J - 1)
45 CONTINUE
F(I) = (F(I) - SUM) / S(I, 1)
47 CONTINUE
RETURN
END
Input file
EXAMPLE 3.6
NN NE NM NDIM NEN NDN
5 2 2 1 2 1
ND NL NCH NPR NMPC
2 1 1 1 2
Node# X-Coordinate
1 0
2 0
3 -4500
4 -3000
5 0
Elem# N1 N2 Mat# Area TempRise (NCH=2 Elem Char: Area, TempRise)
1 1 3 1 1200 0
2 2 4 2 900 0
DOF# Displacement
3 0
4 0
DOF# Load
5 30000
MAT# Prop1
1 200000
2 70000
B1 i B2 j B3 (Multi-point constr. B1*Qi+B2*Qj=B3)
1 1 -.333 5 0
1 2 -.833 5 0
I have a piece of FEA code that came with an old text book I bought that I cannot get to compile. I wrote codes in Fortran 90/95 in college but I have never used Fortran 77. If someone can please help me get this to compile. I have attached the code and the input file that it reads from. The error that it is end-of-file during read. I tried a few things that I thought might work but I cannot get to work. It's probably something stupid that I am over looking.
Thanks for everybodys help
C ----------------------- FEM1D2 ------------------------
DIMENSION X(50),NOC(25,2),F(50),AREA(25),MAT(25),DT(25),
$ PM(10,2),NU(20),U(20),MPC(15,2),BT(15,3)
DIMENSION S(100,55)
C IMAX = FIRST DIMENSION OF THE S-MATRIX
CHARACTER*16 FILE1,FILE2
CHARACTER*81 DUMMY,TITLE
PRINT *, '***************************************'
PRINT *, '* PROGRAM FEM1D2 *'
PRINT *, '* WITH MULTI-POINT CONSTRAINTS *'
PRINT *, '* T.R.Chandrupatla and A.D.Belegundu *'
PRINT *, '***************************************'
IMAX=100
PRINT *,'Input Data File Name'
READ '(A)',FILE1
LINP=10
OPEN(UNIT=10,FILE=FILE1,STATUS='UNKNOWN',ACTION='READ')
IF(IOSTAT /=0)THEN
WRITE(*,*)'INPUT FILE DID NOT OPEN'
ELSE
WRITE(*,*)'IOSTAT', IOSTAT
END IF
PRINT *,'Output Data File Name'
READ(*,*),FILE2
LOUT=11
OPEN(UNIT=11,FILE=FILE2,STATUS='UNKNOWN',ACTION='WRITE')
if(IOSTAT /= 0)then
write(*,*)"OUTPUT file did not open"
else
write(*,*)"iostat",IOSTAT
end if
READ(LINP,*)TITLE
READ(LINP,'(A)')DUMMY
READ(LINP,*) NN, NE, NM, NDIM, NEN, NDN
READ(LINP,'(A)')DUMMY
READ(LINP,*) ND, NL, NCH, NPR, NMPC
write(*,*)'running'
C ----- Coordinates -----
READ(LINP,'(A)'),DUMMY
READ(LINP,*), (N,X(N),I=1,NN)
C ----- Connectivity -----
READ(LINP,'(A)'),DUMMY
READ(LINP,*), (N,NOC(N,1),NOC(N,2),MAT(N),AREA(N),DT(N),I=1,NE)
c ----- Specified Displacements -----
READ(LINP,'(A)')DUMMY
READ(LINP,*)(NU(I),U(I),I=1,ND)
C ----- Component Loads -----
READ(LINP,'(A)')DUMMY
DO 101 I=1,NN
101 F(I)=0.
READ(LINP,*) (N,F(N),I=1,NL)
C ----- Material Properties -----
READ(LINP,'(A)')DUMMY
READ(LINP,*) (N,(PM(N,J),J=1,NPR),I=1,NM)
C ----- Multi-point Constraints B1*Qi+B2*Qj=B0
IF (NMPC .GT. 0) THEN
READ(LINP,'(A)')DUMMY
DO 9 I=1,NMPC
READ(LINP,*)BT(I, 1), MPC(I, 1), BT(I, 2), MPC(I, 2), BT(I, 3)
9 CONTINUE
ENDIF
C print*,'file read'
CLOSE (LINP)
C ----- Bandwidth Evaluation -----
NBW = 0
DO 11 N=1,NE
NABS = IABS(NOC(N, 1) - NOC(N, 2)) + 1
IF (NBW .LT. NABS) NBW = NABS
11 CONTINUE
DO 13 I=1,NMPC
NABS = IABS(MPC(I, 1) - MPC(I, 2)) + 1
IF (NBW .LT. NABS) NBW = NABS
13 CONTINUE
DO 102 I=1,NN
DO 102 J=1,NBW
102 S(I,J)=0.
C ----- Stiffness Matrix -----
DO 25 N=1,NE
N1 = NOC(N, 1)
N2 = NOC(N, 2)
N3 = MAT(N)
X21 = X(N2) - X(N1)
EL = ABS(X21)
EAL = PM(N3, 1) * AREA(N) / EL
IF (NPR .GT. 1) C = PM(N3, 2)
TL = PM(N3, 1) * C * DT(N) * AREA(N) * EL / X21
C ----- Temperature Loads -----
F(N1) = F(N1) - TL
F(N2) = F(N2) + TL
C ----- Element Stiffness in Global Locations -----
S(N1, 1) = S(N1, 1) + EAL
S(N2, 1) = S(N2, 1) + EAL
IR = N1
IF (IR .GT. N2) IR = N2
IC = IABS(N2 - N1) + 1
S(IR, IC) = S(IR, IC) - EAL
25 CONTINUE
C ----- Decide Penalty Parameter CNST -----
CNST = 0
DO 27 I=1,NN
IF (CNST .LT. S(I, 1)) CNST = S(I, 1)
27 CONTINUE
CNST = CNST * 10000
C ----- Modify for Boundary Conditions -----
C --- Displacement BC ---
DO 29 I=1,ND
N = NU(I)
S(N, 1) = S(N, 1) + CNST
F(N) = F(N) + CNST * U(I)
29 CONTINUE
C --- Multi-point Constraints ---
DO 31 I=1,NMPC
I1 = MPC(I, 1)
I2 = MPC(I, 2)
S(I1, 1) = S(I1, 1) + CNST * BT(I, 1) * BT(I, 1)
S(I2, 1) = S(I2, 1) + CNST * BT(I, 2) * BT(I, 2)
IR = I1
IF (IR .GT. I2) IR = I2
IC = IABS(I2 - I1) + 1
S(IR, IC) = S(IR, IC) + CNST * BT(I, 1) * BT(I, 2)
F(I1) = F(I1) + CNST * BT(I, 1) * BT(I, 3)
F(I2) = F(I2) + CNST * BT(I, 2) * BT(I, 3)
31 CONTINUE
C ----- Equation Solving using Band Solver -----
CALL BANSOL(NN,NBW,IMAX,S,F)
WRITE(*,'(A)') TITLE
WRITE(LOUT,'(A)') TITLE
PRINT *,'NODE# DISPLACEMENT'
WRITE(LOUT,*)'NODE# DISPLACEMENT'
DO 33 I=1,NN
WRITE(*,'(I5,E15.5)') I, F(I)
WRITE(LOUT,'(I5,E15.5)') I, F(I)
33 CONTINUE
C ----- Stress Calculation -----
PRINT *, 'ELEM# STRESS'
WRITE(LOUT,*) 'ELEM# STRESS'
DO 35 N=1,NE
N1 = NOC(N, 1)
N2 = NOC(N, 2)
N3 = MAT(N)
EPS = (F(N2) - F(N1)) / (X(N2) - X(N1))
IF (NPR .GT. 1) C = PM(N3, 2)
STRESS = PM(N3, 1) * (EPS - C * DT(N))
WRITE(*,'(I5,E15.5)') N, STRESS
WRITE(LOUT,'(I5,E15.5)') N, STRESS
35 CONTINUE
C ----- Reaction Calculation -----
PRINT *, 'NODE# REACTION'
WRITE(LOUT,*) 'NODE# REACTION'
DO 37 I=1,ND
N = NU(I)
R = CNST * (U(I) - F(N))
WRITE(*,'(I5,E15.5)') N, R
WRITE(LOUT,'(I5,E15.5)') N, R
37 CONTINUE
CLOSE(LOUT)
PRINT *, 'RESULTS ARE IN FILE ', FILE2
END
SUBROUTINE BANSOL(NN,NBW,IMAX,S,F)
DIMENSION S(IMAX,1),F(1)
N = NN
C ----- Forward Elimination -----
DO 39 K=1,N-1
NBK = N - K + 1
IF ((N - K + 1) .GT. NBW) NBK = NBW
DO 43 I=K+1, NBK+K-1
I1 = I - K + 1
C = S(K, I1) / S(K, 1)
DO 41 J=I, NBK+K-1
J1 = J - I + 1
J2 = J - K + 1
S(I, J1) = S(I, J1) - C * S(K, J2)
41 CONTINUE
F(I) = F(I) - C * F(K)
43 CONTINUE
39 CONTINUE
C ----- Back Substitution -----
F(N) = F(N) / S(N, 1)
DO 47 II=1,N-1
I = N - II
NBI = N - I + 1
IF ((N - I + 1) .GT. NBW) NBI = NBW
SUM = 0.
DO 45 J=2,NBI
SUM = SUM + S(I, J) * F(I + J - 1)
45 CONTINUE
F(I) = (F(I) - SUM) / S(I, 1)
47 CONTINUE
RETURN
END
Input file
EXAMPLE 3.6
NN NE NM NDIM NEN NDN
5 2 2 1 2 1
ND NL NCH NPR NMPC
2 1 1 1 2
Node# X-Coordinate
1 0
2 0
3 -4500
4 -3000
5 0
Elem# N1 N2 Mat# Area TempRise (NCH=2 Elem Char: Area, TempRise)
1 1 3 1 1200 0
2 2 4 2 900 0
DOF# Displacement
3 0
4 0
DOF# Load
5 30000
MAT# Prop1
1 200000
2 70000
B1 i B2 j B3 (Multi-point constr. B1*Qi+B2*Qj=B3)
1 1 -.333 5 0
1 2 -.833 5 0
