%%%%%% MAIN Program %%%%%%%%%%%%%%%% % A simple fe program for 2D bar elements %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all clc % Read Input Data FE_Input1 % Calculating length and angle of each elements [NN,NNe]=size(Nodes); [Ne,Nee]=size(Elements); for i=1:Ne Node1=Elements(i,2); xnode1=Nodes(Node1,2); ynode1=Nodes(Node1,3); Node2=Elements(i,3); xnode2=Nodes(Node2,2); ynode2=Nodes(Node2,3); Le(i)=sqrt((xnode2-xnode1)^2+(ynode2-ynode1)^2); Tete(i)=atan((ynode2-ynode1)/(xnode2-xnode1)); end Kglobal=zeros(NN*2,NN*2); % Calculate the connectivity matrix conect=[]; for i=1:Ne conect=[conect Elements(i,2)*2-1 Elements(i,2)*2 Elements(i,3)*2-1 Elements(i,3)*2 ] pause end % Calculate the stiffness matrix for each element (Ke) % and then assemble it to the global stiffness matrix Kglobal for i=1:Ne Ke=Kelement(Eproperties(i,2),Eproperties(i,3), Le(i),Tete(i)); for j1=1:4 for j2=1:4 m1=conect(i,j1); m2=conect(i,j2); Kglobal(m1,m2)=Kglobal(m1,m2)+Ke(j1,j2); end end end % Calculate the index for Boundary conditions [NBC,MBCN]=size(BCN); XD=[]; for i=1:NBC if BCN(i,2)==1 XD=[XD BCN(i,1)*2-1]; end if BCN(i,3)==1 XD=[XD BCN(i,1)*2]; end end % Calculate those BCs which are not fixed XDD=[]; S=[1:NN*2]; for i=1:NN*2 aaa=find(S(i)==XD); if isempty(aaa)==1; XDD=[XDD S(i)]; end end % Calculate the force vector Fglobal Fglobal=zeros(NN*2,1); [LN1,LM1]=size(LoadN); for i=1:LN1 Fglobal(LoadN(i,1)*2-1)=LoadN(i,2); Fglobal(LoadN(i,1)*2)=LoadN(i,3); end % Calculate stiffness matrix and force vector for non-fixed DOFs KglobalR=Kglobal(XDD,XDD) FglobalR=Fglobal(XDD) % Calculate the displacement of non-fixed DOFs U=inv(KglobalR)*FglobalR