Page 4 - Demo

Planar TrussExamplefor ComrelAdd-onRCPConsult, 2021-2025Page 4@Kcc=@K[c,c]@Kcd=@K[c,d1]@Kdc=@K[d1,c]@Kdd=@K[d-1,d-1]@ad=Numo::Linalg.solve(@Kdd,f@Kdc.dot(@ac));@qd=@Kcc.dot(@ac)+@Kcd.dot(@ad)## Final displacements@a=DFloat.zeros(ndof)@q=DFloat.zeros(ndof)@a[c]=@ac@q[c]=@qd@a[d-1]=@ad## Axial loads in bars@NN=DFloatzeros(nfe)forn in0..nfe-1@NN[n]=k[n]*DFloat[-1,0,1,0].dot(@T[n]).dot(@a[@LM[n,0..]-1])end## Visualizationif$StrurelPlotif$StrurelMode==0or$StrurelMode==2plt=Matplotlib::Pyplotcmx=Matplotlib::cmmpl=Matplotlibendif$StrurelMode==0# Deterministic Solutiontit='Planar Truss Deterministic Solution'# 0 intial deterministic solutionelsif$StrurelMode==2# Stochastic Solutiontit='Planar Truss Stochastic Solution'# 2 final stochastic solutionendif$StrurelMode==0or$StrurelMode==2usr32=Fiddle::dlopen(\gsm=Fiddle::Function.new(usr32[\sw=gsm.call(0);sh=gsm.call(1)dpi=Fiddle::Functionnew(usr32[\_LONG)sysDpi=dpicall;px=1.0/sysDpi# X and Y coordinates in m@XY=DFloat[[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]]ien=Int32.cast(@IEN)## Create 3 plotsfornf in1..3# unique name for each plotfig=pltfigure('%s [%d]'%[tit,nf],figsize:[sw/2*px,sh/2*px])# clear everything in plotfig.clear()ifnf ==1plt.text(12.0,8.0,'Initial and deformed states',fontsize:16,ha'center')# deformed coordinates with scale factor 2.@DEF=DFloat.zeros(nnp,ned)bb=@a.reshape(nnp,ned)@DEF=@XY+bb*2# forces with scale factor 0.000025ff=DFloat.zeros(ndof)forn in0@Plength1ff[15+2*n]=@P[n]endff=ff.reshape(nnp,ned)*0.000025# plot initial and deformed states of truss forn in0..nfe1x=@XY[ien[n,0..]1,0]y=@XY[ien[n,0..]-1,1]plt.plot(x.to_a,y.to_a,'b-')# blue color for initial statex=@DEF[ien[n,0..]1,0]y=@DEF[ien[n,0..]-1,1]plt.plot(x.to_a,y.to_a,'r-')# red color for deformed statex=(@XY[ien[n,0]1,0]+@XY[ien[n,1]1,0])/2y=(@XY[ien[n,0]-1,1]+@XY[ien[n,1]-1,1])/2plt.text(x,y,'R%d'%(n+1),fontsize:8,ha:'center')# rod number in the middle of rod lengthendforn in0..nnp-1plt.text(@XY[n,0],@XY[n,1],'N%d'%(n+1),fontsize:8)# nodal numbersendplt.text(@DEF[3,0],@DEF[3,1]-0.3,'U = %.7f'%bb[3,1],fontsize:8)# value of vertical deflection in node N04# plot vertical forcesn=0foriin0..nnp-1ifff[i,1]>0.1n=n+1plt.arrow(@XY[i,0]+ff[i,0],@XY[i,1]+ff[i,1],-ff[i,0],-ff[i,1],length_includes_head:true,head_width:0.2,head_length:0.1)plt.text(@XY[i,0]+0.1+ff[i,0],@XY[i,1]+ff[i,1],'P%d=%.1f N'%[n,ff[i,1]*40000])# force valueendendelsifnf ==2# set a main subtitle of plotplttext(12.0,8.0,'Internal member forces [N]',fontsize:16,ha:'center')# plot internal member force Nfore in0..nfe1l=leng[e]/2n=@NN[e].abs*0.0000005@X=DFloat[*[l,l,l,l]]@Y=DFloat[*[-n,-n,n,n]]c=Math.cos(theta[e])# cosine inclination angles=Math.sin(theta[e])# sinus inclination anglex=(@XY[ien[e,0]-1,0]+@XY[ien[e,1]-1,0])/2y=(@XY[ien[e,0]-1,1]+@XY[ien[e,1]-1,1])/2@Xrot=@X*c@Y*s+x@Yrot=@X*s+@Y*c+yif@NN[e]>=0co=[1.0,0.6,0.5]elseco=[0.5,0.6,1.0]endplt.fill(@Xrot.to_a,@Yrot.to_a,color:co,ec:'0',l0.5)plt.text(x,y,'%0.1f N'%@NN[e],fontsize:7,ha:'center',va:'center')