clc
format long e
clear all
syms  x z et
list={'2 edge clamp','1 edge clamp-free','2 edge simply support'};
[selection,ok] = listdlg('ListString',list,'SelectionMode','single','ListSize',[200 60],'Name','Cross Section Type','PromptString','Desired Beam Cross Section');
if selection==1
   bc=[1 1 1 ;         
    1 1 1]; 
elseif selection==2
  bc=[1 1 1 ;         
    0 0 0]; 
else
 bc=[1 1 0 ;         
    1 1 0];
end
 prompt={'Length of Beam(m)','width of Beam(m)','Height of Beam(m)','Number of Elements','Power of FGM','load(N/M)','Upper Theta','Lower Theta','volt'};
    defans={'.025','0.001','0.001','40','1','10000','0','0','0'};
    answer=inputdlg(prompt,'Beam Characteristics',1,defans);
    a=str2double(cell2mat(answer(1,1)));
    b=str2double(cell2mat(answer(2,1)));
    h=0.5*str2double(cell2mat(answer(3,1)));
    mm=str2double(cell2mat(answer(4,1)));
    n=str2double(cell2mat(answer(5,1)));
    p=str2double(cell2mat(answer(6,1)));
    Tu=str2double(cell2mat(answer(7,1)));
    Tl=str2double(cell2mat(answer(8,1)));
    volt=str2double(cell2mat(answer(9,1)));
     prompt={'Elasticity Module','poison ratio','e31','k33','alpha','piro(z)','density','thermal conductivity'};
    defans={'8.134e10','.3','-5.2','1.177e-8','0.2e-5','2.5e-5','7500','75'};
    answer=inputdlg(prompt,'Upper Layer',1,defans);
    Eu=str2double(cell2mat(answer(1,1)));
    nou=str2double(cell2mat(answer(2,1)));
    e31u=str2double(cell2mat(answer(3,1)));
    k33u=str2double(cell2mat(answer(4,1)));
    alphau=str2double(cell2mat(answer(5,1)));
    pirou=str2double(cell2mat(answer(6,1)));
    Rou=str2double(cell2mat(answer(7,1)));
    lambdau=str2double(cell2mat(answer(8,1)));
    
     prompt={'Elasticity Module','poison ratio','e31','k33','alpha','piro(z)','density','thermal conductivity'};
    defans={'6.069e10','.3','-6.5','2.554e-8','1e-5','.548e-5','7500','75'};
    answer=inputdlg(prompt,'Lower Layer',1,defans);
    El=str2double(cell2mat(answer(1,1)));
    nol=str2double(cell2mat(answer(2,1)));
    e31l=str2double(cell2mat(answer(3,1)));
    k33l=str2double(cell2mat(answer(4,1)));
    alphal=str2double(cell2mat(answer(5,1)));
    pirol=str2double(cell2mat(answer(6,1)));
    Rol=str2double(cell2mat(answer(7,1)));
     lambdal=str2double(cell2mat(answer(8,1)));
    %--------------------------------------------------
    x1=a/(2*mm);delx=2*x1; 
%--------------------------------------------
% %Boundary conditions%
% bc=[1 1 1 ;         
%     1 1 1];
% %--------------------------------------------


%------------------------------------------------------------------
%FGM shodane made%
e31=(e31u-e31l)*((z+h)/(2*h))^n+e31l;
k33=(k33u-k33l)*((z+h)/(2*h))^n+k33l;
no=(nou-nol)*((z+h)/(2*h))^n+nol;
E=(Eu-El)*((z+h)/(2*h))^n+El;
alpha=(alphau-alphal)*((z+h)/(2*h))^n+alphal;
piro=(pirou-pirol)*((z+h)/(2*h))^n+pirol;
Ro=(Rou-Rol)*((z+h)/(2*h))^n+Rol;
%------------------------------------------------------------------
%be dast avardane dama dar rastaie z(dar halate estaticki)
lambdak=lambdau-lambdal;
C=1-(lambdak/((n+1)*lambdal))+((lambdak)^2/(((2*n)+1)*(lambdal^2)))-((lambdak)^3/(((3*n)+1)*(lambdal^3)))+((lambdak)^4/(((4*n)+1)*(lambdal^4)))-((lambdak)^5/(((5*n)+1)*(lambdal^5)));
del_teta=Tl+(((Tu-Tl)/C)*((((2*z)+2*h)/(4*h))-((lambdak/((n+1)*lambdal))*(((2*z)+2*h)/(4*h))^(n+1))+(((lambdak^2)/(((2*n)+1)*(lambdal^2)))*(((2*z)+2*h)/(4*h))^((2*n)+1))-(((lambdak^3)/(((3*n)+1)*(lambdal^3)))*(((2*z)+2*h)/(4*h))^((3*n)+1))+(((lambdak^4)/(((4*n)+1)*(lambdal^4)))*(((2*z)+2*h)/(4*h))^((4*n)+1))-(((lambdak^5)/(((5*n)+1)*(lambdal^5)))*(((2*z)+2*h)/(4*h))^((5*n)+1))));
%--------------------------------------------------------------------------
cm=E;
alph=alpha;
beta=cm*alph;
e_p=e31;
piro_f=piro;
ke=k33;
he=z/h;
phi=volt;
%%-----------------------------------------------------------------------
%shoro Finite element
sai1=(.5*(1-et));
sai2=(.5*(1+et));
bu=[-(1/(2*x1)),((-z*(3*et))/(2*(x1^2))),(-z*((3*et)-1))/(2*(x1)),(1/(2*x1)),(z*3*et/(2*(x1^2))),((-z*((3*et)+1))/(2*(x1)))];
bf=-(1/(2*h));
N=[(.5*(1-et)),0,0,(.5*(1+et)),0,0;0,.25*((1-et)^2)*(2+et),.25*x1*(1-(et^2))*(1-et),0,(.25*((1+et)^2)*(2-et)),(.25*x1*(1+et)*(-1+(et)^2));0,(((3*(et^2))-3)/(4*x1)),(((3*(et^2))-(2*et)-1)/(4)),0,(((-3*(et^2))+3)/(4*x1)),((-1+(2*et)+(3*(et^2)))/4)];
H=[1,0,-z;0,1,0];
%%--------------------------------------------------------------------------
% taranehade hameie matris haie lazem
for i=1:3
    for j=1:6
        Nt(j,i)=N(i,j); %#ok<AGROW>
    end
end

for i=1
    for j=1:6
        but(j,i)=bu(i,j); %#ok<AGROW>
    end
end
for i=1:2
    for j=1:3
        Ht(j,i)=H(i,j); %#ok<AGROW>
    end
end
e_pt=e_p;
bft=bf;
%-----------------------------------------------
muu=Nt*Ht*Ro*H*N;          
kuu=but*cm*bu;
kuf=but*e_pt*bf;        
kfu=bft*e_p*bu;
kff=bft*ke*bf;
% %------------------------------------------------
%integration dar jahate x,y
for i=1:6
    for j=1:6
        muu_f(i,j)=int(muu(i,j)*x1*b,et,-1,1); %#ok<AGROW>
        kuu_f(i,j)=int(kuu(i,j)*x1*b,et,-1,1); %#ok<AGROW>
    end
end
for i=1:6
    kuf_f(i,1)=int(kuf(i,1)*x1*b,et,-1,1); %#ok<AGROW>       ?????????????
end
for j=1:6
    kfu_f(1,j)=int(kfu(1,j)*x1*b,et,-1,1); %#ok<AGROW>
end
kff_f(1,1)=int(kff(1,1)*x1*b,et,-1,1);

for i=1:1
    for j=1:6
        bu_f(i,j)=int(bu(i,j)*x1,et,-1,1); %#ok<AGROW>
    end
end
% %--------------------------------------------
%be dast avardan force ha
Htt=[1,0,-h;0,1,0];
fb=Nt*Ht*[0;0];                     %body force%
fsur=Nt*Htt'*[0;p];               %surface force contains distributed load
f_ut=but*beta*del_teta;
%fm=fb+fsur;                            %summation of thermal and mechanical load(one term is droped
f_tphi=bft*piro_f*del_teta;
for i=1:6
    f_f_b(i)=int(fb(i)*x1*b,et,-1,1); %#ok<AGROW,NASGU>
      f_f_sur(i)=int(fsur(i)*x1*b,et,-1,1); %#ok<AGROW>
end
for i=1:6
    f_f_ut(i)=int(f_ut(i)*x1*b,et,-1,1);           
end
    f_f_tphi(1)=int(f_tphi(1)*x1*b,et,-1,1);      
% %--------------------------------------------------------------------------
%integration dar jahate z
%matix KUU koli
kuu_f1=zeros(6,6);
for i=1:6
    for j=1:6
        f1=zeros(1,13);
        hp=kuu_f(i,j);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
        kuu_f1(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;
  %matix MUU koli    
  %muu_f1=zeros(6,6);
        hp=muu_f(i,j);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
        muu_f1(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; %#ok<AGROW>
        
    end
end
%matix KUF koli 
kuf_f1=zeros(6,1);
for i=1:6
    for j=1:1
        f1=zeros(1,13);
        hp=kuf_f(i,j);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
        kuf_f1(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
%matix KFU koli 
 kfu_f1=zeros(1,6);
for i=1:1
    for j=1:6 
        hp=kfu_f(i,j);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
      kfu_f1(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
%matix KFF koli 
   hp=kff_f(1,1);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
      kff_f1(1,1)=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;    
%matix FUT koli 
f_f_ut1=zeros(1,6);
for i=1:1
    for j=1:6
   hp=f_f_ut(i,j);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
      f_f_ut1(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
%matix FTPhi koli 
   hp=f_f_tphi(1,1);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
      f_f_tphi1(1,1)=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;
%matrix F body
f_f_b=zeros(1,6);
for i=1:1
    for j=1:6
   hp=f_f_b(i,j);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end
      f_f_b(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
%matrix Bu
bu_z=zeros(1,6);
for i=1:1
    for j=1:6
   hp=bu_f(i,j);
        hhh=h/6;
        xx=-h;
        kij2=subs(hp,z,xx);
        f1(1)=kij2;
        for ii=2:13
            xx=-h+(ii-1)*hhh;
            kij2=subs(hp,z,xx);
            f1(ii)=kij2;
        end