Page 4 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2021-2025 Page 4@Kcd=@K[c,d-1]@Kdc=@K[d-1,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=DFloat.zeros(nfe)for n in 0..nfe-1 @NN[n]=k[n]*DFloat[-1, 0, 1, 0].dot(@T[n]).dot(@a[@LM[n,0..]-1])end## Visualizationif $StrurelPlot if $StrurelMode == 0 or $StrurelMode == 2 plt=Matplotlib::Pyplot cmx=Matplotlib::cm end if $StrurelMode == 0 # Deterministic Solution tit='Planar Truss - Deterministic Solution' # 0 - intial deterministic solution elsif $StrurelMode == 2 # Stochastic Solution tit='Planar Truss -Stochastic Solution' # 2 - final stochastic solution end if $StrurelMode == 0 or $StrurelMode == 2 usr32=Fiddle::dlopen(\ gsm=Fiddle::Function.new(usr32[\_LONG],Fiddle::TYPE_LONG) sw=gsm.call(0);sh=gsm.call(1) dpi=Fiddle::Function.new(usr32[\_LONG) sysDpi=dpi.call; 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 plots for nf in 1..3 # unique name for each plot plt.figure('%s [%d]' % [tit,nf],figsize: [sw/2*px, sh/2*px]) if nf == 1 plt.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.000025 ff=DFloat.zeros(ndof) for n in 0..@P.length-1 ff[15+2*n]=@P[n] end ff=ff.reshape(nnp,ned)*0.000025 # plot initial and deformed states of truss for n in 0..nfe-1 x=@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 state x=@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 state x=(@XY[ien[n,0]-1,0]+@XY[ien[n,1]-1,0])/2 y=(@XY[ien[n,0]-1,1]+@XY[ien[n,1]-1,1])/2 plt.text(x, y, 'R%d' % (n+1),fontsize: 8, ha: 'center') # rod number in the middle of rod length end for n in 0..nnp-1 plt.text(@XY[n,0], @XY[n,1], 'N%d' % (n+1),fontsize: 8) # nodal numbers end plt.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 forces n=0 for i in 0..nnp-1 if ff[i,1] > 0.1 n=n+1 plt.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 value end end elsif nf == 2 # set a main subtitle of plot plt.text(12.0, 8.0, 'Internal member forces [N]',fontsize: 16,ha: 'center') # plot internal member force N for e in 0..nfe-1 l=leng[e]/2 n=@NN[e].abs*0.0000005 @X=DFloat[*[-l, l, l, -l]] @Y=DFloat[*[-n, -n, n, n]] c=Math.cos(theta[e]) # cosine inclination angle s=Math.sin(theta[e]) # sinus inclination angle x=(@XY[ien[e,0]-1,0]+@XY[ien[e,1]-1,0])/2 y=(@XY[ien[e,0]-1,1]+@XY[ien[e,1]-1,1])/2 @Xrot=@X*c-@Y*s+x @Yrot=@X*s+@Y*c+y if @NN[e] >= 0 co=[1.0, 0.6, 0.5] else co=[0.5, 0.6, 1.0] end plt.fill(@Xrot.to_a, @Yrot.to_a, color: co, ec: '0', l 0.5) plt.text(x, y, '%0.1f N' % @NN[e],fontsize: 7,ha: 'center',va: 'center') end elsif nf == 3 # set a main subtitle of plot plt.text(12.0, 8.0, 'Internal member von Mises stress $\\sigma_{v}$ [MPa]',fontsize: 16,ha: 'center')