Page 4 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2021-2025 Page 4 (*! coordinate transformation matrix for the bar n *) T[[n]]={{c, s, 0, 0},{-s, c, 0, 0},{0, 0, c, s},{0, 0, -s, c}}; (*! local stiffness matrix for the bar n *) Kloc=k[[n]]*{{1, 0, -1, 0},{0, 0, 0, 0},{-1, 0, 1, 0},{0, 0, 0, 0}}; (*! global stiffness matrix assembly *) K[[LM[[n]],LM[[n]]]]+=Transpose[T[[n]]].Kloc.T[[n]];];(*# Solve the system of equations *)ac={0,0,0};Kcc=K[[cc,cc]];Kcd=K[[cc,dd]];Kdc=K[[dd,cc]];Kdd=K[[dd,dd]];ad=LinearSolve[Kdd,fc-Kdc.ac];qd=Kcc.ac+Kcd.ad;(*# Final displacements *)a=Array[0&,ndof];q=Array[0&,ndof];a[[cc]]=ac;q[[cc]]=qd;a[[dd]]=ad;(*# Axial loads in bars *)NN=Array[0&,nfe];For[n=1,n<=nfe,n++, NN[[n]]=k[[n]]*{-1, 0, 1, 0}.T[[n]].a[[LM[[n]]]];];(*# Visualisation *)If[StrurelPlot, (*! begin of plot block *) If[StrurelMode == 0 || StrurelMode == 2, (*# No Front End in MathKernel - use Java graphics instead *) Needs[\ (*# Get a screen resolution *) gm=DefineDLLFunction[\ (*! Plot title *) If[StrurelMode == 0, title=\Deterministic Solution\ title=\Stochastic Solution\ (*! X and Y coordinates in m *) XY={{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}}; For[p=1,p<=3,p++, Switch[p, 1, (*# plot set #1 *) tit=Text[Style[Row[{Style[\Blue],\Red],\18], {14,6}, Center]; BarL=Array[0&,nfe]; (*! rods as lines *) BarR=Array[0&,nfe]; (*! rod numbers *) NodN=Array[0&,nnp]; (*! nodal numbers *) (*! plot initial state of truss *) For[i=1,i<=nfe,i++, 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}}]; (*! rod line *) BarR[[i]]=Text[\{(x1+x2)/2, (y1+y2)/2}]; (*! rod number in the middle of rod length *) ]; For[i=1,i<=nnp,i++, NodN[[i]]=Text[\{XY[[i,1]], XY[[i,2]]}]; (*! nodal numbers *) ]; DefU=Text[Style[Row[{\{XY[[4,1]], XY[[4,2]]-0.4}]; (*! value of vertical deflection in node N04 *) BarD=Array[0&,nfe]; (*! plot deformed state of truss - scale factor 2 *) For[i=1,i<=nfe,i++, 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}}]; (*! deformed rod line *) ]; (*! plot vertical forces - scale factor 0.00002 *) fc=Array[0&,ndof]; For[n=1,n<=Length[P],n++,fc[[14+2*n]]=-P[[n]]*0.00002;]; AP=Array[0&,Length[P]]; AV=Array[0&,Length[P]]; n=0; For[i=1,i<=nnp,i++, If[Abs[fc[[i*2]]] > 0.001, n=n+1; AP[[n]]=Arrow[{{XY[[i,1]],XY[[i,2]]-fc[[i*2]]},{XY[[i,1]],XY[[i,2]]}}]; (*! single vector *) AV[[n]]=Text[Style[Row[{\{XY[[i,1]],XY[[i,2]]-fc[[i*2]]}, Left]; (*! force value *) ]; ]; base={Blue,BarL,Black,BarR,NodN, (*! initial state in blue *) DefU,Red,BarD, (*! deformed state in red *) Black,Arrowheads[0.01],AP,AV}; (*! Loading in black *) ep={};, 2, (*# plot set #2 *) tit=Text[Style[\18], {14,6}, Center]; BarN=Array[0&,nfe]; (*! polygon for each bar *) AN=Array[0&,nfe]; (*! text for each bar *) (*! plot internal member force N *) For[e=1,e<=nfe,e++, l=leng[[e]]/2; n=Abs[NN[[e]]]*0.0000005;