% Axisymmetric mug 
% r–z plane, rotation around r = 0
clear; clc; close all;

%% General values that we use in the entire script 
maxIter = 5000; 
tol = 1e-6;
%Task 3) Thermal conductivity 
lambda_wall  = 2.0;    % ceramic RANDOM NUMBER
lambda_fluid = 0.6;    % water RANDOM NUMBER
lambda_air   = 0.025;  % air RANDOM NUMBER
%Task 4) 
alpha = 10;      % heat transfer coefficient
u_out = 18;      % ambient air temperature (°C)
%Task 6) Heat capacities (example values)
c_wall  = 900;   % ceramic RANDOM NUMBER
c_fluid = 4180;  % water RANDOM NUMBER
c_air   = 1005;  % air RANDOM NUMBER

% This numbers don't simulate real life because in real life we have
% volumetric heat capacities. this means we are having a reduced factore of
% approx 10^4. 

%% Task 1: Mesh definition
%  Create PDE model
model = createpde();

%  Points (meters)
A = [0,     0];
B = [0.055, 0];
C = [0.083, 0.105];
H = [0.078, 0.105];
F = [0.050, 0.005];
E = [0,     0.005];
G = [0.067, 0.066];
I = [0,     0.066];
D = [0,     0.105];

%  Geometry matrix (edges)
% Axis (split by material)
g1  = [2; A(1); E(1); A(2); E(2); 1; 0]; % ceramic
g2  = [2; E(1); I(1); E(2); I(2); 2; 0]; % fluid
g3  = [2; I(1); D(1); I(2); D(2); 3; 0]; % air 
% Outer ceramic
g4  = [2; A(1); B(1); A(2); B(2); 1; 0];
g5  = [2; B(1); C(1); B(2); C(2); 1; 0];
% Top rim: C -> H is ceramic-air (ONLY)
g6  = [2; C(1); H(1); C(2); H(2); 1; 3];
% Inner ceramic wall (H -> F) ceramic-air
g7  = [2; H(1); F(1); H(2); F(2); 1; 3];
% Inner ceramic bottom (F -> E) ceramic-fluid
g8  = [2; F(1); E(1); F(2); E(2); 1; 2];
% Fluid surface (F -> G) fluid-air
g9  = [2; F(1); G(1); F(2); G(2); 2; 3];
% Fluid surface (G -> I) fluid-air
g10 = [2; G(1); I(1); G(2); I(2); 2; 3];
% Air top boundary: D -> H (air-outside)
g11 = [2; D(1); H(1); D(2); H(2); 3; 0];

% Assemble geometry
g = [g1 g2 g3 g4 g5 g6 g7 g8 g9 g10 g11];

%  Load geometry into PDE model
geometryFromEdges(model, g);

%  Plot geometry with labels
% figure(1);
% pdegplot(model, 'EdgeLabels','on', 'FaceLabels','on');
% axis equal;
% title('Geometry with edge and face labels');

%  Generate mesh
mesh = generateMesh(model, 'Hmax', 0.006, 'GeometricOrder','linear');

%  Plot mesh
% figure(2);
% pdemesh(model);
% axis equal;
% title('Generated mesh');

%% Task 2: Jacobi solver with constant lambda
nodes = mesh.Nodes;         % coordinates of all mesh nodes
elements = mesh.Elements;   % which nodes make up each triangle element

Nnodes = size(nodes,2);
Nelems = size(elements,2);

% Define material properties (for simplicity, lambda = 1 everywhere)
lambda = ones(Nelems,1);    % thermal conductivity

% Initialize global stiffness matrix and RHS
K = sparse(Nnodes, Nnodes);
F = zeros(Nnodes,1);

% Assemble K and F
for e = 1:Nelems %Loop over each triangle element 
    vert = elements(:,e);   %nodes of element
    x = nodes(1,vert);
    y = nodes(2,vert);

    Ae = polyarea(x,y);     % Compute area of the triangle 

    % Linear triangle gradients
    b = [y(2)-y(3); y(3)-y(1); y(1)-y(2)]; % derivative with respect to x 
    c = [x(3)-x(2); x(1)-x(3); x(2)-x(1)]; % derivative with respect to y

    % Element stiffness matrix
    Ke = (lambda(e)/(4*Ae)) * (b*b.' + c*c.');

    % Assemble
    K(vert,vert) = K(vert,vert) + Ke;

    % Element load vector (f=0)
    F(vert) = F(vert) + zeros(3,1);
end

% Find boundary nodes (couldn't find a better way)
edgesAll = [elements([1 2],:), elements([2 3],:), elements([3 1],:)]; % all edges -> 2x(3*Nelems) awway since each triangle has 3 edges 
edgesSorted = sort(edgesAll,1); % sort nodes of each edge, ensure that [i,j] and [j,i] are the same
[~,~,ic] = unique(edgesSorted','rows'); % identifies unique edges and assignes indices ic 
counts = accumarray(ic,1); % counts how many times each edge appears in the mesh
boundaryEdges = find(counts==1); %edges belonging to only 1 element: hence on the boundary 
boundaryNodes = unique(edgesSorted(:,boundaryEdges)); % nodes belonging to these boundary edges 

% Jacobi iteration
D = diag(K);
R = K - diag(D);
u = zeros(Nnodes,1);

for iter = 1:maxIter
    u_new = (F - R*u) ./ D;
    u_new(boundaryNodes) = 0; % Enforce Dirichlet BC
    if norm(u_new - u, inf) < tol % Convergence check
        fprintf('Task 2) Jacobi converged in %d iterations.\n', iter);
        break;
    end
    u = u_new;
end

% % Plot solution
% figure(3)
% pdeplot(model, 'XYData', u, 'Mesh','on');
% axis equal;
% title('Stationary Dirichlet solution via Jacobi');
% colorbar;

%% Task 3: Laplace with multiple lambdas
[K, F] = CalculateLaplace_mult(model, lambda_wall, lambda_fluid, lambda_air);

D = diag(K);
R = K - diag(D);
u = zeros(Nnodes,1);
for iter = 1:maxIter
    u_new = (F - R*u) ./ D;
    u_new(boundaryNodes) = 0;   % Dirichlet BC

    if norm(u_new - u, inf) < tol
        fprintf('Task 3) Jacobi converged in %d iterations.\n', iter);
        break;
    end
    u = u_new;
end

% figure(4)
% pdeplot(model, 'XYData', u, 'Mesh','on');
% axis equal
% title('Task 3: Stationary solution with multiple conductivities');
% colorbar

%% Task 4: Robin boundary condition 
[K, F] = ApplyRobinBC_mult(model, K, F, alpha, u_out);
D = diag(K);
R = K - diag(D);
u = zeros(Nnodes,1);

for iter = 1:maxIter
    u_new = (F - R*u) ./ D;

    if norm(u_new - u, inf) < tol
        fprintf('Task 4) Jacobi (Robin BC) converged in %d iterations.\n', iter);
        break;
    end
    u = u_new;
end

figure(5)
pdeplot(model, 'XYData', u, 'Mesh','on');
axis equal
title('Task 4: Stationary solution with Robin BC');
colorbar

%% Task 5: Axisymmetric Laplace + Robin BC 
[K, F] = CalculateLaplace_mult_rot(model, lambda_wall, lambda_fluid, lambda_air);
[K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out);

D = diag(K);
R = K - diag(D);
u = zeros(Nnodes,1);

for iter = 1:maxIter
    u_new = (F - R*u) ./ D;
    if norm(u_new - u, inf) < tol
        fprintf('Task 5) Jacobi (axisymmetric Robin) converged in %d iterations.\n', iter);
        break;
    end
    u = u_new;
end

% Plot solution
% figure(6)
% pdeplot(model, 'XYData', u, 'Mesh','on');
% axis equal
% title('Task 5: Axisymmetric stationary solution with Robin BC');
% colorbar

%To see it in 3D paste here the code in "AdditionalPlotCodes.txt". 

%% Task 6: Axisymmetric mass matrix
M = sparse(Nnodes, Nnodes);
M = AddMass_mult_rot(model, M, c_wall, c_fluid, c_air);

%% Task 7: Initial solution
u0 = Init_Solution_mult(model, 18, 80, 18);

% figure(7)
% pdeplot(model, 'XYData', u0, 'Mesh','on');
% axis equal
% title('Initial temperature distribution');
% colorbar

%% Task 8: Time-dependent simulation (explicit scheme)
tau = 0.5;              % time step in seconds 
T_end = 400;            % total simulation time (seconds)
Nt = ceil(T_end/tau);   % number of time steps

A = (1/tau)*M+K; % Left-hand side matrix
[L,U,P,Q] = lu(A); % A is constant -> factorize it once: PAQ=LU

% Initialize solution
u = u0;

for k = 1:Nt
    b = (1/tau)*M*u + F;   % F is the load vector, F=0 
    
    % Solve for next time step
    u_next = Q*(U\(L\(P*b))); % efficient solution using LU factors 
    
    u = u_next; % Update
    
    if mod(k,20) == 0
        % figure(8)
        % pdeplot(model, 'XYData', u, 'Mesh','on');
        % axis equal
        % title(['Temperature at t = ', num2str(k*tau), ' s']);
        % colorbar
        % drawnow
    end
end

%To see the 9 snapshots paste here the codes in "AdditionalPlotCodes.txt"

%% Task 9 (i): Heating time using inner ceramic wall temperature
T_target = 67;   % [°C]

innerWallNodes = findNodes(model.Mesh,'region','Edge',8); % Edge 8 = ceramic–fluid interface 
u = u0;

% Storage
timeVec = (0:Nt-1)' * tau;
innerWallTemp = zeros(Nt,1);
Twarm = NaN;

for k = 1:Nt
    b = (1/tau)*M*u + F;

    u = Q*(U\(L\(P*b)));

    % Average inner wall temperature
    innerWallTemp(k) = mean(u(innerWallNodes));

    % Check heating criterion
    if innerWallTemp(k) >= T_target
        Twarm = k * tau;
        fprintf('Task 9 (i): Inner wall reaches %.1f°C at T = %.1f s\n', ...
                T_target, Twarm);
        break
    end
end

%% --- Plot inner wall temperature evolution
figure(9)
plot(timeVec(1:k), innerWallTemp(1:k), 'LineWidth', 2)
hold on
yline(T_target,'r--','67°C','LineWidth',1.5)
xlabel('Time [s]')
ylabel('Average inner wall temperature [°C]')
title('Heating of the inner ceramic wall')
grid on