py_fem_truss.py
cat << EOT > inp_test1.txt
5 7 1 2 2
1.0 1.0 0.0 0.0 0.0 0.0
1 2 1
1 3 1
2 3 1
2 4 1
3 4 1
3 5 1
4 5 1
0.0 0.0000 0.0
4.0 0.0000 0.0
4.0 2.3094 0.0
8.0 0.0000 0.0
8.0 4.6188 0.0
4 1 0 0.0 0.0
5 1 1 0.0 0.0
1 0.0 -10.0
3 0.0 -10.0
EOT
python3 py_fem_truss.py inp_test1.txt out_test1.txt
| ################################### | |
| # 2D Truss Analysis | |
| ################################### | |
| import numpy as np | |
| import sys | |
| import time | |
| def PRINP_TRS(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node): | |
| fout=open(fnameW,'w') | |
| # print out of input data | |
| print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout) | |
| print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout) | |
| print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}' | |
| .format('sec','E','A','alpha','gamma','gkh','gkv'),file=fout) | |
| for i in range(0,nsec): | |
| print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}' | |
| .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i]),file=fout) | |
| print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s}' | |
| .format('node','x','y','fx','fy','deltaT','kox','koy'),file=fout) | |
| for i in range(0,npoin): | |
| lp=i+1 | |
| print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d}' | |
| .format(lp,x[0,i],x[1,i],fp[2*lp-2],fp[2*lp-1],deltaT[i],mpfix[0,i],mpfix[1,i]),file=fout) | |
| print('{0:>5s} {1:>5s} {2:>5s} {3:>15s} {4:>15s}' | |
| .format('node','kox','koy','rdis_x','rdis_y'),file=fout) | |
| for i in range(0,npoin): | |
| if 0<mpfix[0,i]+mpfix[1,i]: | |
| lp=i+1 | |
| print('{0:5d} {1:5d} {2:5d} {3:15.7e} {4:15.7e}' | |
| .format(lp,mpfix[0,i],mpfix[1,i],rdis[0,i],rdis[1,i]),file=fout) | |
| print('{0:>5s} {1:>5s} {2:>5s} {3:>5s}'.format('elem','i','j','sec'),file=fout) | |
| for ne in range(0,nele): | |
| print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout) | |
| fout.close() | |
| def PROUT_TRS(fnameW,npoin,nele,disg,fsec): | |
| fout=open(fnameW,'a') | |
| # displacement | |
| print('{0:>5s} {1:>15s} {2:>15s}'.format('node','dis-x','dis-y'),file=fout) | |
| for i in range(0,npoin): | |
| lp=i+1 | |
| print('{0:5d} {1:15.7e} {2:15.7e}'.format(lp,disg[2*lp-2],disg[2*lp-1]),file=fout) | |
| # section force | |
| print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s}' | |
| .format('elem','N_i','S_i','N_j','S_j'),file=fout) | |
| for ne in range(0,nele): | |
| print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e}' | |
| .format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne]),file=fout) | |
| fout.close() | |
| def SM_TRS(ne,node,x,ae): | |
| ek=np.zeros([4,4],dtype=np.float64) # local stiffness matrix | |
| tt=np.zeros([4,4],dtype=np.float64) # transformation matrix | |
| i=node[0,ne]-1 | |
| j=node[1,ne]-1 | |
| x1=x[0,i] | |
| y1=x[1,i] | |
| x2=x[0,j] | |
| y2=x[1,j] | |
| xx=x2-x1 | |
| yy=y2-y1 | |
| el=np.sqrt(xx**2+yy**2) | |
| m=node[2,ne]-1 | |
| ee=ae[0,m] | |
| aa=ae[1,m] | |
| EA=ee*aa | |
| ek[0,0]= EA/el; ek[0,1]=0.0; ek[0,2]=-EA/el; ek[0,3]=0.0 | |
| ek[1,0]= 0.0; ek[1,1]=0.0; ek[1,2]= 0.0; ek[1,3]=0.0 | |
| ek[2,0]=-EA/el; ek[2,1]=0.0; ek[2,2]= EA/el; ek[2,3]=0.0 | |
| ek[3,0]= 0.0; ek[3,1]=0.0; ek[3,2]= 0.0; ek[3,3]=0.0 | |
| s=yy/el | |
| c=xx/el | |
| tt[0,0]= c; tt[0,1]= s; tt[0,2]=0.0; tt[0,3]=0.0 | |
| tt[1,0]= -s; tt[1,1]= c; tt[1,2]=0.0; tt[1,3]=0.0 | |
| tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]= c; tt[2,3]= s | |
| tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]= -s; tt[3,3]= c | |
| return ek,tt | |
| def BFVEC_TRS(ne,node,x,ae): | |
| # Body force vector in global coordinate system | |
| bfe_g=np.zeros(4,dtype=np.float64) | |
| i=node[0,ne]-1 | |
| j=node[1,ne]-1 | |
| m=node[2,ne]-1 | |
| AA =ae[1,m] | |
| gamma=ae[3,m] | |
| gkh =ae[4,m] | |
| gkv =ae[5,m] | |
| el=np.sqrt((x[0,j]-x[0,i])**2+(x[1,j]-x[1,i])**2) | |
| bfe_g[0]=0.5*gamma*AA*el*gkh | |
| bfe_g[1]=0.5*gamma*AA*el*gkv | |
| bfe_g[2]=0.5*gamma*AA*el*gkh | |
| bfe_g[3]=0.5*gamma*AA*el*gkv | |
| return bfe_g | |
| def TFVEC_TRS(ne,node,x,ae,deltaT): | |
| # Thermal load vector in local coordinate system | |
| tfe_l=np.zeros(4,dtype=np.float64) | |
| i=node[0,ne]-1 | |
| j=node[1,ne]-1 | |
| m=node[2,ne]-1 | |
| E =ae[0,m] | |
| AA =ae[1,m] | |
| alpha=ae[2,m] | |
| tempe=0.5*(deltaT[i]+deltaT[j]) | |
| tfe_l[0]=-E*AA*alpha*tempe | |
| tfe_l[1]=0.0 | |
| tfe_l[2]= E*AA*alpha*tempe | |
| tfe_l[3]=0.0 | |
| return tfe_l | |
| # Main routine | |
| start=time.time() | |
| args = sys.argv | |
| fnameR=args[1] # input data file | |
| fnameW=args[2] # output data file | |
| f=open(fnameR,'r') | |
| text=f.readline() | |
| text=text.strip() | |
| text=text.split() | |
| npoin=int(text[0]) # Number of nodes | |
| nele =int(text[1]) # Number of elements | |
| nsec =int(text[2]) # Number of sections | |
| npfix=int(text[3]) # Number of restricted nodes | |
| nlod =int(text[4]) # Number of loaded nodes | |
| nod=2 # Number of nodes per element | |
| nfree=2 # Degree of freedom per node | |
| n=nfree*npoin | |
| x =np.zeros([2,npoin],dtype=np.float64) # Coordinates of nodes | |
| deltaT=np.zeros(npoin,dtype=np.float64) # Temperature change of nodes | |
| ae =np.zeros([6,nsec],dtype=np.float64) # Section characteristics | |
| node =np.zeros([nod+1,nele],dtype=np.int) # Node-element relationship | |
| fp =np.zeros(nfree*npoin,dtype=np.float64) # External force vector | |
| mpfix =np.zeros([nfree,npoin],dtype=np.int) # Ristrict conditions | |
| rdis =np.zeros([nfree,npoin],dtype=np.float64) # Ristricted displacement | |
| ir =np.zeros(nod*nfree,dtype=np.int) # Work vector for matrix assembly | |
| gk =np.zeros([n,n],dtype=np.float64) # Global stiffness matrix | |
| fsec =np.zeros([nod*nfree,nele],dtype=np.float64) # Section force vector | |
| work =np.zeros(nod*nfree,dtype=np.float64) # work vector for section force calculation | |
| # section characteristics | |
| for i in range(0,nsec): | |
| text=f.readline() | |
| text=text.strip() | |
| text=text.split() | |
| ae[0,i]=float(text[0]) # E : Elastic modulus | |
| ae[1,i]=float(text[1]) # AA : Section area | |
| ae[2,i]=float(text[2]) # alpha : Thermal expansion coefficient | |
| ae[3,i]=float(text[3]) # gamma : Unit weight | |
| ae[4,i]=float(text[4]) # gkh : Acceleration in horizontal direction | |
| ae[5,i]=float(text[5]) # gkv : Acceleration in vertical direction | |
| # element-node | |
| for i in range(0,nele): | |
| text=f.readline() | |
| text=text.strip() | |
| text=text.split() | |
| node[0,i]=int(text[0]) #node_1 | |
| node[1,i]=int(text[1]) #node_2 | |
| node[2,i]=int(text[2]) #section characteristic number | |
| # node coordinates | |
| for i in range(0,npoin): | |
| text=f.readline() | |
| text=text.strip() | |
| text=text.split() | |
| x[0,i]=float(text[0]) # x-coordinate | |
| x[1,i]=float(text[1]) # y-coordinate | |
| deltaT[i]=float(text[2]) # Temperature change | |
| # boundary conditions (0:free, 1:restricted) | |
| for i in range(0,npfix): | |
| text=f.readline() | |
| text=text.strip() | |
| text=text.split() | |
| lp=int(text[0]) #fixed node | |
| mpfix[0,lp-1]=int(text[1]) #fixed in x-direction | |
| mpfix[1,lp-1]=int(text[2]) #fixed in y-direction | |
| rdis[0,lp-1]=float(text[3]) #fixed displacement in x-direction | |
| rdis[1,lp-1]=float(text[4]) #fixed displacement in y-direction | |
| # load | |
| if 0<nlod: | |
| for i in range(0,nlod): | |
| text=f.readline() | |
| text=text.strip() | |
| text=text.split() | |
| lp=int(text[0]) #loaded node | |
| fp[2*lp-2]=float(text[1]) #load in x-direction | |
| fp[2*lp-1]=float(text[2]) #load in y-direction | |
| f.close() | |
| PRINP_TRS(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node) | |
| # assembly of stiffness matrix & load vectors | |
| for ne in range(0,nele): | |
| ek,tt=SM_TRS(ne,node,x,ae) # Stiffness matrix in local coordinate | |
| ck =np.dot(np.dot(tt.T,ek),tt) # Stiffness matrix in global coordinate | |
| tfe_l=TFVEC_TRS(ne,node,x,ae,deltaT) # Thermal load vector in local coordinate | |
| tfe =np.dot(tt.T,tfe_l) # Thermal load vector in global coordinate | |
| bfe =BFVEC_TRS(ne,node,x,ae) # Body force vector in global coordinate | |
| i=node[0,ne]-1 | |
| j=node[1,ne]-1 | |
| ir[3]=2*j+1; ir[2]=ir[3]-1 | |
| ir[1]=2*i+1; ir[0]=ir[1]-1 | |
| for i in range(0,nod*nfree): | |
| it=ir[i] | |
| fp[it]=fp[it]+tfe[i]+bfe[i] | |
| for j in range(0,nod*nfree): | |
| jt=ir[j] | |
| gk[it,jt]=gk[it,jt]+ck[i,j] | |
| # treatment of boundary conditions | |
| for i in range(0,npoin): | |
| for j in range(0,nfree): | |
| if mpfix[j,i]==1: | |
| iz=i*nfree+j | |
| fp[iz]=0.0 | |
| for k in range(0,n): | |
| fp[k]=fp[k]-rdis[j,i]*gk[k,iz] | |
| gk[k,iz]=0.0 | |
| gk[iz,iz]=1.0 | |
| # solution of simultaneous linear equations | |
| disg = np.linalg.solve(gk, fp) | |
| # recovery of restricted displacements | |
| for i in range(0,npoin): | |
| for j in range(0,nfree): | |
| if mpfix[j,i]==1: | |
| iz=i*nfree+j | |
| for k in range(0,n): | |
| disg[iz]=rdis[j,i] | |
| # calculation of section force | |
| for ne in range(0,nele): | |
| ek,tt=SM_TRS(ne,node,x,ae) | |
| i=node[0,ne]-1 | |
| j=node[1,ne]-1 | |
| m=node[2,ne]-1 | |
| E =ae[0,m] | |
| AA =ae[1,m] | |
| alpha=ae[2,m] | |
| tempe=0.5*(deltaT[i]+deltaT[j]) | |
| work[0]=disg[2*i]; work[1]=disg[2*i+1]; | |
| work[2]=disg[2*j]; work[3]=disg[2*j+1]; | |
| fsec[:,ne]=np.dot(ek,np.dot(tt,work)) | |
| fsec[0,ne]=fsec[0,ne]+E*AA*alpha*tempe | |
| fsec[2,ne]=fsec[2,ne]-E*AA*alpha*tempe | |
| PROUT_TRS(fnameW,npoin,nele,disg,fsec) | |
| dtime=time.time()-start | |
| print('n={0} time={1:.3f}'.format(n,dtime)+' sec') | |
| fout=open(fnameW,'a') | |
| print('n={0} time={1:.3f}'.format(n,dtime)+' sec',file=fout) | |
| fout.close() |