Page 4 - Demo

Planar Truss Examplefor ComrelAdd-onRCPConsult, 2021-2025Page 4T[[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]]].KlocT[[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,fcKdc.ac];qd=Kccac+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*21]]*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;X={l,l,l,l};