fgpmbeam2.m

skyline method to inverse K matrix in [K]{X}={F}

      bu_z(i,j)=hhh*(f1(1)+4*f1(2)+2*f1(3)+4*f1(4)+2*f1(5)+4*f1(6)+2*f1(7)+4*f1(8)+2*f1(9)+4*f1(10)+2*f1(11)+4*f1(12)+f1(13))/3;
    end
end
% %--------------------------------------------------------------------------
%tashkil matrix F koli
f_f_m=f_f_b+f_f_sur;
f_f_ut2=f_f_ut1';
 f_f_tphi2= f_f_tphi1';
 f_f_m1=f_f_m';
 FK=f_f_m1+f_f_ut2-( kuf_f1*kff_f1'*f_f_tphi2);      %F koli ke shamele F mekaniki+ F(u phi) - ebarate marbot be F (phi theta)
%tashkile matrix K koli
KK=kuu_f1+(kuf_f1*kff_f1'*kfu_f1);
%tashkile matrix M koli
MK=muu_f1;
%--------------------------------------------------------------------------
%dadan mokhtasat be har nod
node=zeros((mm+1),1);
for i=1:(mm+1)
    node(i,1)=(i-1)*delx;
end
%--------------------------------------------------------------------------
%moshakhas kardane elemanha ba nodhaie atraf%
element=zeros(mm,1);
for i=1:mm
    element(i,1)=i;
    element(i,2)=i;
    element(i,3)=i+1;
end
%--------------------------------------------------------------------------
%3 daraje azadi darim ...matrix haie sefr ra ijad mikonim%
ka=zeros(3*(mm+1),3*(mm+1));
keij=zeros(3*(mm+1),3*(mm+1));
ma=zeros(3*(mm+1),3*(mm+1));
maij=zeros(3*(mm+1),3*(mm+1));
da=zeros(3*(mm+1),3*(mm+1));
daij=zeros(3*(mm+1),3*(mm+1));
 fa=zeros(1,3*(mm+1));
 v0=zeros(1,3*(mm+1));
 d0=zeros(1,3*(mm+1));
for i=1:3*(mm+1)
    for j=1:3*(mm+1)
        ka(i,j)=0;
        keij(i,j)=0;
        ma(i,j)=0;
        maij(i,j)=0;
        da(i,j)=0;
        daij(i,j)=0;
    end
    fa(i)=0;
    d0(i)=0;
    v0(i)=0;
end
%--------------------------------------------------------------------------
d0=d0'; v0=v0'; fa=fa'; 
fe=[0;0;0;0;0;0]; 
FK=sym2poly(FK);
   fe=FK+fe;        
for ii=1:mm
    i=element(ii,2);
    j=element(ii,3); 
    xi=node(ii,1);
%%%%%%%%%%%%%%%%%
    for v=1:2
        for w=1:2
            keij(3*element(ii,v+1)-2:3*element(ii,v+1),3*element(ii,w+1)-2:3*element(ii,w+1))=KK(3*v-2:3*v,3*w-2:3*w);
            ka=ka+keij;
            keij=zeros(3*(mm+1),3*(mm+1));
            maij(3*element(ii,v+1)-2:3*element(ii,v+1),3*element(ii,w+1)-2:3*element(ii,w+1))=MK(3*v-2:3*v,3*w-2:3*w);
            ma=ma+maij;
            maij=zeros(3*(mm+1),3*(mm+1));
        end
    end
%%%%%%%%%%%%%%%
    fa(3*i-2)=fe(1)+fa(3*i-2);
    fa(3*i-1)=fe(2)+fa(3*i-1);
    fa(3*i)=fe(3)+fa(3*i);  
    fa(3*j-2)=fe(4)+fa(3*j-2);
    fa(3*j-1)=fe(5)+fa(3*j-1);
    fa(3*j)=fe(6)+fa(3*j);
end
faa=fa;   
 kaa=ka;
 maa=ma;
 %-----------------------------------------------

    %---------------------------------------------------------
% %boundary condition 
%finding of nodes that should be constrained
%ll=left line of beam; rl=right line of beam
jj=1;
for i=1:3
    if bc(1,i)==1
            del(jj)=i; %#ok<AGROW>
            jj=jj+1;
    end
    if bc(2,i)==1
            del(jj)=3*(mm)+i; %#ok<AGROW>
            jj=jj+1;
    end
end
del1=del;
del=sort(del);
bc_size=size(del);
%hazfe derayeheye tekrariye del
l=1;
gg=0;
% %boundary condition without elimineting method
j=1;
for i=1:3*(mm+1)
    if i==del(j)          
            ka(i,:)=0;
             ka(:,i)=0;
            ka(i,i)=1;
           
            ma(i,:)=0;
            ma(:,i)=0;
            ma(i,i)=1;
            
            fa(i)=0;
         j=j+1;
         if j==bc_size(2)+1
             j=j-1;
         end
    end
end
%boundary condition with rows and colomns elimineting method
for i=1:3*(mm+1)
    k=0;
    for j=1:bc_size(2)
        if i==del(j);
            k=1;
        end
    end
    if k==0
        red_ka1(l,:)=ka(i,:); %#ok<AGROW>
        red_fa(l)=fa(i); %#ok<AGROW>
        red_ma1(l,:)=ma(i,:); %#ok<AGROW>
        l=l+1;
        gg=gg+1;
    end
end
ggg=0;
l=1;
for i=1:3*(mm+1)
    k=0;
    for j=1:bc_size(2)
        if i==del(j)
            k=1;         
        end
    end
    if k==0
        red_ka(:,l)=red_ka1(:,i); %#ok<AGROW>
        red_ma(:,l)=red_ma1(:,i); %#ok<AGROW>
        l=l+1;
        ggg=ggg+1;
    end
end

%bedast avardan javab ha az ravesh boundary aval
red_fa=red_fa';
res1=inv(ka)*fa;
res1=res1/(b);
s=0;
rest2=inv(red_ka)*red_fa;
def_x=zeros(1,(mm+1));
def_z=zeros(1,(mm+1));
slip_x=zeros(1,(mm+1));
for i=1:3:3*(mm+1)
    s=s+1;
    def_x(s)=res1(i);
    def_z(s)=res1(i+1);
    slip_x(s)=res1(i+2);
end



eme1=(inv(ma))*ka;
omeg=eig(eme1);
omeghaa=sort(sqrt(omeg));
%Dimentionless_omegha=omeghaa*2*h*(sqrt(Rol/El));
%be dast avardan ferekans az ravashe boundary dovom
E=subs(E,z,0);
Ro=subs(Ro,z,0);
A=b*2*h;
I=(1/12)*b*(2*h)^3;
Dimentionless_omegha=omeghaa*a*(sqrt(I/A));
W_moghavemat_2sarsimply=5*p*(a^4)/(384*E*I);                 %do sar gir dar
Omegh_Rao1_2sargir=(4.7123^2)*(E*I/(Ro*A*(a^4)))^.5  ;  % do sar gir dar
Omegh_Rao2_2sargir=(7.8539^2)*(E*I/(Ro*A*(a^4)))^.5 ;
Omegh_Rao3_2sargir=(10.9955^2)*(E*I/(Ro*A*(a^4)))^.5 ;
Omegh_Rao1_2sarsimply=(3.1416^2)*(E*I/(Ro*A*(a^4)))^.5 ;   % 2sar simply
Omegh_Rao2_2sarsimply=(6.2832^2)*(E*I/(Ro*A*(a^4)))^.5 ;   % 2sar simply
Omegh_Rao3_2sarsimply=(9.4247^2)*(E*I/(Ro*A*(a^4)))^.5 ;
W_moghavemat_1sargir=p*(a^4)/(8*E*I)  ;               %do sar gir dar
W_moghavemat_2sargir=p*(a^4)/(384*E*I)  ;               %do sar gir dar
theta_1sargirdar=p*(a^3)/(6*E*I) ;
theta_2sarsimply=p*(a^3)/(24*E*I) ;
Dimentionless_omegha_Rao1_2sargir=Omegh_Rao1_2sargir*a*(sqrt(I/A));
Dimentionless_omegha_Rao2_2sargir=Omegh_Rao2_2sargir*a*(sqrt(I/A));
Dimentionless_omegha_Rao3_2sargir=Omegh_Rao3_2sargir*a*(sqrt(I/A));
Dimentionless_omegha_Rao1_2sarsimply=Omegh_Rao1_2sarsimply*a*(sqrt(I/A));
Dimentionless_omegha_Rao2_2sarsimply=Omegh_Rao2_2sarsimply*a*(sqrt(I/A));
Dimentionless_omegha_Rao3_2sarsimply=Omegh_Rao3_2sarsimply*a*(sqrt(I/A));
%-----------------------------------------------------------------
%tarif WW ha
nn=a/mm;
x=0:nn:a;
%W_1sargirdar
i=1;
ww_1s=zeros(1,(mm+1));
for j=1:mm+1
      ww_1s(i,j)=((-4*p*a*(x(j)^3)+(p*(x(j)^4))+6*p*(x(j)^2)*(a^2))/(24*E*I));
end
ww_1s=subs(ww_1s,z,0);            % deflection tare vasat ra mikhahim 
%-----------------%
  %W_2sargirdar
  i=1;
  ww_2s=zeros(1,(mm+1));
for j=1:mm+1
      ww_2s(i,j)=((p*(x(j)^4))+(p*(a^2)*(x(j)^2))-(2*p*(x(j)^3)*a))/(24*E*I);
end
ww_2s=subs(ww_2s,z,0);             % deflection tare vasat ra mikhahim 
%------------------%
  %W_2sarsimply
  i=1;
  ww_2si=zeros(1,mm+1);
for j=1:mm+1
    ww_2si(i,j)=(p*x(j)*(a^3-(2*a*(x(j)^2))+(x(j)^3))/(24*E*I));
end
ww_2si=subs(ww_2si,z,0);            % deflection tare vasat ra mikhahim 
%--------------------------------------------------------------------------
%Result 1
x=x/a;
% if bc(4)==1
%     if bc(6)==1
%         plot(x,def_z,'*',x,ww_2s)
%         xlabel('Lenght of Beam')
%          ylabel('Deflection of Beam')
%          hold on
%     else
%         plot(x,def_z,'*',x,ww_2si)
%         xlabel('Lenght of Beam')
%         ylabel('Deflection of Beam')
%        hold on
%     end
% else
%      plot(x,def_z,'*',x,ww_1s)
%         xlabel('Lenght of Beam')
%          ylabel('Deflection of Beam')
%          hold on
% end

if bc(4)==1
    if bc(6)==1
        plot(x,def_z)
        xlabel('Lenght of Beam')
         ylabel('Deflection of Beam')
         hold on
    else
        plot(x,def_z)
        xlabel('Lenght of Beam')
        ylabel('Deflection of Beam')
       hold on
    end
else
     plot(x,def_z)
        xlabel('Lenght of Beam')
         ylabel('Deflection of Beam')
         hold on
end

%--------------------------------------------------------------------------

%be dast avardane stress va E electricite
EE=-bf*phi;              %bordare meidane electriki
EK=zeros(1,mm);
for i=1:mm
    EK(i)=EE;
end
j=1;
for i=3:3:3*(mm)   
    UU(j,:)=[res1(i-2),res1(i-1),res1(i),res1(i+1),res1(i+2),res1(i+3)]; %#ok<AGROW>
     j=j+1;
end

E=subs(E,z,0);
e31=subs(e31,z,0);
k33=subs(k33,z,0);
epsilon=bu_z*UU';
Stress=(E*epsilon)+(-e31*EK);
DZ=(e31*epsilon)+(k33*EK);
%Stress=sym2poly(Stress);
%DZ=sym2poly(DZ);