Skip to content

Instantly share code, notes, and snippets.

@damyarou

damyarou/0_python_FEM_gfrm.md

Last active December 27, 2016 03:41
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save damyarou/bce74ff69c690ecdd25449feb4abf3e3 to your computer and use it in GitHub Desktop.
Save damyarou/bce74ff69c690ecdd25449feb4abf3e3 to your computer and use it in GitHub Desktop.
Download ZIP
Geometrically Nonlinear Analysis of 2D frame
Raw
0_python_FEM_gfrm.md

Uploaded codes

  • py_fem_gfrmAL.py
  • inp_gfrm_canti.txt
  • inp_gfrm_cable.txt
  • inp_gfrm_arch.txt
  • inp_gfrm_lee.txt
  • py_fig_gfrm.py
  • py_fig_gfrm_mode.py
Raw
1_py_fem_gfrmAL.py
##############################################
# 2D Frame Geometrical Non-linear Analysis
##############################################
import numpy as np
import sys
import time
def PRINP_FRM(fnameW,npoin,nele,nsec,npfix,nlod,nnmax,ae,x,df,mpfix,node):
fout=open(fnameW,'w')
# print out of input data
print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'.format('npoin','nele','nsec','npfix','nlod','nnmax'),file=fout)
print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}'.format(npoin,nele,nsec,npfix,nlod,nnmax),file=fout)
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'
.format('sec','E','A','I'),file=fout)
for i in range(0,nsec):
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'
.format(i+1,ae[0,i],ae[1,i],ae[2,i]),file=fout)
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s} {8:>5s}'
.format('node','x','y','fx','fy','fr','kox','koy','kor'),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} {8:5d}'
.format(lp,x[0,i],x[1,i],df[3*lp-3],df[3*lp-2],df[3*lp-1],mpfix[0,i],mpfix[1,i],mpfix[2,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_FRM(fnameW,nnn,iii,lam,npoin,nele,fp,dist,fsec,dr):
fout=open(fnameW,'a')
print('* nnn={0:} iii={1:} lam={2:15.7e}'.format(nnn,iii,lam),file=fout)
# displacement
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s} {8:>15s} {9:>15s}'.format('node','fp-x','fp-y','fp-r','dis-x','dis-y','dis-r','dr-x','dr-y','dr-r'),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:15.7e} {7:15.7e} {8:15.7e} {9:15.7e}'.format(lp,fp[3*lp-3],fp[3*lp-2],fp[3*lp-1],dist[3*lp-3],dist[3*lp-2],dist[3*lp-1],dr[3*lp-3],dr[3*lp-2],dr[3*lp-1]),file=fout)
# section force
print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
.format('elem','N_i','S_i','M_i','N_j','S_j','M_j'),file=fout)
for ne in range(0,nele):
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
.format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
fout.close()
def SM_GFRM(ne,node,ae,x1,y1,x2,y2):
ek=np.zeros([6,6],dtype=np.float64) # local stiffness matrix
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]
ai=ae[2,m]
EA=ee*aa
EI=ee*ai
ek[0,0]= EA/el; ek[0,1]= 0.0; ek[0,2]= 0.0; ek[0,3]=-EA/el; ek[0,4]= 0.0; ek[0,5]= 0.0
ek[1,0]= 0.0; ek[1,1]= 12*EI/el**3; ek[1,2]= 6*EI/el**2; ek[1,3]= 0.0; ek[1,4]=-12*EI/el**3; ek[1,5]= 6*EI/el**2
ek[2,0]= 0.0; ek[2,1]= 6*EI/el**2; ek[2,2]= 4*EI/el ; ek[2,3]= 0.0; ek[2,4]= -6*EI/el**2; ek[2,5]= 2*EI/el
ek[3,0]=-EA/el; ek[3,1]= 0.0; ek[3,2]= 0.0; ek[3,3]= EA/el; ek[3,4]= 0.0; ek[3,5]= 0.0
ek[4,0]= 0.0; ek[4,1]=-12*EI/el**3; ek[4,2]=-6*EI/el**2; ek[4,3]= 0.0; ek[4,4]= 12*EI/el**3; ek[4,5]=-6*EI/el**2
ek[5,0]= 0.0; ek[5,1]= 6*EI/el**2; ek[5,2]= 2*EI/el ; ek[5,3]= 0.0; ek[5,4]= -6*EI/el**2; ek[5,5]= 4*EI/el
return ek
def SG_GFRM(x1,y1,x2,y2):
eg=np.zeros([6,6],dtype=np.float64) # transformation matrix
xx=x2-x1
yy=y2-y1
el=np.sqrt(xx**2+yy**2)
eg[0,0]= 1/el; eg[0,1]= 0.0; eg[0,2]= 0.0; eg[0,3]=-1/el; eg[0,4]= 0.0; eg[0,5]= 0.0
eg[1,0]= 0.0; eg[1,1]= 6/5/el; eg[1,2]= 1/10; eg[1,3]= 0.0; eg[1,4]=-6/5/el; eg[1,5]= 1/10
eg[2,0]= 0.0; eg[2,1]= 1/10; eg[2,2]=2*el/15; eg[2,3]= 0.0; eg[2,4]= -1/10; eg[2,5]= -el/30
eg[3,0]=-1/el; eg[3,1]= 0.0; eg[3,2]= 0.0; eg[3,3]= 1/el; eg[3,4]= 0.0; eg[3,5]= 0.0
eg[4,0]= 0.0; eg[4,1]=-6/5/el; eg[4,2]= -1/10; eg[4,3]= 0.0; eg[4,4]= 6/5/el; eg[4,5]= -1/10
eg[5,0]= 0.0; eg[5,1]= 1/10; eg[5,2]= -el/30; eg[5,3]= 0.0; eg[5,4]= -1/10; eg[5,5]=2*el/15
return eg
def TM_GFRM(x1,y1,x2,y2):
tt=np.zeros([6,6],dtype=np.float64) # transformation matrix
xx=x2-x1
yy=y2-y1
el=np.sqrt(xx**2+yy**2)
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[0,4]=0.0; tt[0,5]=0.0
tt[1,0]= -s; tt[1,1]= c; tt[1,2]=0.0; tt[1,3]=0.0; tt[1,4]=0.0; tt[1,5]=0.0
tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]=1.0; tt[2,3]=0.0; tt[2,4]=0.0; tt[2,5]=0.0
tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]=0.0; tt[3,3]= c; tt[3,4]= s; tt[3,5]=0.0
tt[4,0]=0.0; tt[4,1]=0.0; tt[4,2]=0.0; tt[4,3]= -s; tt[4,4]= c; tt[4,5]=0.0
tt[5,0]=0.0; tt[5,1]=0.0; tt[5,2]=0.0; tt[5,3]=0.0; tt[5,4]=0.0; tt[5,5]=1.0
return tt
# Main routine
start=time.time()
args = sys.argv
fnameR=args[1] # Input data file
fnameW=args[2] # Output data file
nnmax=int(args[3]) # Number of loading step
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=3 # Degree of freedom per node
n=nfree*npoin
x =np.zeros([2,npoin],dtype=np.float64) # Coordinates of nodes
ae =np.zeros([3,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
df =np.zeros(nfree*npoin,dtype=np.float64) # Increment of External force vector
mpfix =np.zeros([nfree,npoin],dtype=np.int) # Ristrict conditions
ir =np.zeros(nod*nfree,dtype=np.int) # Work vector for matrix assembly
dist =np.zeros(n,dtype=np.float64) # Total displacement in global coordinate system
fsec =np.zeros([nod*nfree,nele],dtype=np.float64) # Section force vector
dr =np.zeros(n,dtype=np.float64) # Un-balanced force in global coordinate system
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]) # AI : Moment of inertia
# 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
# 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
mpfix[2,lp-1]=int(text[3]) #fixed in rotation
# 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
df[3*lp-3]=float(text[1]) #load increment in x-direction
df[3*lp-2]=float(text[2]) #load increment in y-direction
df[3*lp-1]=float(text[3]) #load increment n rotation
f.close()
PRINP_FRM(fnameW,npoin,nele,nsec,npfix,nlod,nnmax,ae,x,df,mpfix,node)
itmax=50
alpha=1.0
spara=np.sqrt(alpha/np.sum(df*df))
dis_ref=np.zeros(n,dtype=np.float64)
PROUT_FRM(fnameW,0,0,0.0,npoin,nele,fp,dist,fsec,dr)
#===========================
# Loading step
for nnn in range(1,nnmax):
#===========================
lam=0.0
gk=np.zeros([n,n],dtype=np.float64) # Global stiffness matrix
for ne in range(0,nele):
i=node[0,ne]-1
j=node[1,ne]-1
x1=x[0,i]+dist[3*i]
y1=x[1,i]+dist[3*i+1]
x2=x[0,j]+dist[3*j]
y2=x[1,j]+dist[3*j+1]
P=0.5*(fsec[3,ne]-fsec[0,ne]) # Axial force for non-linear part of stiffness matrix
ek=SM_GFRM(ne,node,ae,x1,y1,x2,y2) # Linear part of Stiffness matrix
eg=P*SG_GFRM(x1,y1,x2,y2) # Non-linear part of stiffness matrix
sm=ek+eg
tt=TM_GFRM(x1,y1,x2,y2) # Transformation matrix
ck=np.dot(np.dot(tt.T,sm),tt) # Stiffness matrix in global coordinate
ir[5]=3*j+2; ir[4]=ir[5]-1; ir[3]=ir[4]-1
ir[2]=3*i+2; ir[1]=ir[2]-1; ir[0]=ir[1]-1
for i in range(0,nod*nfree):
it=ir[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
df[iz]=0.0
for k in range(0,n):
gk[k,iz]=0.0
gk[iz,iz]=1.0
# solution of simultaneous linear equations
# (calculation of incremental displacement)
dis0 = np.linalg.solve(gk, df)
del gk
# 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):
dis0[iz]=0.0
# Initial parameter setting for Arc-length method
if nnn==1:
ds0=np.sqrt(np.sum(dis0*dis0)+spara*spara*np.sum(df*df))
ds=ds0
dlam0=1.0
dlam=np.sqrt(ds*ds/(np.sum(dis0*dis0)+spara*spara*np.sum(df*df)))
if np.abs(dlam)<0.1*np.abs(dlam0): ds=ds*1.2
if np.abs(dlam0)<np.abs(dlam): ds=ds
dlam=np.sqrt(ds*ds/(np.sum(dis0*dis0)+spara*spara*np.sum(df*df)))
if np.sum(dis_ref*dis0)<0.0: dlam=-dlam
del dis_ref
# Initial displacement for iterative calculation
dis=dlam*dis0
dis_ref=np.zeros(n,dtype=np.float64)
#===========================
# Iterative calculation
for iii in range(1,itmax):
#===========================
lam=lam+dlam
dist=dist+dis
dis_ref=dis_ref+dis
#========================================================
# Calculation of section force & internal force vector
#========================================================
rvec=np.zeros(n,dtype=np.float64) # Internal force vector
for ne in range(0,nele):
#Elimination of rigid roration
i=node[0,ne]-1
j=node[1,ne]-1
#Initial element length
x1=x[0,i]
y1=x[1,i]
x2=x[0,j]
y2=x[1,j]
CL=np.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
tt=TM_GFRM(x1,y1,x2,y2) # Transformation matrix
#Accumulated displacement at previous step in global coordinate system
work[0]=dist[3*i]-dis[3*i]; work[1]=dist[3*i+1]-dis[3*i+1]; work[2]=dist[3*i+2]-dis[3*i+2]
work[3]=dist[3*j]-dis[3*j]; work[4]=dist[3*j+1]-dis[3*j+1]; work[5]=dist[3*j+2]-dis[3*j+2]
wde=np.dot(tt,work)
#Rotation angle at previous step in local coordinate system
Ti=np.tan(wde[2])
Tj=np.tan(wde[5])
theta1i=((CL+wde[3]-wde[0])*Ti-(wde[4]-wde[1]))/((CL+wde[3]-wde[0])+(wde[4]-wde[1])*Ti)
theta1j=((CL+wde[3]-wde[0])*Tj-(wde[4]-wde[1]))/((CL+wde[3]-wde[0])+(wde[4]-wde[1])*Tj)
#Accumulated displacement at current step in global coordinate system
work[0]=dist[3*i]; work[1]=dist[3*i+1]; work[2]=dist[3*i+2]
work[3]=dist[3*j]; work[4]=dist[3*j+1]; work[5]=dist[3*j+2]
wde=np.dot(tt,work)
#Rotation angle at current step in local coordinate system
Ti=np.tan(wde[2])
Tj=np.tan(wde[5])
theta2i=((CL+wde[3]-wde[0])*Ti-(wde[4]-wde[1]))/((CL+wde[3]-wde[0])+(wde[4]-wde[1])*Ti)
theta2j=((CL+wde[3]-wde[0])*Tj-(wde[4]-wde[1]))/((CL+wde[3]-wde[0])+(wde[4]-wde[1])*Tj)
#Rotation angle increment for calculation of internal force
deltai=theta2i-theta1i
deltaj=theta2j-theta1j
# Calculation of section force & internal force vector
#Previous element length
x1=x[0,i]+dist[3*i] -dis[3*i]
y1=x[1,i]+dist[3*i+1]-dis[3*i+1]
x2=x[0,j]+dist[3*j] -dis[3*j]
y2=x[1,j]+dist[3*j+1]-dis[3*j+1]
BL=np.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
#Previous linear stiffness matrix
ek=SM_GFRM(ne,node,ae,x1,y1,x2,y2) # Stiffness matrix in local coordinate
#Current element length
x1=x[0,i]+dist[3*i]
y1=x[1,i]+dist[3*i+1]
x2=x[0,j]+dist[3*j]
y2=x[1,j]+dist[3*j+1]
AL=np.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
work[0]=0.0
work[1]=0.0
work[2]=deltai
work[3]=AL-BL
work[4]=0.0
work[5]=deltaj
#Internal force increment using previous linear stiffness
dfs=np.dot(ek,work)
#Accumulated section force
fsec[:,ne]=fsec[:,ne]+dfs[:]
#Accumulated internal force
ir[5]=3*j+2; ir[4]=ir[5]-1; ir[3]=ir[4]-1
ir[2]=3*i+2; ir[1]=ir[2]-1; ir[0]=ir[1]-1
tt=TM_GFRM(x1,y1,x2,y2) # Transformation matrix
work=np.dot(tt.T,fsec[:,ne]) # Internal force vector
for i in range(0,nod*nfree):
it=ir[i]
rvec[it]=rvec[it]+work[i]
# Un-balanced force
dr=fp+lam*df-rvec
del rvec
#===========================
# Judgement of convergence
#===========================
ic=0
for i in range(0,n):
if 0.0<np.abs(dist[i]):
if 1.0e-3<np.abs(dis[i]/dist[i]): ic=ic+1
if ic==0: break
#===========================
# Update of dlam & dis
#===========================
gk=np.zeros([n,n],dtype=np.float64) # Global stiffness matrix
for ne in range(0,nele):
i=node[0,ne]-1
j=node[1,ne]-1
x1=x[0,i]+dist[3*i]
y1=x[1,i]+dist[3*i+1]
x2=x[0,j]+dist[3*j]
y2=x[1,j]+dist[3*j+1]
P=0.5*(fsec[3,ne]-fsec[0,ne]) # Axial force for non-linear part of stiffness matrix
ek=SM_GFRM(ne,node,ae,x1,y1,x2,y2) # Linear part of Stiffness matrix
eg=P*SG_GFRM(x1,y1,x2,y2) # Non-linear part of stiffness matrix
sm=ek+eg
tt=TM_GFRM(x1,y1,x2,y2) # Transformation matrix
ck=np.dot(np.dot(tt.T,sm),tt) # Stiffness matrix in global coordinate
ir[5]=3*j+2; ir[4]=ir[5]-1; ir[3]=ir[4]-1
ir[2]=3*i+2; ir[1]=ir[2]-1; ir[0]=ir[1]-1
for i in range(0,nod*nfree):
it=ir[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
df[iz]=0.0
dr[iz]=0.0
for k in range(0,n):
gk[k,iz]=0.0
gk[iz,iz]=1.0
# solution of simultaneous linear equations
# (calculation of incremental displacement)
dis0 = np.linalg.solve(gk, df)
disR = np.linalg.solve(gk, dr)
del gk
# 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):
dis0[iz]=0.0
disR[iz]=0.0
dlam=-np.sum(dis0*disR)/(np.sum(dis0*dis0)+spara*spara*np.sum(df*df))
dis=dlam*dis0+disR
#============================
# End of iteration loop
#============================
fp=fp+lam*df
PROUT_FRM(fnameW,nnn,iii,lam,npoin,nele,fp,dist,fsec,dr)
print(nnn,iii,ic,lam)
#============================
# End of loading loop
#============================
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()
Raw
2_inp_gfrm_canti.txt
11 10 1 1 1
200000 100 833
1 2 1
2 3 1
3 4 1
4 5 1
5 6 1
6 7 1
7 8 1
8 9 1
9 10 1
10 11 1
0.000000000 0.0
0.012311659 100.0
0.048943484 200.0
0.108993476 300.0
0.190983006 400.0
0.292893219 500.0
0.412214748 600.0
0.546009500 700.0
0.690983006 800.0
0.843565535 900.0
1.000000000 1000.0
1 1 1 1
11 0.0 -40.0 0.0
Raw
3_inp_gfrm_cable.txt
12 11 2 2 1
2.00E+05 100 833
2.00E+05 28.3 0
1 2 1
2 3 1
3 4 1
4 5 1
5 6 1
6 7 1
7 8 1
8 9 1
9 10 1
10 11 1
11 12 2
0.000000000 0
0.061558295 100
0.244717420 200
0.544967380 300
0.954915030 400
1.464466095 500
2.061073740 600
2.730047500 700
3.454915030 800
4.217827675 900
5.000000000 1000
0.000000000 0
1 1 1 1
12 1 0 1
12 0 -100 0
Raw
4_inp_gfrm_arch.txt
41 40 1 2 1
2.00E+05 100 833
1 2 1
2 3 1
3 4 1
4 5 1
5 6 1
6 7 1
7 8 1
8 9 1
9 10 1
10 11 1
11 12 1
12 13 1
13 14 1
14 15 1
15 16 1
16 17 1
17 18 1
18 19 1
19 20 1
20 21 1
21 22 1
22 23 1
23 24 1
24 25 1
25 26 1
26 27 1
27 28 1
28 29 1
29 30 1
30 31 1
31 32 1
32 33 1
33 34 1
34 35 1
35 36 1
36 37 1
37 38 1
38 39 1
39 40 1
40 41 1
476.8584754 -150.3528998
488.8458402 -105.0225907
496.5342285 -58.76869873
499.8560276 -11.99798689
498.7820251 34.87823687
493.321666 81.4477367
483.5229695 127.3009741
469.4721064 172.0347095
451.2926422 215.2555484
429.1444493 256.5834009
403.2223021 295.6548242
373.7541634 332.126219
340.99918 365.6768508
305.2454037 396.0116709
266.807258 422.8639109
226.0227731 445.9974283
183.2506134 465.208784
138.8669229 480.3290307
93.262018 491.2251989
46.83695426 497.8014662
0 500
-46.83695426 497.8014662
-93.262018 491.2251989
-138.8669229 480.3290307
-183.2506134 465.208784
-226.0227731 445.9974283
-266.807258 422.8639109
-305.2454037 396.0116709
-340.99918 365.6768508
-373.7541634 332.126219
-403.2223021 295.6548242
-429.1444493 256.5834009
-451.2926422 215.2555484
-469.4721064 172.0347095
-483.5229695 127.3009741
-493.321666 81.4477367
-498.7820251 34.87823687
-499.8560276 -11.99798689
-496.5342285 -58.76869873
-488.8458402 -105.0225907
-476.8584754 -150.3528998
1 1 1 1
41 1 1 0
21 0 -100 0
Raw
5_inp_gfrm_lee.txt
21 20 1 2 1
2.00E+05 100 833
1 2 1
2 3 1
3 4 1
4 5 1
5 6 1
6 7 1
7 8 1
8 9 1
9 10 1
10 11 1
11 12 1
12 13 1
13 14 1
14 15 1
15 16 1
16 17 1
17 18 1
18 19 1
19 20 1
20 21 1
0 0
0 100
0 200
0 300
0 400
0 500
0 600
0 700
0 800
0 900
0 1000
100 1000
200 1000
300 1000
400 1000
500 1000
600 1000
700 1000
800 1000
900 1000
1000 1000
1 1 1 0
21 1 1 0
13 0 -100 0
Raw
6_py_fig_gfrm.py
##############################################
# 2D GFRM fig
##############################################
import matplotlib.pyplot as plt
import numpy as np
import sys
# Main routine
args = sys.argv
fnameR=args[1] # input data file
pnode=int(args[2])
dnode=int(args[3])
nodimP=float(args[4])
nodimX=float(args[5])
nodimY=float(args[6])
leg_x=args[7]
leg_y=args[8]
lax_x=args[9]
lax_y=args[10]
lloc=args[11]
f=open(fnameR,'r')
text=f.readline()
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
nnmax=int(text[5])
nod=2 # Number of nodes per element
nfree=3 # Degree of freedom per node
n=nfree*npoin
ip=2
x =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
y =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
node =np.zeros([nod+1,nele],dtype=np.int) # Node-element relationship
fp =np.zeros([nnmax,n],dtype=np.float64) # External force vector
dist =np.zeros([nnmax,n],dtype=np.float64) # Total displacement in global coordinate system
p_dis =np.zeros([ip,nnmax],dtype=np.float64) # Total displacement in global coordinate system
p_lod =np.zeros(nnmax,dtype=np.float64) # Total displacement in global coordinate system
# section characteristics
text=f.readline()
for i in range(0,nsec):
text=f.readline()
# node coordinates
text=f.readline()
for i in range(0,npoin):
text=f.readline()
text=text.strip()
text=text.split()
x[i]=float(text[1]) # x-coordinate
y[i]=float(text[2]) # y-coordinate
# element-node
text=f.readline()
for i in range(0,nele):
text=f.readline()
text=text.strip()
text=text.split()
node[0,i]=int(text[1]) #node_1
node[1,i]=int(text[2]) #node_2
node[2,i]=int(text[3]) #section characteristic number
# Loading step
for nnn in range(0,nnmax):
text=f.readline()
text=f.readline()
for i in range(0,npoin):
text=f.readline()
text=text.strip()
text=text.split()
lp=int(text[0]) #loaded node
fp[nnn,3*lp-3]=float(text[1]) #load in x-direction
fp[nnn,3*lp-2]=float(text[2]) #load in y-direction
fp[nnn,3*lp-1]=float(text[3]) #load in rotation
dist[nnn,3*lp-3]=float(text[4]) #displacement in x-direction
dist[nnn,3*lp-2]=float(text[5]) #displacement in y-direction
dist[nnn,3*lp-1]=float(text[6]) #displacement in rotation
text=f.readline()
for ne in range(0,nele):
text=f.readline()
f.close()
for nnn in range(0,nnmax):
p_dis[0,nnn]=dist[nnn,3*dnode-3]/nodimX # displacement in x-direction
p_dis[1,nnn]=dist[nnn,3*dnode-2]/nodimY # displacement in y-direction
p_lod[nnn]=fp[nnn,3*pnode-2]/nodimP # Load
fnameF=(fnameR.replace('.txt','.png')).replace('out','fig')
if lloc=='UR': slloc='upper right'
if lloc=='UL': slloc='upper left'
if lloc=='LR': slloc='lower right'
if lloc=='LL': slloc='lower left'
plt.plot(p_dis[0,:],p_lod[:],'o-',label=leg_x)
plt.plot(p_dis[1,:],p_lod[:],'^-',label=leg_y)
plt.legend(loc=slloc)
plt.xlabel('Displacement '+lax_x)
plt.ylabel('Load '+lax_y)
plt.grid(which='major',color='#aaaaaa',linestyle='-')
plt.savefig(fnameF, bbox_inches="tight", pad_inches=0.2)
plt.show()
Raw
7_py_fig_gfrm_mode.py
##############################################
# 2D GFRM fig mode
##############################################
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import numpy as np
# Main routine
for icase in range(0,4):
if icase==0: fnameR='out_gfrm_arch.txt' ; mmm=200; lnod=21
if icase==1: fnameR='out_gfrm_cable.txt'; mmm=270; lnod=12
if icase==2: fnameR='out_gfrm_canti.txt'; mmm= 50; lnod=11
if icase==3: fnameR='out_gfrm_lee.txt' ; mmm=100; lnod=13
f=open(fnameR,'r')
text=f.readline()
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
nnmax=int(text[5])
nod=2 # Number of nodes per element
nfree=3 # Degree of freedom per node
n=nfree*npoin
x =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
y =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
node =np.zeros([nod+1,nele],dtype=np.int) # Node-element relationship
fp =np.zeros([nnmax,n],dtype=np.float64) # External force vector
dist =np.zeros([nnmax,n],dtype=np.float64) # Total displacement in global coordinate system
xx =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes after deformation
yy =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes after deformation
# section characteristics
text=f.readline()
for i in range(0,nsec):
text=f.readline()
# node coordinates
text=f.readline()
for i in range(0,npoin):
text=f.readline()
text=text.strip()
text=text.split()
x[i]=float(text[1]) # x-coordinate
y[i]=float(text[2]) # y-coordinate
# element-node
text=f.readline()
for i in range(0,nele):
text=f.readline()
text=text.strip()
text=text.split()
node[0,i]=int(text[1]) #node_1
node[1,i]=int(text[2]) #node_2
node[2,i]=int(text[3]) #section characteristic number
# Loading step
for nnn in range(0,nnmax):
text=f.readline()
text=f.readline()
for i in range(0,npoin):
text=f.readline()
text=text.strip()
text=text.split()
lp=int(text[0]) #loaded node
fp[nnn,3*lp-3]=float(text[1]) #load in x-direction
fp[nnn,3*lp-2]=float(text[2]) #load in y-direction
fp[nnn,3*lp-1]=float(text[3]) #load in rotation
dist[nnn,3*lp-3]=float(text[4]) #displacement in x-direction
dist[nnn,3*lp-2]=float(text[5]) #displacement in y-direction
dist[nnn,3*lp-1]=float(text[6]) #displacement in rotation
text=f.readline()
for ne in range(0,nele):
text=f.readline()
f.close()
# Coordinates after deformation
for i in range(0,npoin):
xx[i]=x[i]+dist[mmm,3*i]
yy[i]=y[i]+dist[mmm,3*i+1]
xmin=np.min([np.min(x),np.min(xx)])
xmax=np.max([np.max(x),np.max(xx)])
ymin=np.min([np.min(y),np.min(yy)])
ymax=np.max([np.max(y),np.max(yy)])
midx=0.5*(xmax+xmin)
midy=0.5*(ymax+ymin)
xmin=midx-(midx-xmin)*1.3
xmax=midx+(xmax-midx)*1.3
ymin=midy-(midy-ymin)*1.3
ymax=midy+(ymax-midy)*1.3
fnameF=(fnameR.replace('.txt','.png')).replace('out_','fig_mode_')
ax=plt.subplot(111)
ax.set_xlim([xmin,xmax])
ax.set_ylim([ymin,ymax])
ax.set_xlabel('x-direction (mm)')
ax.set_ylabel('y-direction (mm)')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
aspect = (ymax-ymin)/(xmax-xmin)*(ax.get_xlim()[1] - ax.get_xlim()[0]) / (ax.get_ylim()[1] - ax.get_ylim()[0])
ax.set_aspect(aspect)
ssxv=np.array([0,-0.5,0.5,0])
ssyv=np.array([0,-1,-1,0])
ssxh=np.array([0,1,1,0])
ssyh=np.array([0,-0.5,0.5,0])
# load for original model
uu=0.0
vv=(ymax-ymin)*0.1
if icase==1:
x1=x[lnod-1]
y1=y[lnod-1]
else:
x1=x[lnod-1]
y1=y[lnod-1]+vv
hl=vv*0.4
hw=hl*0.5
ax.arrow(x1,y1,uu,-(vv-hl), lw=2.0,head_width=hw, head_length=hl, fc='#aaaaaa', ec='#aaaaaa')
# load after deformation
uu=0.0
vv=(ymax-ymin)*0.1
if icase==1:
x1=xx[lnod-1]
y1=yy[lnod-1]
else:
x1=xx[lnod-1]
y1=yy[lnod-1]+vv
hl=vv*0.4
hw=hl*0.5
ax.arrow(x1,y1,uu,-(vv-hl), lw=2.0,head_width=hw, head_length=hl, fc='#555555', ec='#555555')
# Model
for ne in range(0,nele):
i=node[0,ne]-1
j=node[1,ne]-1
x1=x[i]; y1=y[i]
x2=x[j]; y2=y[j]
ax.plot([x1,x2],[y1,y2],color='#aaaaaa',linestyle='-',linewidth=1.0)
ax.plot(x,y,color='#aaaaaa',marker='o',markersize=6)
# displacement
for ne in range(0,nele):
i=node[0,ne]-1
j=node[1,ne]-1
x1=xx[i]; y1=yy[i]
x2=xx[j]; y2=yy[j]
ax.plot([x1,x2],[y1,y2],color='#000000',linestyle='-',linewidth=1.0)
ax.plot(xx,yy,color='#000000',marker='o',markersize=6)
# Boundary conditions
scl=50.0
if icase==0:
lp=41; x1=x[lp-1]; y1=y[lp-1]
ax.plot(x1+ssxv*scl,y1+ssyv*scl,color='#000000',linestyle='-',linewidth=1.5)
lp=1; x1=x[lp-1]; y1=y[lp-1]
px=[x1-2*scl,x1+2*scl,x1+2*scl,x1-2*scl]
py=[y1,y1,y1-1*scl,y1-1*scl]
ax.fill(px,py,fill=False,linewidth=0.0,hatch='///')
ax.plot([x1-2*scl,x1+2*scl],[y1,y1],color='#000000',linestyle='-',linewidth=1.5)
if icase==1:
lp=1; x1=x[lp-1]; y1=y[lp-1]
px=[x1-2*scl,x1+2*scl,x1+2*scl,x1-2*scl]
py=[y1,y1,y1-1*scl,y1-1*scl]
ax.fill(px,py,fill=False,linewidth=0.0,hatch='///')
ax.plot([x1-2*scl,x1+2*scl],[y1,y1],color='#000000',linestyle='-',linewidth=1.5)
if icase==2:
lp=1; x1=x[lp-1]; y1=y[lp-1]
px=[x1-2*scl,x1+2*scl,x1+2*scl,x1-2*scl]
py=[y1,y1,y1-1*scl,y1-1*scl]
ax.fill(px,py,fill=False,linewidth=0.0,hatch='///')
ax.plot([x1-2*scl,x1+2*scl],[y1,y1],color='#000000',linestyle='-',linewidth=1.5)
if icase==3:
lp=1; x1=x[lp-1]; y1=y[lp-1]
ax.plot(x1+ssxv*scl,y1+ssyv*scl,color='#000000',linestyle='-',linewidth=1.5)
lp=21; x1=x[lp-1]; y1=y[lp-1]
ax.plot(x1+ssxh*scl,y1+ssyh*scl,color='#000000',linestyle='-',linewidth=1.5)
plt.savefig(fnameF, bbox_inches="tight", pad_inches=0.2)
plt.show()
plt.clf()
del x,y,node,fp,dist,xx,yy
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment