Page 9 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2023-2025 Page 9end do!! Material propertiesdo n=1,nfe ke(n)=E(n)*Area(n)/leng(n)end do!! Global stiffness matrixcall MATNUL(K,ndof,ndof) ! cleando n=1,nfe rad=theta(n)*PI/180.d0 ! sin and cos of the inclination c=dcos(rad) s=dsin(rad)! coordinate transformation matrix for the bar n call MATNUL(T,nbd,nbd) T(1,1)= c; T(1,2)=s T(2,1)=-s; T(2,2)=c T(3,3)= c; T(3,4)=s T(4,3)=-s; T(4,4)=c! local stiffness matrix Kloc for the bar n call MATNUL(KL,nbd,nbd) KL(1,1)= ke(n); KL(1,3)=-ke(n) KL(3,1)=-ke(n); KL(3,3)= ke(n)! transpose T in TT call MATRAN(T,TT,nbd,nbd)! get product T'*Kloc*T in KG call MATMLT(TT,KL,KT,nbd,nbd,nbd) call MATMLT(KT, T,KG,nbd,nbd,nbd)! collect coefficients in a global stiffness matrix K do i=1,nbd do j=1,nbd K(LM(i,n),LM(j,n))=K(LM(i,n),LM(j,n))+KG(i,j) end do end doend do!! Boundary conditions! Applied loadscall MATNUL(fc,ndof,1)do i=8,nnp fc(i*ned)=-P(i-7) ! negative vertical forces [N]end do! Supports and free nodes! Remove supports from force vectordo i=1,nc fc(bc(i))=0.d0end do! Remove supports from K matrix - all coefficents on row and col to 0. and diagonal to 1.do i=1,nc do j=1,ndof K(bc(i),j)=0.d0 ! clean a row members end do do j=1,ndof K(j,bc(i))=0.d0 ! clean a col members end do K(bc(i),bc(i))=1.d0 ! diagonal memberend do!! Solve the system of equations K*ac=fccall LINSOL(K,ac,fc,ndof)!! Axial loads in barsdo n=1,nfe ken(1)=-ke(n); ken(2)=0.d0; ken(3)=ke(n); ken(4)=0.d0 rad=theta(n)*PI/180.d0 ! sin and cos of the inclination c=dcos(rad) s=dsin(rad)! coordinate transformation matrix for the bar n call MATNUL(T,nbd,nbd) T(1,1)= c; T(1,2)=s T(2,1)=-s; T(2,2)=c T(3,3)= c; T(3,4)=s T(4,3)=-s; T(4,4)=c! transpose T in TT call MATRAN(T,TT,nbd,nbd) do i=1,nbd keT(i)=0.d0 do j=1,nbd keT(i)=keT(i)+TT(i,j)*ken(j) end do end do do i=1,nbd de(i)=ac(LM(i,n)) end do NN(n)=0.d0 do i=1,nbd NN(n)=NN(n)+keT(i)*de(i) end doend do!!* Visualizationif(StrurelPlot) then if(StrurelMode.eq.0.or.StrurelMode.eq.2) then! Local variables: width, heigth, ouput window number, etc. wo=-2.d0; we=26.d0; ho=-2.d0; he=10.d0 title=\ if(StrurelMode.eq.0) then ! 0 - initial deterministic state state=\ else ! 2 - final stochastic solution state=\ end if! Set of 3 plots do np=1,3! Number of figure if(StrurelMode.eq.0) n=1 nf=np*2-n! Initialize an Gnuplot environment Init=GnuplotInit() if(Init) then! Create an output window call GnuplotWindowRelative(nf,40,40,.TRUE.)! Set a title call GnuplotWindowTitle(nf, trim(title)//c_) write(buffer,\',a,'''')\trim(title),trim(state); call GnuplotPut(nf,trim(buffer)//c_)! Define range for axes write(buffer,\wo,we; call GnuplotPut(nf,trim(buffer)//c_) write(buffer,\ho,he; call GnuplotPut(nf, trim(buffer)//c_) call GnuplotPut(nf, \_) call GnuplotPut(nf,\_)! Plot a main subtitle of plot write(buffer,\& trim(StrurelName),trim(StrurelModule),trim(StrurelVersion); call GnuplotPut(nf,trim(buffer)//c_)