Page 10 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2023-2025 Page 10 if(StrurelMode.eq.0) then! Input data at state of mean values call GnuplotPut(nf, \_) else! Final data of stochastic solution write(buffer,\& trim(StrurelIMET),trim(StrurelIOPT); call GnuplotPut(nf,trim(buffer)//c_) write(buffer,\& StrurelBeta,StrurelPf; call GnuplotPut(nf,trim(buffer)//c_) end if select case(np) case(1) !! plot set #1! Set a main subtitle of plot call GnuplotPut(nf,\_)! Plot initial state of truss sc=2.d0 do n=1,nfe v(1,1)=XY(1,IEN(1,n)); v(2,1)=XY(2,IEN(1,n)) v(1,2)=XY(1,IEN(2,n)); v(2,2)=XY(2,IEN(2,n))! blue color for initial state call GnuplotPlotData(nf,\_,v,2,2)! rod number in the middle of rod length write(buffer,\& n,(v(1,1)+v(1,2))/2,(v(2,1)+v(2,2))/2; call GnuplotPut(nf,trim(buffer)//c_)! red color for deformed state v(1,1)=v(1,1)+ac(IEN(1,n)*ned-1)*sc v(2,1)=v(2,1)+ac(IEN(1,n)*ned )*sc v(1,2)=v(1,2)+ac(IEN(2,n)*ned-1)*sc v(2,2)=v(2,2)+ac(IEN(2,n)*ned )*sc call GnuplotPlotData(nf,\_,v,2,2) end do! nodal numbers do i=1,nnp write(buffer,\& i,XY(1,i),XY(2,i); call GnuplotPut(nf,trim(buffer)//c_) end do! Forces with scale factor 0.000025 call MATNUL(fp,ned,nnp) do i=8,nnp fp(2,i)=P(i-7)*0.000025 v(1,1)=XY(1,i); v(2,1)=XY(2,i)+fp(2,i) v(1,2)=0.d0; v(2,2)=-fp(2,i)! Single vector call GnuplotPlotData(nf,\_,v,4,1)! Force value write(buffer,\ i-7,fp(2,i)*40000.d0,XY(1,i)+0.1d0,XY(2,i)+fp(2,i); call GnuplotPut(nf,trim(buffer)//c_) end do! Value of vertical deflection in node N04 write(buffer,\& ac(8),XY(1,4),XY(2,4)+ac(8)*5.d0; call GnuplotPut(nf,trim(buffer)//c_) case(2) !! plot set #2! Set a main subtitle of plot call GnuplotPut(nf,\_)! Plot internal member force N do n=1,nfe l1=leng(n)/2.d0 ! length of rectangle l2=dabs(NN(n))*0.0000005d0 ! height of rectangle X(1)=-l1; X(2)=l1; X(3)=l1; X(4)=-l1 Y(1)=-l2; Y(2)=-l2; Y(3)=l2; Y(4)=l2 rad=theta(n)*PI/180.d0 c=dcos(rad) s=dsin(rad) xc=(XY(1,IEN(1,n))+XY(1,IEN(2,n)))/2 yc=(XY(2,IEN(1,n))+XY(2,IEN(2,n)))/2 do i=1,nbd XYrot(1,i)=X(i)*c-Y(i)*s+xc XYrot(2,i)=X(i)*s+Y(i)*c+yc end do! Need to be closed XYrot(1,nbd+1)=XYrot(1,1) XYrot(2,nbd+1)=XYrot(2,1)! Color of polygon if(NN(n).gt.0.d0) then co=\ else co=\ end if write(buffer,\ notitle')\co call GnuplotPlotData(nf, trim(buffer)//c_,XYrot,ned,nbd+1) write(buffer,\& NN(n),xc,yc; call GnuplotPut(nf,trim(buffer)//c_) end do case(3) !! plot set #3! Set a main subtitle of plot call GnuplotPut(nf, & \_v [MPa]' at 12.,8. center font 'Arial,16'\! Prepare color mapping - take Jet palette call GnuplotPut(nf, & \ & 5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')\_)! Find minimal and maximal values of von Mises stress minA=1.d6; maxA=-1.d6 do n=1,nfe vs=dabs(NN(n)/Area(n)*1.d-6) minA=dmin1(minA, vs) maxA=dmax1(maxA, vs) end do ddA=maxA-minA! Prepare colorbar data call GnuplotPut(nf,\_) write(buffer,\minA,maxA; call GnuplotPut(nf, trim(buffer)//c_)! Plot of von Mises stress for internal member do n=1,nfe Sv=dabs(NN(n)/Area(n)*1.d-6) l1=leng(n)/2.d0 ! length of rectangle l2=Sv/ddA*0.2d0 ! height of rectangle X(1)=-l1; X(2)=l1; X(3)=l1; X(4)=-l1 Y(1)=-l2; Y(2)=-l2; Y(3)=l2; Y(4)=l2 rad=theta(n)*PI/180.d0 c=dcos(rad) s=dsin(rad) xc=(XY(1,IEN(1,n))+XY(1,IEN(2,n)))/2 yc=(XY(2,IEN(1,n))+XY(2,IEN(2,n)))/2 do i=1,nbd XYten(1,i)=X(i)*c-Y(i)*s+xc XYten(2,i)=X(i)*s+Y(i)*c+yc