Years ago I wrote a program in Fortran for a structural analysis class. At the time I had very little programming experience and mottled my way through the assignment. Fast forward to now and I'm attempting to teach myself Python 3. I figured I'd would try to convert the old program from Fortran to Python. I've made some good progress but I've hit a wall. I've got the program running to the point where all the element stiffness matrices are defined and combined into the structure stiffness matrix using the code number technique. The structure stiffness matrix is triangular. At this point in the Fortran code, it calls a subroutine which I believe solves for d in [K]d = p using Gauss Elimination. Here's the Fortran code:
DIMENSION BGK(100,30),XL(60),COSX(60),COSY(60),COSZ(60),SK(6,6),
1AR(60),NC(60,6),QL(100),P(60)
OPEN(UNIT = 2, FILE = '3DTRUSS.OUT', STATUS = 'UNKNOWN')
PRINT 4
4 FORMAT(///" ***2D OR 3D TRUSS ANALYSIS***"/,
1/" INPUT IS IN FREE FORM."/)
PRINT 8
8 FORMAT(" TYPE THE NUMBER OF MEMBERS, THE NUMBER OF DEGREES ",
1"OF FREEDOM."/" AND THE VALUE OF THE MODULUS E.")
READ*,NM,ND,E
PRINT 12
12 FORMAT(" TYPE THE MEMBER END COORDINATES, MEMBER 1, 2, ETC.,"/
1" IN THE ORDER:"/" XI YI ZI XJ YJ ZJ")
DO 16 I=1,NM
READ*,XI,YI,ZI,XJ,YJ,ZJ
XL(I)=SQRT((ZJ-ZI)**2+(YJ-YI)**2+(XJ-XI)**2)
COSX(I)=(XJ-XI)/XL(I)
COSY(I)=(YJ-YI)/XL(I)
16 COSZ(I)=(ZJ-ZI)/XL(I)
PRINT 20
20 FORMAT(" TYPE THE MEMBER AREAS - A(1), A(2), ETC.")
READ*,(AR(I),I=1,NM)
PRINT 24
24 FORMAT(" TYPE THE LOCATION VECTORS; MEMBER 1, 2, ETC.")
DO 28 I=1,NM
28 READ*,(NC(I,J),J=1,6)
PRINT 32
32 FORMAT(" TYPE THE EXTERNAL LOADS(Q).")
READ*,(QL(I),I=1,ND)
WRITE (2,36)
36 FORMAT(//" INPUT VALUES:"//" MEMBER",30X,"LOCATION VECTORS"/,
1" NUMBER LENGTH AREA NUMBERS"/)
DO 40 I=1,NM
40 WRITE (2,44) I,XL(I),AR(I),(NC(I,J),J=1,6)
44 FORMAT(2X,I4,2E13.5,1X,6I3)
WRITE (2,48) E
48 FORMAT(/" THE VALUE OF E IS ",E11.5)
WRITE (2,52)
52 FORMAT(1X/" ID NO. LOADS APPLIED "/)
DO 56 I=1,ND
56 WRITE (2,60) I,QL(I)
60 FORMAT(2X,I3,4X,E13.5)
IBND=0
DO 64 I=1,NM
DO 64 J=1,5
JJ=J+1
DO 64 K=JJ,6
IF(NC(I,J).EQ.0.OR.NC(I,K).EQ.0)GO TO 64
M=IABS(NC(I,J)-NC(I,K))
IF(M.GT.IBND)IBND=M
64 CONTINUE
IBND=IBND+1
WRITE (2,68) IBND
68 FORMAT(/" THE BANDWIDTH IS ",I2)
IF(IBND.LE.30)GO TO 72
WRITE (2,70)
70 FORMAT(/" BANDWIDTH EXCEEDS DIMENSION."/)
STOP
72 DO 76 I=1,ND
DO 76 J=1,IBND
76 BGK(I,J)=0.
DO 84 N=1,NM
SK(1,1)=COSX(N)**2
SK(4,4)=COSX(N)**2
SK(1,2)=COSX(N)*COSY(N)
SK(4,5)=COSX(N)*COSY(N)
SK(1,3)=COSX(N)*COSZ(N)
SK(4,6)=COSX(N)*COSZ(N)
SK(1,4)=-SK(1,1)
SK(1,5)=-SK(1,2)
SK(2,4)=-SK(1,2)
SK(1,6)=-SK(1,3)
SK(3,4)=-SK(1,3)
SK(2,2)=COSY(N)**2
SK(5,5)=COSY(N)**2
SK(2,3)=COSY(N)*COSZ(N)
SK(5,6)=COSY(N)*COSZ(N)
SK(2,5)=-SK(2,2)
SK(2,6)=-SK(2,3)
SK(3,5)=-SK(2,3)
SK(3,3)=COSZ(N)**2
SK(6,6)=COSZ(N)**2
SK(3,6)=-SK(3,3)
DO 80 I=1,6
DO 80 J=I,6
80 SK(I,J)=(E*AR(N)/XL(N))*SK(I,J)
DO 84 I=1,6
DO 84 J=I,6
K=NC(N,I)
L=NC(N,J)
IF(K.EQ.0.OR.L.EQ.0)GO TO 84
IF(K.LE.L)GO TO 82
IT=K
K=L
L=IT
82 IPOS=L-K+1
BGK(K,IPOS)=BGK(K,IPOS)+SK(I,J)
84 CONTINUE
CALL BNDSOL(BGK,QL,ND,IBND)
DO 88 N=1,NM
C=E*AR(N)/XL(N)
SK(1,4)=COSX(N)
SK(1,1)=-SK(1,4)
SK(1,5)=COSY(N)
SK(1,2)=-SK(1,5)
SK(1,6)=COSZ(N)
SK(1,3)=-SK(1,6)
P(N)=0.
DO 88 J=1,6
K=NC(N,J)
IF(K.EQ.0)GO TO 88
P(N)=P(N)+C*SK(1,J)*QL(K)
88 CONTINUE
WRITE (2,92)
92 FORMAT(//" OUTPUT VALUES "//" ID NO. DEFLECTIONS",3X,
1"AXIAL FORCE"/)
IF(NM-ND)96,106,106
96 DO 98 I=1,NM
98 WRITE (2,100) I,QL(I),P(I)
100 FORMAT(3X,I3,2X,2E14.5)
J=NM+1
DO 102 I=J,ND
102 WRITE (2,104) I,QL(I)
104 FORMAT(3X,I3,2X,E14.5)
GO TO 120
106 DO 108 I=1,ND
108 WRITE (2,100) I,QL(I),P(I)
IF(NM.EQ.ND)GO TO 120
J=ND+1
DO 110 I=J,NM
110 WRITE (2,112) I,P(I)
112 FORMAT(3X,I3,16X,E14.5)
120 WRITE (2,124)
124 FORMAT(/" ANALYSIS COMPLETE."///)
STOP
END
SUBROUTINE BNDSOL(BGK,Q,NDIS,MB)
DIMENSION BGK(100,30),Q(100),F(30)
N=0
500 N=N+1
Q(N)=Q(N)/BGK(N,1)
IF(N-NDIS)550,700,550
550 DO 600 K=2,MB
F(K)=BGK(N,K)
600 BGK(N,K)=BGK(N,K)/BGK(N,1)
DO 660 L=2,MB
I=N+L-1
IF(NDIS-I)660,640,640
640 J=0
DO 650 K=L,MB
J=J+1
650 BGK(I,J)=BGK(I,J)-F(L)*BGK(N,K)
Q(I)=Q(I)-F(L)*Q(N)
660 CONTINUE
GO TO 500
700 N=N-1
IF(N)750,900,750
750 DO 800 K=2,MB
L=N+K-1
IF(NDIS-L)800,770,770
770 Q(N)=Q(N)-BGK(N,K)*Q(L)
800 CONTINUE
GO TO 700
900 RETURN
END
I make it the continue statement after line 84 and my head fogs up. I would appreciate it immensely if someone could help me understand the rest of the code so I can convert it to python. I've attached my python code and am open to any suggestions on how I could make it better. This is my first attempt at writing in python.
import math
import os
from numpy import array, zeros
from numpy.linalg import solve
def buildMember(temp,coords):
mbr = []
mno = temp[0]
jti = temp[1]
jtj = temp[2]
mbr.append(mno)
for row in range(len(coords)):
if jti == coords[row][0]:
xyzi = coords[row][1:]
mbr.append(xyzi)
for row in range(len(coords)):
if jtj == coords[row][0]:
xyzj = coords[row][1:]
mbr.append(xyzj)
return mbr
def codeMatrix(temp,table):
codeno = []
n = temp[0]
pt1 = temp[1]
pt2 = temp[2]
codeno.append(n)
for row in range(len(table)):
if pt1 == table[row][0]:
x1 = table[row][1]
y1 = table[row][2]
z1 = table[row][3]
codeno.append(x1)
codeno.append(y1)
codeno.append(z1)
for row in range(len(table)):
if pt2 == table[row][0]:
x2 = table[row][1]
y2 = table[row][2]
z2 = table[row][3]
codeno.append(x2)
codeno.append(y2)
codeno.append(z2)
B = []
for i in range(1,len(codeno)):
if codeno[i] != 0:
B.append(codeno[i])
band = max(B)-min(B)+1
return codeno, band
os.chdir(r'c:/users/colby/documents/dropbox/dropbox/3dtruss')
#os.chdir('c:/pyzo2013c/projects')
print(" ***2D or 3D Truss Analysis*** ")
infile = open(input("Enter input file name; name.txt:"),'r')
output = open(input("Enter output file name; name.txt:"), "w+")
tempfile = open('tempfile.txt','w+')
NM, DOF, E = eval(infile.readline())
#NM, DOF, E = eval(input("Enter the number of members, the number of degrees"+\
#"of freedom and \nthe value of the modulus E. Separate with commas:"))
# Read in joint coordinates and create point table
#jts = eval(input("Enter number of joints:"))
jts = eval(infile.readline())
output.write('{0:^8}{1:^6}{2:^6}{3:^6}\n'.format('Joint','X','Y','Z'))
points =[]
for i in range(jts):
#tempjt = eval(input("Enter joint {} coordinates; [jt no.,x,y,z]:".format(i+1)))
tempjt = eval(infile.readline())
points.append(tempjt)
output.write('{0:^8d}{1:^6.1f}{2:^6.1f}{3:^6.1f}\n'.format(tempjt[0],tempjt[1],tempjt[2],tempjt[3]))
print(points)
dtable = []
for i in range(jts):
#tempd = eval(input("Enter the code-number for point {}; [pt no,x,y,z]:".format(i+1)))
tempd = eval(infile.readline())
dtable.append(tempd)
print(dtable)
# Assign joints to members
members = []
area = []
codes = []
xl=[]
cosx=[]
cosy=[]
cosz=[]
bndwth = []
output.write('\n{0:^8} {1:^9} {2:^8} {3:^18}\n'.format('MEMBER','LENGTH','AREA','CODE NUMBERS'))
for i in range(1,NM+1):
#tempmbr = eval(input("Enter member {} data; [no, i, j] where i & j are pt numbers:".format(i)))
tempmbr = eval(infile.readline())
m = buildMember(tempmbr, points)
members.append(m)
codenumbers, bw = codeMatrix(tempmbr,dtable)
codes.append(codenumbers)
xlen = math.sqrt((m[2][2]-m[1][2])**2+(m[2][1]-m[1][1])**2+(m[2][0]-m[1][0])**2)
cos_x = (m[2][0]-m[1][0])/xlen
cos_y = (m[2][1]-m[1][1])/xlen
cos_z = (m[2][2]-m[1][2])/xlen
#a = eval(input("Enter area of member {}:".format(i)))
a = eval(infile.readline())
area.append(a)
xl.append(xlen)
cosx.append(cos_x)
cosy.append(cos_y)
cosz.append(cos_z)
print('{0:8.4f},{1:7.3f},{2:7.3f},{3:7.3f}'.format(xlen,cos_x,cos_y,cos_z))
print('Bandwidth is {}\n'.format(bw))
#code_str = str(codenumbers[1:])
output.write('{0:^8} {1:>8.4f} {2:>8.4f} {3:^}\n'.format(i,xlen,a,str(codenumbers[1:])))
bndwth.append(bw)
print(members)
print(codes)
#q = eval(input("Enter the external loads. Loads must be ordered per code numbering:"))
q = eval(infile.readline())
Q = array(q)
print(q)
output.write('The value of modulus E is {:10.1f}\n\n'.format(E))
output.write('Applied loads corresponding to structure node labels:\n{}\n\n'.format(q))
B = max(bndwth) # defines bandwith
print('Bandwidth is {}\n'.format(B))
output.write('The bandwidth is {}\n\n'.format(B))
k = zeros((6,6))
K = zeros((DOF,B))
print(k)
print(K)
for n in range(NM):
k[0,0]=cosx[n]**2
k[3,3]=cosx[n]**2
k[0,1]=cosx[n]*cosy[n]
k[3,4]=cosx[n]*cosy[n]
k[0,2]=cosx[n]*cosz[n]
k[3,5]=cosx[n]*cosz[n]
k[0,3]=-k[0,0]
k[0,4]=-k[0,1]
k[1,3]=-k[0,1]
k[0,5]=-k[0,2]
k[2,3]=-k[0,2]
k[1,1]=cosy[n]**2
k[4,4]=cosy[n]**2
k[1,2]=cosy[n]*cosz[n]
k[4,5]=cosy[n]*cosz[n]
k[1,4]=-k[1,1]
k[1,5]=-k[1,2]
k[2,4]=-k[1,2]
k[2,2]=cosz[n]**2
k[5,5]=cosz[n]**2
k[2,5]=-k[2,2]
for i in range(6):
for j in range(6):
k[i,j] = (E*area[n]/xl[n])*k[i,j]
print('k[{},{}] = {:7.3f}'.format(i+1,j+1,k[i,j]))
print('n i j m l K k',end='\n')
for i in range(DOF):
for j in range(B):
m = codes[n][i+1]
l = codes[n][j+1]
if m != 0:
if l != 0:
K[m-1,l-1] = K[m-1,l-1]+k[i,j]
print(n,i,j,m,l,K[m-1,l-1],k[i,j],end='\n')
tempfile.write('{} {} {} {} {} {} {}\n'.format(n+1,i+1,j+1,m,l,K[m-1,l-1],k[i,j]))
print(K)
output.write('Global Stiffness Matrix\n{}'.format(K))
D = solve(K,Q)
print(D)
infile.close()
output.close()
tempfile.close()
Thanks
