%Me-2040 Numerical methods 
% Steady State 2D heat conduction example - finite difference Example
% using Gauss-Seidel method
% Heated metal notch
clc;
%Set Domain Size and grid size
nx=22;              %number of grid point in x direction
ny=22;              %number of grid points in y direction
Lx = 1;             %length of domain in x-direction
Ly = 1;             %length of domain in y-direction
dx=Lx/(nx-1);       %grid spacing in x-direction
dy=Ly/(ny-1);       %grid spacing in y-direction
C=dx/dy;            %ratio of x and y grid spacing
disp('nx       ny          Lx         Ly        dx          dy')
disp([nx,ny,Lx,Ly,dx,dy])

%Set boundary conditions
T(1:nx,1:ny)=0;
T(1:nx,nx)=0;
T(1,1:ny)=0;
T(ny,1:ny)=1;
T(1:nx,1)=1;

for i=1:nx
    for j=1:ny
    x(i,j)=dx*(i-1);
    y(i,j)=dy*(j-1);
    if x(i,j)>=0.65 & y(i,j)<=0.25
        T(i,j)=1;
    end
    end
end

%Set iteration information
e_s=0.1;
flag =1;
iter=0;

%Main GS loop with error check
while flag > 0
e_check_temp = 0;
e_check = 0;
for i=2:nx-1
    for j=2:ny-1
        %Use Gauss-Seidel to calculate new temperature value
    Tn(i,j)=(T(i+1,j)+T(i-1,j)+C*(T(i,j+1)+T(i,j-1)))/(2*(1+C));
    if Tn(i,j) > 0
    e_check_temp = abs((Tn(i,j)-T(i,j))/Tn(i,j));
    end
    if e_check_temp >= e_check
        e_check = e_check_temp;
    end
    
%Update T value
    T(i,j)=Tn(i,j);
        if x(i,j)>=0.8 & y(i,j)<=0.2
        T(i,j)=1;
        end
    iter=iter+1;
    
    end %end for j
    end %end for i
    
    if e_check <= e_s
        flag= -1;
    end
    
end     %end while loop

disp(['Iterations to converge: ', num2str(iter)]);

contourf(x,y,T);
xlabel('x');
ylabel('y');
zlabel('T');
title('Temperature Contours from GS Calculation');
figure;
surf(x,y,T);
xlabel('x');
ylabel('y');
zlabel('T');
title('Temperature Surface from GS Calculation');