% =========================================================
%  贴体网格与几何变换
% =========================================================

clear; clc; close all;

% ---------------- 网格 ----------------
Nx = 80; Ny = 60;
xi  = linspace(0,1,Nx);
eta = linspace(0,1,Ny);
[XI, ETA] = meshgrid(xi, eta);

mode = 2;  % 1: 指数拉伸  % 2: 圆柱绕流


% 变换
switch mode

    case 1  % 指数拉伸（边界层）
    a = 3;
    X = XI;
    Y = (exp(a*ETA)-1)/(exp(a)-1);

    title_str = '指数拉伸（边界层网格）';

    case 2  % 圆柱绕流（极坐标）
    r_min = 1;
    r_max = 4;

    r = r_min + (r_max - r_min)*ETA;
    theta = 2*pi*XI;

    X = r .* cos(theta);
    Y = r .* sin(theta);

    title_str = '圆柱绕流网格（极坐标）';

end

% 可视化
figure('Color','w','Position',[100 100 1200 500])

subplot(1,2,1)
plot(XI, ETA, 'k'); hold on
plot(XI', ETA', 'k')
title('计算空间')
axis equal tight
grid on

subplot(1,2,2)
plot(X, Y, 'b'); hold on
plot(X', Y', 'b')
title(title_str)
axis equal tight
grid on

% 动态演化
figure('Color','w')

for t_anim = linspace(0,1,60)

    switch mode

        case 1
            a = 3*t_anim;
            Y = (exp(a*ETA)-1)/(exp(a+1e-6)-1);
            X = XI;

        case 2
            r = r_min + (r_max - r_min)*ETA;
            theta = 2*pi*XI;

            X = (1-t_anim)*XI + t_anim*(r.*cos(theta));
            Y = (1-t_anim)*ETA + t_anim*(r.*sin(theta));

    end

    clf
    plot(X, Y, 'r'); hold on
    plot(X', Y', 'r')

    title('网格生成过程（变换视角）')
    axis equal tight
    grid on

    pause(0.05)
end

%% =========================================================
% Jacobian
%% =========================================================
dx_dxi  = gradient(X, xi, eta);
dx_deta = gradient(X', eta, xi)';
dy_dxi  = gradient(Y, xi, eta);
dy_deta = gradient(Y', eta, xi)';

J = dx_dxi .* dy_deta - dx_deta .* dy_dxi;

figure('Color','w')
surf(X, Y, J, 'EdgeColor','none')
view(2); colorbar
title('Jacobian（网格质量/变形强度）')