Page 4 - Demo

Planar TrussExamplefor ComrelAdd-onRCPConsult, 2019-2025Page 4T =nfe*[None]# transformation matrixfore inrange(nfe):idx =np.r_[dof[IEN[e,0],:],dof[IEN[e,1],:]]# extract the dofs of the elementc =np.cos(theta[e])# cosine inclination angles =np.sin(theta[e])# sinus inclination angleT[e]=np.array([[c,s,0,0],# coordinate transformation matrix[s,c,0,0],[0,0,c,s],[0,0,-s,c]])Kloc =k[e]*np.array([[1,0,1,0],# local stiffness matrix[0,0,0,0],[1,0,1,0],[0,0,0,0]])K[np.ix_(idx,idx)]+=T[e].T.dot(Kloc).dot(T[e])# add to K global# Solve systems of equations# f = vector of equivalent nodal forces# q = vector of equilibrium nodal forces# a = displacements#| qd | | Kcc Kcd || ac | | fd |#| | = | || | | |#| qc | | Kdc Kdd || ad | | fc |Kcc =K[cc,:][:,cc]Kcd =K[cc,:][:,dd]Kdc =K[dd,:][:,cc]Kdd =K[dd,:][:,dd]ac =np.array([[0],[0],[0]])# displacements for fixed/support nodesad =np.linalg.solve(Kdd,fc Kdc.dot(ac))# solutionqd =Kcc.dot(ac)+Kcd.dot(ad)## Assemble vectors of displacements (a) and forces (q)a =np.zeros((ndof,1))q =np.zeros((ndof,1))a[cc]=acq[cc]=qda[dd]=ad # q[dd] = qc = 0## Compute axial loadsN =np.zeros(nfe)fore inrange(nfe):idx =np.r_[dof[IEN[e,0],:],dof[IEN[e,1],:]]N[e]=k[e]*np.array([1,0,1,0]).dot(T[e]).dot(a[idx])[0]## Output vertical deflection at bottom center (N04)uout =a[7]## VisualisationifStrurelPlot:# begin of plot blockifStrurelMode==0:tit='Planar Truss Deterministic Solution'# 0 initial deterministic solutionelifStrurelMode==2:tit='Planar Truss Stochastic Solution'# 2 final stochastic solutionifStrurelMode==0orStrurelMode==2:# Get screen resolution and DPIuser32 =ctypes.windll.user32user32.SetProcessDPIAware()screensize =user32.GetSystemMetrics(0),user32.GetSystemMetrics(1)sysDpi =user32.GetDpiForSystem()px =1./sysDpi# X and Y coordinates in mXY =np.array([[0,0],[4,0],[8,0],[12,0],[16,0],[20,0],[24,0],[22,2],[18,2],[14,2],[10,2],[6,2],[2,2]])## Create 3 plotsfornf inrange(3):# Create figure with unique name for each plotplt.figure('%s [%d]'%(tit,nf+1),figsize=[screensize[0]/2*px,screensize[1]/2*px])ifnf ==0:# Set a main subtitle of plotplt.text(12.0,8.0,'Initial and deformed states',fontsize=16,ha='center');bb=np.reshape(a,(nnp,ned))# deflection vectorDEF =XY +bb*2# deformed coordinates with scale factor 2.fd =np.zeros(ndof)fd[15::2]=np.array(P)fd=np.reshape(fd,(nnp,ned))*0.000025# forces with scale factor 0.000025# Plot initial and deformed states of truss with loadsfore inrange(nfe):plt.plot(XY[IEN[e,:],0],XY[IEN[e,:],1],\\# blue color for initial stateplt.plot(DEF[IEN[e,:],0],DEF[IEN[e,:],1],\\# red color for deformed statex=(XY[IEN[e,0],0]+XY[IEN[e,1],0])/2y=(XY[IEN[e,0],1]+XY[IEN[e,1],1])/2plt.text(x,y,'R'+str(e+1),fontsize=8)# rod number in the middle of rod lengthfore inrange(nnp):plt.text(XY[e,0],XY[e,1],'N'+str(e+1),fontsize=8)# nodal numbersplt.text(DEF[3,0],DEF[3,1]0.2,'U = %g m'%bb[3,1],fontsize=8)# value of vertical deflection in node N04n=0fore inrange(nnp):if(fd[e,1]>0.1):n=n+1;plt.arrow(XY[e,0]+fd[e,0],XY[e,1]+fd[e,1],-fd[e,0],-fd[e,1],length_includes_head=True,head_width=0.2,head_length=0.1)# single vectorplt.text(XY[e,0]+0.1+fd[e,0],XY[e,1]+fd[e,1],'P%d=%g N'%(n,fd[e,1]*40000),fontsize=8)# force valueelifnf ==1:# Set a main subtitle of plot