function FEM_shrink(a)

% solve the finite element model for shrink fitting a shaft into a rotor
% in the plane stress case

s_d=input('shaft diameter (inches)=');
r_i=input('rotor inner diameter (inches)=');
r_o=input('rotor outer diameter (inches)=');

e_s=input('shaft modulus of elasticity(psi)=');
e_r=input('rotor modulus of elasticity(psi)=');

nu_s=input('shaft Poissons ratio=');
nu_r=input('rotor Poissons ratio=');

alpha_s=input('shaft coefficient of thermal expansion(in/in/degree F)=');
alpha_r=input('rotor coefficient of thermal expansion(in/in/degree F)=');

delta_t=input('thermal change(degree F)=');

s_elem=input('number of elements in the shaft=');
r_elem=input('number of elements in the rotor=');

c1_s=e_s/(1-nu_s^2);
c1_r=e_r/(1-nu_r^2);
c2_s=c1_s*nu_s;
c2_r=c1_r*nu_r;
c3_s=e_s/(1-nu_s);
c3_r=e_r/(1-nu_r);

f1_s=c3_s*alpha_s*delta_t;
f1_r=c3_r*alpha_r*delta_t;

syms a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mesh generation

h_s=s_d/(2*s_elem);
h_r=(r_o-r_i)/(2*r_elem);

elem_x_s=[0:h_s:(s_elem-1)*h_s];% global coordinate of the first node of each element in the shaft
elem_x_r=[r_i/2:h_r:r_o/2-h_r];% global coordinate of the first node of each element in the rotor

elem_n_s=[0:h_s/2:r_i/2-h_s/2];% global coordinate of the nodes in the shaft
elem_n_r=[r_i/2:h_r/2:r_o/2];% global coordinate of the nodes in the rotor


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% obtain the assembled stiffness matrix for quadratic elements

assemb_stiff=zeros((s_elem+r_elem)*2+1);

% compute the element stiffness matrices and assembly in the shaft region

for j=1:s_elem
    
    i=(j-1)*2+1;
    k=stiff(c1_s,c2_s,elem_x_s(j),h_s,a);
    % k is a symbolic object and has to convert to 'double' format
    
    assemb_stiff(i:i+2,i:i+2)=assemb_stiff(i:i+2,i:i+2)+double(k);
    
    
end

% continue computation and assembly in the rotor

for j=1:r_elem
    
    i=(j-1)*2+((s_elem-1)*2+3);
    k=stiff(c1_r,c2_r,elem_x_r(j),h_r,a);
    assemb_stiff(i:i+2,i:i+2)=assemb_stiff(i:i+2,i:i+2)+double(k);
  
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% obtain the assembled load vector 

assemb_load=zeros((s_elem+r_elem)*2+1,1);

% compute the element load vectors and assembly in the shaft region

for j=1:s_elem
    
    i=(j-1)*2+1;
    l=thermal_load(f1_s,elem_x_s(j),h_s,a);
    assemb_load(i:i+2)=assemb_load(i:i+2)+double(l);
    
end

% continue computing and assembly in the rotor

for j=1:r_elem
    
    i=(j-1)*2+((s_elem-1)*2+3);
    l=thermal_load(f1_r,elem_x_r(j),h_r,a);
    assemb_load(i:i+2)=assemb_load(i:i+2)+double(l);
    
end

% the load vector so computed has not contained the stresses due to the boundary
% conditions yet
% the first element of this vector should be added to -sigma_i which is
% the radial stress at the center of the assembly ( i.e the center of the
% shaft ) and is unknown
% also the last element of this vector should be added to sigma_o which is 
% the radial stress at the periphery of the assembly and is equal to zero

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for the unknown (radial) displacements

% the first displacement ( at the center of the shaft ) is equal to zero
% and at the other hand the first element of the load vector is unknown too
% so the inversion is done on first-element-eliminated load and
% displacement vectors ( and correspondingly the first row and column of
% the stiffness matrix are eliminated )

radial_disp=zeros((s_elem+r_elem)*2+1,1);
radial_disp(2:end)=inv(assemb_stiff(2:end,2:end))*assemb_load(2:end);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% obtain the radial and tangential stresses at the nodal points

radial_stress=zeros((s_elem+r_elem)*2+1,1);
tangential_stress=zeros((s_elem+r_elem)*2+1,1);

et_s=alpha_s*delta_t;
et_r=alpha_r*delta_t;

radial_stress(2:s_elem*2)=c1_s*diff(radial_disp(2:s_elem*2+1))/(h_s/2)+c2_s*radial_disp(2:s_elem*2)./elem_n_s(2:end)'-c3_s*et_s;
radial_stress(s_elem*2+1:end-1)=c1_r*diff(radial_disp(s_elem*2+1:end))/(h_r/2)+c2_r*radial_disp(s_elem*2+1:end-1)./elem_n_r(1:end-1)'-c3_r*et_r;

tangential_stress(2:s_elem*2)=c2_s*diff(radial_disp(2:s_elem*2+1))/(h_s/2)+c1_s*radial_disp(2:s_elem*2)./elem_n_s(2:end)'-c3_s*et_s;
tangential_stress(s_elem*2+1:end-1)=c2_r*diff(radial_disp(s_elem*2+1:end))/(h_r/2)+c1_r*radial_disp(s_elem*2+1:end-1)./elem_n_r(1:end-1)'-c3_r*et_r;

% remove the singulairty of the radial and tangential stresses at x=0 ( where
% radial disp. is equal to zero too ) one can use the hopital rule i.e. to
% compute the ratio of diff(radial_disp)/diff(x)
radial_stress(1)=c1_s*diff(radial_disp(1:2))/(h_s/2)+c2_s*diff(radial_disp(1:2))./diff(elem_n_s(1:2))'-c3_s*et_s;
tangential_stress(1)=c2_s*diff(radial_disp(1:2))/(h_s/2)+c1_s*diff(radial_disp(1:2))./diff(elem_n_s(1:2))'-c3_s*et_s;


% compute the tangential stress at the periphery,the derivative of
% radial disp. should be extrapolated using the values of the derivative at 
% the previous nodes 
diff_disp=diff(radial_disp);
diff_disp_length=length(diff_disp);
diff_disp(diff_disp_length+1)=spline([1:diff_disp_length],diff_disp,diff_disp_length+1);
tangential_stress(end)=c2_r*diff_disp(end)/(h_r/2)+c1_r*radial_disp(end)/elem_n_r(end)-c3_r*et_r;

% the radial stress at the periphery is zero

subplot(2,1,1)
plot(-radial_stress)
title('radial stress due to thermal gradient')
ylabel('stress(psi)')
xlabel('displacement from center of assembly(inch)')

subplot(2,1,2)
plot(-tangential_stress)
title('tangential stress due to thermal gradient')
ylabel('stress(psi)')
xlabel('displacement from center of assembly(inch)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ke=stiff(c1,c2,glob_x,h,a)
% compute the quadratic element stiffness matrix 

syms a

shapef_1=1/2*a*(a-1);
shapef_2=1-a^2;
shapef_3=1/2*a*(a+1);

dshf_1=diff(shapef_1,a);
dshf_2=diff(shapef_2,a);
dshf_3=diff(shapef_3,a);

integrand_1_1=(c1*((h/2)*(a+1)+glob_x)*(2/h)^2*dshf_1*dshf_1+c2*(2/h)*(dshf_1*shapef_1+dshf_1*shapef_1)+c1*shapef_1*shapef_1/(glob_x+(h/2)*(a+1)))*(h/2);
integrand_1_2=(c1*((h/2)*(a+1)+glob_x)*(2/h)^2*dshf_1*dshf_2+c2*(2/h)*(dshf_1*shapef_2+dshf_2*shapef_1)+c1*shapef_1*shapef_2/(glob_x+(h/2)*(a+1)))*(h/2);
integrand_1_3=(c1*((h/2)*(a+1)+glob_x)*(2/h)^2*dshf_1*dshf_3+c2*(2/h)*(dshf_1*shapef_3+dshf_3*shapef_1)+c1*shapef_1*shapef_3/(glob_x+(h/2)*(a+1)))*(h/2);
integrand_2_2=(c1*((h/2)*(a+1)+glob_x)*(2/h)^2*dshf_2*dshf_2+c2*(2/h)*(dshf_2*shapef_2+dshf_2*shapef_2)+c1*shapef_2*shapef_2/(glob_x+(h/2)*(a+1)))*(h/2);
integrand_2_3=(c1*((h/2)*(a+1)+glob_x)*(2/h)^2*dshf_2*dshf_3+c2*(2/h)*(dshf_2*shapef_3+dshf_3*shapef_2)+c1*shapef_2*shapef_3/(glob_x+(h/2)*(a+1)))*(h/2);
integrand_3_3=(c1*((h/2)*(a+1)+glob_x)*(2/h)^2*dshf_3*dshf_3+c2*(2/h)*(dshf_3*shapef_3+dshf_3*shapef_3)+c1*shapef_3*shapef_3/(glob_x+(h/2)*(a+1)))*(h/2);

ke(1,1)=int(integrand_1_1,a,-1,1);
ke(1,2)=int(integrand_1_2,a,-1,1);
ke(1,3)=int(integrand_1_3,a,-1,1);
ke(2,2)=int(integrand_2_2,a,-1,1);
ke(2,3)=int(integrand_2_3,a,-1,1);
ke(3,3)=int(integrand_3_3,a,-1,1);

ke(2,1)=ke(1,2);
ke(3,1)=ke(1,3);
ke(3,2)=ke(2,3);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function le=thermal_load(f1,glob_x,h,a)
% compute the quadratic element load vector due to thermal gradient

syms a

shapef_1=1/2*a*(a-1);
shapef_2=1-a^2;
shapef_3=1/2*a*(a+1);

dshf_1=diff(shapef_1,a);
dshf_2=diff(shapef_2,a);
dshf_3=diff(shapef_3,a);

integrand_1=(shapef_1*f1+dshf_1*(2/h)*f1*(glob_x+(h/2)*(a+1)))*(h/2);
integrand_2=(shapef_2*f1+dshf_2*(2/h)*f1*(glob_x+(h/2)*(a+1)))*(h/2);
integrand_3=(shapef_3*f1+dshf_3*(2/h)*f1*(glob_x+(h/2)*(a+1)))*(h/2);

le(1,1)=int(integrand_1,a,-1,1);
le(2,1)=int(integrand_2,a,-1,1);
le(3,1)=int(integrand_3,a,-1,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%