% ------------------------------------------------------------------------
% 
% Problem: 2D steady-state conduction, discontinuous conducibility
%
% Domain: square [0,1] x [0,1]
%
% Grid: structured, cartesian, constant spacing, dx = dy
%
% Boundary conditions (BC):
%   Dirichet ( T           = T_BC ) on the lower side ( y = 0 )
%   Neumann  ( Gamma * T_n = Jd_BC  ) on the other sides
%
% Discretization method: Finite Volume (FV)
%
% Solution methodology: direct (sparse coefficient matrix)
%
% Date: 10/10/2024
% Author: R. Zamolo
%
% ------------------------------------------------------------------------
%
% Domain, coordinates systems and boundary conditions:
%
%         ^ y( j ) 
%         |
%         |
%         |    dT/dy=0
%       1 +---------------+
%         |               |
%         |               |
%      T1 |               | T2
%         |               |
%         |               |
%       0 +---------------+------> x( i )
%         0    dT/dy=0    1
%
% Spatial indexing: T( i , j )
%   i: index along x: i = 1 , ... , Ni
%   j: index along y: j = 1 , ... , Ni
%
% Linear indexing: T( k )
%   k = i + (j-1) * Ni
%
% The number of unknowns is N = Ni * Ni
% 
%   k_max = i_max + (j_max-1) * Ni =
%           Ni    + (Ni   -1) * Ni = Ni * Ni = N
%
% ------------------------------------------------------------------------
% Properties
GammaA = 10 ;
GammaB = 1 ;
T1 = 1 ;
T2 = 0 ;

% Domain square side length
L = 2 ;

% FV number for each dimension and FV total number
Ni = 600 ;
N  = Ni * Ni ;

% FV side length
dx = L / Ni ;

% FV 'volume' (surface)
dA = dx * dx ;

% FV centroids coordinates
X = dx * ( (1:Ni) - 0.5 )' ;
Y = X ;

% FV equation:
%     A_P*T_P + A_E*T_E + A_W*T_W + A_N*T_N + A_S*T_S = S_P
%
% Final system, in compact (matrix) form:
%                           A * T = S
%
%   A: sparse coefficient matrix
% T: solution vector
%   S: source term vector (RHS)

% Diagonals for Gamma (discontinuous)
D = ones( Ni , Ni ) ;
D_Gamma = GammaA * D ;
h = L/2 ;
D_Gamma(X>h, Y<h) = GammaB ;

% Preparation of the 5 diagonals of A, now stored as Ni x Ni matrices
% (spatial indexing)
iM = [ 1 1:(Ni-1) ] ;
iP = [ 2:Ni Ni ] ;
D_W = 2 ./ ( 1 ./ D_Gamma + 1 ./ D_Gamma(iM,:) ) ;
D_E = 2 ./ ( 1 ./ D_Gamma + 1 ./ D_Gamma(iP,:) ) ;
D_N = 2 ./ ( 1 ./ D_Gamma + 1 ./ D_Gamma(:,iP) ) ;
D_S = 2 ./ ( 1 ./ D_Gamma + 1 ./ D_Gamma(:,iM) ) ;
D_P = -(D_W + D_E + D_N + D_S) ;

% RHS preparation (midpoint second order integration)
% stored as Ni x Ni matrix (spatial indexing)
S = 0*D ;

% Diagonals correction in order to impose BC
%   South side ( y=0, j=1 ) (Neumann)
D_P( : , 1 ) = D_P( : , 1 ) + D_S( : , 1 ) ;
%  S( : , 1 ) =   S( : , 1 ) - D_S( : , 1 ) * dx * Jd_BC / Gamma ;
D_S( : , 1 ) = 0 ;

%   North side ( y=L, j=Ni ) (Neumann)
D_P( : , Ni ) = D_P( : , Ni ) + D_N( : , Ni ) ;
%  S( : , Ni ) =   S( : , Ni ) - D_N( : , Ni ) * dx * Jd_BC / Gamma ;
D_N( : , Ni ) = 0 ;

%   West side ( x=0, i=1 ) (Dirichlet)
D_P( 1 , : ) = D_P( 1 , : ) -   D_W( 1 , : ) ;
  S( 1 , : ) =   S( 1 , : ) - 2*D_W( 1 , : ) * T1 ;
D_W( 1 , : ) = 0 ;

%   East side ( x=L, i=Ni ) (Dirichlet)
D_P( Ni , : ) = D_P( Ni , : ) -   D_E( Ni , : ) ;
  S( Ni , : ) =   S( Ni , : ) - 2*D_E( Ni , : ) * T2 ;
D_E( Ni , : ) = 0 ;

% Diagonals (and RHS) reshape to get vectors (linear indexing)
D_P = reshape( D_P , N , 1 ) ;
D_W = reshape( D_W , N , 1 ) ;
D_E = reshape( D_E , N , 1 ) ;
D_N = reshape( D_N , N , 1 ) ;
D_S = reshape( D_S , N , 1 ) ;
S   = reshape(   S , N , 1 ) ;

% Lower diagonals (South/West) and upper diagonals (East/North) shift.
% This shift operation is needed because of the way the following spdiags()
% command operates.
% (The circular property is not strictly needed here...)
D_S = circshift( D_S , -Ni ) ;
D_W = circshift( D_W ,  -1 ) ;
D_E = circshift( D_E ,   1 ) ;
D_N = circshift( D_N ,  Ni ) ;

% Coefficient matrix A creation from diagonals (sparse form)
A = spdiags( [ D_S D_W D_P D_E D_N ] , [ -Ni -1 0 1 Ni ] , N , N ) ;

% Direct solution (\ backslash)
T = A\S ;

% Solution reshape to get again the spatial indexing
T = reshape( T , [ Ni , Ni ] ) ;

flux = sum( 2 * D_Gamma(Ni,:) .* (T(Ni,:) - T2) ) ;
Rtot_FV = ( T1 - T2 ) / flux ;
fprintf('Rtot_FV     = %.6f m^2K/W (valore "esatto")\n', Rtot_FV) ;

% Reti elettriche equivalenti
R  = @(width, H, conductivity) width / (H*conductivity) ;
Rp = @(R1, R2) 1 / ( 1/R1 + 1/R2 ) ;

% Parallelo (superficie orizzontale adiabatica)
halfL = L/2 ;
Rtot_parallel = Rp( R(L,halfL,GammaA) , R(halfL,halfL,GammaA) + R(halfL,halfL,GammaB) ) ;
fprintf('Rtot_parall = %.4f   m^2K/W\n', Rtot_parallel) ;

% Serie (superficie verticale isoterma)
Rtot_serial = R(halfL,L,GammaA) + Rp( R(halfL,halfL,GammaA) , R(halfL,halfL,GammaB) ) ;
fprintf('Rtot_serial = %.4f   m^2K/W\n', Rtot_serial) ;

% Media
R_avg = ( Rtot_parallel + Rtot_serial ) / 2 ;
fprintf('Rtot_avg    = %.4f   m^2K/W\n', R_avg) ;

% Medie alternative
C_parallel = 1 / Rtot_parallel ;
C_serial   = 1 / Rtot_serial ;
C_avg = ( C_parallel + C_serial ) / 2 ;
C_avgR = 1 / R_avg ;
C_avgavg = ( C_avg + C_avgR ) / 2 ;
R_avgavg = 1 / C_avgavg ;
fprintf('Rtot_avgavg = %.6f m^2K/W\n', R_avgavg) ;

% Contour line plot (orthogonal to the diffusive flux) on figure 1
figure( 1 ) ;
contourf( X , Y , T', 20 ) ;

% Naming axes
xlabel( 'x' ) ;
ylabel( 'y' ) ;

ax1 = gca ;
ax1.XLim = [0 L] ;
ax1.YLim = [0 L] ;

return

% Surface plot on figure 2
figure( 2 ) ;
Superf = surf( X , Y , T' ) ;

% Naming axes
xlabel( 'x' ) ;
ylabel( 'y' ) ;
zlabel( 'T' ) ;

% Pretty surface properties
Superf.EdgeColor    = 'none'    ; % FV sides not drawn
Superf.FaceLighting = 'gouraud' ; % Lighting model
Superf.FaceColor    = 'interp'  ; % Color interpolation inside FV

% Adding a light
Light_1          = light ;
Light_1.Position = [ -2 0 2 ] ;
Light_1.Color    = [ .7 .7 .7 ] ;

ax2 = gca ;
ax2.XLim = [0 L] ;
ax2.YLim = [0 L] ;