Page 4 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2021-2025 Page 4 K[idx,idx]:=K[idx,idx]+Array(evalm(transpose(T[n]).Kloc.T[n])):end do:## Solve the system of equationsec := [entries(cc, 'nolist')]:ed := [entries(dd, 'nolist')]:ac := Vector([0, 0, 0]):Kcc := SubMatrix(K, ec, ec):Kcd := SubMatrix(K, ec, ed):Kdc := SubMatrix(K, ed, ec):Kdd := SubMatrix(K, ed, ed):ad := LinearSolve(Kdd, fc - Multiply(Kdc, ac)):qd := Multiply(Kcc, ac) + Multiply(Kcd, ad):## Final displacementsa := Vector(ndof):q := Vector(ndof):a[ec] := ac:q[ec] := qd:a[ed] := ad:N := Vector(nfe):## Axial loads in barsfor n to nfe do N[n] := Multiply(Multiply(Multiply(k[n], Array([-1, 0, 1, 0])), T[n]), a(LM[n])):end do:## Visualisationif StrurelPlot then # begin of plot block if StrurelMode = 0 or StrurelMode = 2 then if StrurelMode = 0 then tit:=\Deterministic Solution\# 0 - initial deterministic solution elif StrurelMode = 2 then tit:=\Stochastic Solution\# 2 - final stochastic solution end if: # Get a screen resolution getmetric:=define_external(\t::integer[4], RETURN::integer[4], LIB=\ ## Plot Function to be used in Maplet TrussPlot := proc(np) # X and Y coordinates in m XY := Matrix([[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]]): if np = 1 then PlotT:=textplot([12.5,6.,\ BarL := Array(1..nfe): # rods as lines BarR := Array(1..nfe): # rod numbers NodN := Array(1..nnp): # nodal numbers # Plot initial state of truss for i to nfe do x1:=XY[IEN[i,1],1]: y1:=XY[IEN[i,1],2]: x2:=XY[IEN[i,2],1]: y2:=XY[IEN[i,2],2]: BarL[i]:=line([x1,y1],[x2,y2],color=blue): # blue color for initial state BarR[i]:=textplot([(x1+x2)/2,(y1+y2)/2,sprintf(\# rod number in the middle of rod length end do: for i to nnp do NodN[i]:=textplot([XY[i,1],XY[i,2],sprintf(\# nodal numbers end do: BarD := Array(1..nfe): # deformed rods as lines # Plot deformed state of truss - scale factor 2 for i to nfe do n1:=IEN[i,1]: n2:=IEN[i,2]: x1:=XY[n1,1]+a[n1*2-1]*2: y1:=XY[n1,2]+a[n1*2]*2: x2:=XY[n2,1]+a[n2*2-1]*2: y2:=XY[n2,2]+a[n2*2]*2: BarD[i]:=line([x1,y1],[x2,y2],color=red): # red color for deformed state end do: DefU:=textplot([XY[4,1],XY[4,2]-0.4,`` || (sprintf(\|| ``]): # value of vertical deflection in node N04 AP := Array(1..nP): AV := Array(1..nP): # Plot vertical forces - scale factor 0.00002 fc := Vector(1..ndof): for n to nP do fc[14 + 2*n] := -p[n]*0.00002: end do: n:=0: for i to nnp do if abs(fc[i*2]) > 0.001 then : n:=n+1: AP[n]:=arrow({[XY[i,1],XY[i,2]-fc[i*2]]},{[0,fc[i*2]]},shape=arrow): # single vector AV[n]:=textplot([XY[i,1],XY[i,2]-fc[i*2],`` || (sprintf(\|| ``],align={above,right}): # force value end if: end do: # Plot set #1 PlotS:=[seq(BarL[i],i=1..nfe),seq(BarR[i],i=1..nfe),seq(NodN[i],i=1..nnp), seq(BarD[i],i=1..nfe),DefU, seq(AP[i],i=1..nP), seq(AV[i],i=1..nP)]: elif np = 2 then PlotT:=textplot([12.5,6.,\ # Plot internal member force N BarN:=Array(1..nfe): # polygon for each bar AN:=Array(1..nfe): # text for each bar pp:=Matrix(2,4): # polygon points for e to nfe do l:=leng[e]/2: n:=abs(N[e])*0.0000005: X:=Array([-l, l, l, -l]): Y:=Array([-n, -n, n, n]): c:=cos(theta[e]): s:=sin(theta[e]): Xrot:= X*c - Y*s: Yrot:= X*s + Y*c: x:=(XY[IEN[e,1],1]+XY[IEN[e,2],1])/2: y:=(XY[IEN[e,1],2]+XY[IEN[e,2],2])/2: pp[1,1..4]:=Xrot+x: pp[2,1..4]:=Yrot+y: if N[e] >= 0 then co:=[1.0,0.6,0.5]: