106 lines
2.3 KiB
Text
106 lines
2.3 KiB
Text
%% TASK 5: 3D plot and video
|
|
% Nodes and solution from Task 5
|
|
r = nodes(1,:)'; % r-coordinate
|
|
z = nodes(2,:)'; % z-coordinate
|
|
u_sol = u; % solution vector
|
|
|
|
% Create a fine phi grid for rotation
|
|
phi = linspace(0, 2*pi, 50);
|
|
|
|
% Create 3D coordinates
|
|
[Phi, R] = meshgrid(phi, r);
|
|
Z = repmat(z, 1, length(phi));
|
|
X = R .* cos(Phi);
|
|
Y = R .* sin(Phi);
|
|
|
|
% Interpolate u onto this grid (simple nearest neighbor)
|
|
U = repmat(u_sol, 1, length(phi));
|
|
|
|
% Plot
|
|
figure
|
|
surf(X, Y, Z, U, 'EdgeColor','none');
|
|
axis equal
|
|
xlabel('x'); ylabel('y'); zlabel('z');
|
|
title('Axisymmetric 3D solution');
|
|
colorbar
|
|
view(30,30)
|
|
|
|
% Assume you already have the 3D plot (X,Y,Z,U)
|
|
figure
|
|
h = surf(X,Y,Z,U, 'EdgeColor','none');
|
|
axis equal
|
|
colorbar
|
|
xlabel('x'); ylabel('y'); zlabel('z');
|
|
title('Axisymmetric 3D solution');
|
|
|
|
v = VideoWriter('mug_rotation.avi'); % create video file
|
|
open(v);
|
|
|
|
for angle = 0:2:360
|
|
view(angle, 30); % rotate camera around z-axis
|
|
drawnow;
|
|
frame = getframe(gcf); % capture current frame
|
|
writeVideo(v, frame);
|
|
end
|
|
close(v);
|
|
|
|
|
|
%% TASK 8: 9 snapshot of the solution over time:
|
|
|
|
% --- Snapshot times (seconds)
|
|
snapshot_times = [0 50 100 150 200 250 300 350 400];
|
|
snapshot_steps = round(snapshot_times / tau) + 1;
|
|
|
|
Nsnap = length(snapshot_steps);
|
|
snapshots = zeros(Nnodes, Nsnap);
|
|
|
|
% Initial condition (t = 0)
|
|
snapshots(:,1) = u0;
|
|
snapshot_counter = 2;
|
|
|
|
% --- Time stepping parameters
|
|
T_end = snapshot_times(end);
|
|
Nt = ceil(T_end / tau);
|
|
|
|
% --- System matrix
|
|
A = (1/tau)*M + K;
|
|
|
|
% LU factorization (constant matrix)
|
|
[L,U,P,Q] = lu(A);
|
|
|
|
% --- Initialize solution
|
|
u = u0;
|
|
|
|
% --- Time stepping loop
|
|
for k = 1:Nt
|
|
% Right-hand side
|
|
b = (1/tau)*M*u + F;
|
|
|
|
% Solve system
|
|
u_next = Q * (U \ (L \ (P * b)));
|
|
u = u_next;
|
|
|
|
% Store snapshots at requested times
|
|
if snapshot_counter <= Nsnap
|
|
if k+1 == snapshot_steps(snapshot_counter)
|
|
snapshots(:,snapshot_counter) = u;
|
|
snapshot_counter = snapshot_counter + 1;
|
|
end
|
|
end
|
|
end
|
|
|
|
%% --- 3x3 grid plot (overview figure)
|
|
figure(20)
|
|
for i = 1:Nsnap
|
|
subplot(3,3,i)
|
|
pdeplot(model, 'XYData', snapshots(:,i), 'Mesh','off');
|
|
axis equal
|
|
title(['T = ', num2str(snapshot_times(i)), ' s'])
|
|
end
|
|
colorbar
|
|
sgtitle('Transient temperature evolution')
|
|
|
|
% Save combined figure
|
|
saveas(gcf, 'Task8_Transient_3x3.png');
|
|
|
|
|