% Submitted by Abhishek Goswami (ME, Materials Engineering, 2012-2014) as
% an assignment; Slightly modified later
%All units are in FPS - American System
%Length and width of the sample is 2 inches and thickness is 0.5 inch and
%the kx = ky = 20 Btu/(hr-ft-F) and h between the sample and the air is 50 Btu/(hr-ft-F)

clc; clear all;
disp('Discretization of domain will conducted in n x n form');
D = input ('Number of elements along width of domain, n: ');

% Known inputs, length/width, thickness, thermal conductivity, convection
% coeffieicent and boundary temperatures

L = 2/12;                  % converting the dimensions in feet
t = 0.5/12;

kx = 20;
ky = 20;
h = 50;

To = 180;   % temperature of the heat source
T_air = 68;   % temperature of air


%Solution
iterations = D*D;          % calculate the total number of discrete elements 
counter = 0;
del_r = 1;
del_c = D+1;
r = 0;
c = 1;

K=zeros((D+1)^2,(D+1)^2);
f=zeros((D+1)^2,1);
T=zeros((D+1)^2,(D+1)^2);
w=zeros(iterations,4);
z=zeros(1,iterations);

% Conductance matrices - from lecture notes
Kh  = (1/iterations)*[0.30864	0.15432	0.07716	0.15432     % convective contribution to the conductance matrix - from lecture notes
0.15432	0.30864	0.15432	0.07716
0.07716	0.15432	0.30864	0.15432
0.15432	0.07716	0.15432	0.30864];

Ke  = [0.55554	-0.13888	-0.27779	-0.13888            % elemental conductance matrix  - from lecture notes
-0.13888	0.55554	-0.13888	-0.27779
-0.27779	-0.13888	0.55554	-0.13888
-0.13888	-0.27779	-0.13888	0.55554];

Kh1 = (2/D)*[0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0];                          % convection matricies which take into account the convection occuring along the edges
Kh2 = (2/D)*[0 0 0 0;0 0 0 0;0 0 0.0579 0.0290;0 0 0.0290 0.0579];         % convection matricies which take into account the convection occuring along the edges
Kh3 = (2/D)*[0 0 0 0; 0 0.0579	0.0290 0 ;0 0.029 0.0579 0; 0 0 0 0];      % convection matricies which take into account the convection occuring along the edges
Kh4 = (2/D)*[0.0579 0.0290 0 0;0.0290 0.0579 0 0;0 0 0 0;0 0 0 0];         % convection matricies which take into account the convection occuring along the edges
fe = (4/iterations)*[11.8056 11.8056 11.8056 11.8056];                     % Lateral convection
f1 = (2/D)*[0 0 0 0];
f2 = (2/D)*[0 0 5.9028 5.9028];
f3 = (2/D)*[0 5.9028 5.9028 0];                                            % element force vectors
f4 = (2/D)*[5.9028 5.9028 0 0];

for j = 1:(D)^2
    if (counter == D)
        r = 1;
        c = c+1;
        counter = 1;
    else
        counter = counter + 1;                                             % assigning the global co-ordinate contribution to the individual elements (relating the local co-ordinates with global co-ordinates)
        r = r+1;
    end
  w(j,1) = 1 + (r-1)*del_r + (c-1)*del_c ;
  w(j,2) = w(j,1) + del_c;
  w(j,3) = w(j,1) + del_c + del_r;
  w(j,4) = w(j,1) + del_r;
   z(1,j)= 0;
   
end

for x = 1:D
    if (x == 1)
        z(1,x) = 1;
    else if (x == D)
          z(1,x) = 2;
        else
            z(1,x) = -1;
        end
    end
    
end

for x = D+1: D :(iterations - D + 1)
    z(1,x) = -4;
    if (x == iterations-D+1 )
        z(1,x)=4;
    end
end


for x = (iterations - D+1): 1: (iterations)
    z(1,x)=-3;
    if (x == iterations)
        z(1,x) = 3;
    end
    if (x == iterations-D+1)
        z(1,x) = 4;
    end
    
end

 for x = D :D: (iterations)
    z(1,x)= -2;
    if (x == iterations)
        z(1,x) = 3;
    end
    if (x == D)
        z(1,x) = 2;
    end
 end   

for j = 1:(D+1)^2
    f(j,1)=0;
    if (j<=(D+1))
        T(j,1)=To;
    else T(j,1)=0;
    end
    for v = 1:(D+1)^2
        K(j,v)=0;
    end
end

for j = 1:iterations
  Keff = Ke + Kh;
  feff = fe;
  if (z(1,j)== 1)
      Keff = Keff + Kh1 + Kh4;
      feff = feff + f1 + f4;                                           %  calculates the effective conductance matrix for each element with every iteration
  else if (z(1,j)== 2)
          Keff = Keff + Kh1 + Kh2;
          feff = feff + f1 + f2;
      else if (z(1,j)== 3)
              Keff = Keff + Kh2 + Kh3;
              feff = feff + f2 +f3;
          else if (z(1,j)== 4)
                  Keff = Keff + Kh3 + Kh4;
                  feff = feff + f3 + f4;
              else if (z(1,j)== -1)
                       Keff = Keff + Kh1;
                       feff = feff + f1;
                  else if (z(1,j)== -2)
                          Keff = Keff + Kh2;
                          feff = feff + f2;
                      else if (z(1,j)== -3)
                              Keff = Keff + Kh3;
                              feff = feff + f3;
                          else if (z(1,j)== -4);
                                  Keff = Keff + Kh4;
                                  feff = feff + f4;
                              end
                          end
                      end
                  end
              end
          end
      end
  end
for x = 1:4
    f (w(j,x),1)= f(w(j,x),1)+ feff(1,x);
   for y = 1:4
       K (w(j,x),w(j,y))= K (w(j,x),w(j,y))+ Keff (x,y);
   end
end
end

Tc = T(1:(D+1),1);
K_act = K((D+2):(D+1)^2,(D+2):(D+1)^2);
f_act = f((D+2):(D+1)^2,1)- K((D+2):(D+1)^2,1:(D+1))*Tc ;                  % relates the elemental matricies to the global conductance matrix
T_act= inv(K_act)*f_act;
T= [Tc ; T_act];

% Plotting and display
Temp = zeros(D+1,D+1);        %Temperature
X = zeros(D+1,D+1);        %x-cordinate
Y = zeros(D+1,D+1);        %y-cordinate

 m=1;
    for j=1:(D+1)
        for n=1:(D+1)                                              
            Temp(n,j)= T(m,1);
            X(n,j)=2/D*(j-1);
            Y(n,j)=2/D*(n-1);
           m=m+1;
        end
    end
    
disp (Temp);     % displays the temperature profile 

xlabel('X (inches)');
ylabel('Y (inches)');

surf(X,Y,Temp);   % plots the contour of the temperature profile
xlabel('X (inch)');
ylabel('Y (inch)');
zlabel('Temperature (Deg F)');