%Jeremy Cuddihy
%Successive Over Relaxation Method
%Assignment 12
%ME 504 Computational Fluid Dynamics
%____________________________________________________________________
clc; clear all; close all;
%____________________________________________________________________  
%%

%Grid Specifications
%If Grid Spacing Is Uniform, Ratio is Equal To 1
m = 30; n = 30; %m = rows, n=columns
delta_x = 1;
delta_y = 1;
ratio=delta_x/delta_y;
i_max = m;
j_max = n;

%____________________________________________________________________  
%%

%Setting Temperature Constraints 
top =300;             %Temperature of top of Plate
left=1000;             %Temperature of Left Side Of Plate
right=300;            %Temperature of Right Side of Plate
bottom =500;          %Temperature of bottom of plate
%____________________________________________________________________
%%

%Solving For "Q" Matrix
N = (m-2)*(n-2);    %Outside Boundaries Of Plate Are Defined
                    %Eliminate Outside Boundaries
                    %Q is product of rows*columns with eliminated
Q = zeros(N,1);     %Preallocates Q Vector For Iteration Speed
for i = 2:(m-3)     %Left and Right Side Plate Grid Points
    Q(i,1) = -left ;
    Q(i+((m-2)*(n-3)),1) = -right ;
end
for i = 1:(n-4)     %Top and Bottom Grid Points
    Q((i*(m-2)+1),1) = -bottom ;
    Q(((i+1)*(m-2)),1) = -top ;
end
Q(1,1) = -(left+bottom) ;   %Bottom Left Corner
Q(m-2,1) = -(left+top) ;    %Top Left Corner
Q(((m-2)*(n-3)+1),1) = -(right+bottom) ;    %Bottom Right Corner
Q(((m-2)*(n-2)),1) = -(right+top) ; %Top Right Corner

%____________________________________________________________________
%%

%Solving For "A" Coefficient Matrix

%If Square Grid Matrix (i.e. if m=n)
 I = speye(m-2,n-2);
 E = sparse(2:m-2,1:(n-2)-1,1,m-2,n-2);
 D = E+E'-2*I;
 a = kron(D,I)+kron(I,D);
 d = zeros((m-2)^2,(n-2)^2);
 A = d+a;
 
%If Not Square Grid Matrix (i.e. m and n aren't equal)

% A = zeros(N,N) ;      %preallocates A Matrix
% for i = 1:N
%     j = i ;
%     A(i,j) = -4 ;     %coefficient for grid point of interest
%     if i <= (N-(m-2)) %Coefficient for grid pt. to right
%         A(i,j+(m-2)) = 1 ;
%     end
%     if i > (m-2)      %Coefficient for grid pt. to left
%         A(i,j-(m-2)) = 1 ;
%     end
%     if i < N
%         A(i+1,j) = 1 ;
%         A(i,j+1) = 1 ;
%     end  
%     A;
% end
% for i = 1:(n-3)
%     j = i ;   
%     A(i*(m-2),j*(m-2)+1) = 0 ;  %Get Rid of Top Coefficients
%     A(i*(m-2)+1,j*(m-2)) = 0 ;  %Get Rid of Bottom Coefficients
% end
% A;

%_____________________________________________________________________
%%

%%Preallocate Temperatures and Residuals
T1 = zeros(N,1) ;
R = zeros(N,1) ;

%Setting Max Iteration Bounds
k_max = 1000 ;

%Set Error Bounds Based On Max Iterations
err=1/k_max;

%Calculating omega_opt for Reduced Iterations;
E = ((cos(pi/(i_max-1))+ratio^2*cos(pi/(j_max-1)))/(1+ratio^2))^2;
omega_opt = 2*((1-sqrt(1-E))/E);

for k = 1:k_max
    for i = 1:N
        sum1 = 0 ;
        sum2 = 0 ;
        for j = 1:(i-1)
            sum1 = sum1+A(i,j)*T1(j,k+1) ;
        end
        for j = i:N
            sum2 = sum2+A(i,j)*T1(j,k) ;
        end
        % Residuals
        R(i,k) = Q(i)-sum1-sum2 ;
        % New temperature, add a column for each iteration of T
        T1(i,k+1) = T1(i,k)+omega_opt*(R(i,k)/A(i,i)) ;
    end
    % Convergence criteria for |Delta_T|_max
    epsilon = max(abs(T1(:,k+1)-T1(:,k))) ;
    if epsilon < err
        break
    end
end

%_____________________________________________________________________
%%

% Print out the iteration history
T1 = T1(:,k) ;
iter = k ;
T(1:m,1) = left ;
T(1:m,n) = right ;
T(1,1:n) = bottom ;
T(m,1:n) = top ;
for i = 1:(n-2)
    T((2:(m-1)),i+1) = T1(((i-1)*(m-2)+1):i*(m-2)) ;
end
T ;

%_____________________________________________________________________
%%

%Plot Statements
figure(1)
pcolor(T)
shading('interp')
Top = num2str(top);
title({'2D Plate With Constant Temp Inputs';['T top =',num2str(top)];...
    ['T left =',num2str(left)];['T right =',num2str(right)];...
    ['T bottom =',num2str(bottom)]})
xlabel('X Grid Pt. #')
ylabel('Y Grid Pt. #')
axis image
hold on
contour(T,'k','LineWidth',2)