非構造格子の有限体積法で２次元拡散方程式を解く【MATLABコード付き】

clearvars;

%% パラメータ
D = 0.01;
dt = 1e-3;
t_end = 5.0;
steps = ceil(t_end / dt);
frame_rate = 20;  % 任意のフレームレートに変更可能

%% メッシュ生成
model = createpde();
geometryFromEdges(model,@squareg);
mesh = generateMesh(model,'Hmax',0.05);
p = mesh.Nodes;
t = mesh.Elements(1:3,:);%　要素を構成する3つのノード番号
[Np, Ne] = deal(size(p,2), size(t,2));

%% メッシュ可視化
figure(1);
pdemesh(model);
% 背景を白に
fig=gcf;
fig.Color='white';

% 軸ラベルやフォント設定
xlabel('$$x$$','Interpreter','latex','FontSize',12);
ylabel('$$y$$','Interpreter','latex','FontSize',12);
set(gca, 'FontName','Times New Roman', 'FontSize',12);

% 軸のアスペクト比・表示範囲の設定
axis('equal');
axis('tight');

% 'TickLabelFormat'を利用してX軸の桁数を指定
ax = gca;
ax.XAxis.TickLabelFormat = '%.1f'; % X軸の桁数を小数点以下2桁に設定
ax.YAxis.TickLabelFormat = '%.1f'; % X軸の桁数を小数点以下2桁に設定
drawnow;

%% 要素重心と面積
centroids = (p(:,t(1,:)) + p(:,t(2,:)) + p(:,t(3,:))) / 3;
area = zeros(Ne,1);
for e = 1:Ne
    coords = p(:,t(:,e));% 要素eの3つのノードの座標
    area(e) = polyarea(coords(1,:), coords(2,:));
end

%% 初期条件（ある領域の中 = 1, それ以外 = 0）
x_min = -0.2; x_max = 0.2;
y_min = -0.2; y_max = 0.2;
u = double(centroids(1,:) >= x_min & centroids(1,:) <= x_max & ...
    centroids(2,:) >= y_min & centroids(2,:) <= y_max)';

%% エッジ情報と面法線ベクトル
[edges, edge2elems, Sf_vec] = build_edges_and_normals(t, p, centroids);

%% ステンシル係数行列 A
A = sparse(Ne, Ne);
for f = 1:size(edges,1)% エッジのループ
    e1 = edge2elems(f,1);% エッジに接する要素１
    e2 = edge2elems(f,2);% エッジに接する要素２
    if e2 == 0, continue; end% エッジに接している要素が一つしかない場合、計算をスキップする。フラックス0となりノイマン条件を満たす。
    d_vec = centroids(:,e2) - centroids(:,e1);
    len2 = norm(d_vec)^2;
    k = D * dot(Sf_vec(:,f), d_vec) / len2;
    A(e1,e1) = A(e1,e1) + k;
    A(e2,e2) = A(e2,e2) + k;
    A(e1,e2) = A(e1,e2) - k;
    A(e2,e1) = A(e2,e1) - k;
end

%% 時間積分行列系: (M + dt*A) * u^{n+1} = M * u^n
M = spdiags(area, 0, Ne, Ne);
LHS = M + dt * A;

%% 可視化用（ノード補間）
W = sparse(Np, Ne);
for e = 1:Ne
    W(t(:,e), e) = W(t(:,e), e) + 1;
end
node_weights = sum(W,2); node_weights(node_weights==0) = 1;

%% MP4動画保存の初期化
video_filename = 'diffusion_result.mp4';
v = VideoWriter(video_filename, 'MPEG-4');
v.FrameRate = frame_rate;
open(v);

%% 時間発展（陰解法）
draw_interval = 10;
for step = 1:steps
    RHS = M * u;
    u = LHS \ RHS;

    % 可視化（補間 + clip）
    if mod(step, draw_interval)==0 || step==steps
        figure(2);
        u_node = (W * u) ./ node_weights;
        u_node = max(u_node, 0);
        trisurf(t', p(1,:), p(2,:), u_node(:), 'EdgeColor', 'none');
        view(2); colorbar; caxis([0 1]);
        title(sprintf('t = %.3f', step*dt));

        % 背景を白に
        fig=gcf;
        fig.Color='white';

        % 軸ラベルやフォント設定
        xlabel('$$x$$','Interpreter','latex','FontSize',12);
        ylabel('$$y$$','Interpreter','latex','FontSize',12);
        set(gca, 'FontName','Times New Roman', 'FontSize',12);

        % 軸のアスペクト比・表示範囲の設定
        axis('equal');
        axis('tight');

        % 'TickLabelFormat'を利用してX軸の桁数を指定
        ax = gca;
        ax.XAxis.TickLabelFormat = '%.1f'; % X軸の桁数を小数点以下2桁に設定
        ax.YAxis.TickLabelFormat = '%.1f'; % X軸の桁数を小数点以下2桁に設定
        drawnow;

        frame = getframe(gcf);
        writeVideo(v, frame);
    end
end

close(v);
fprintf('✅ 動画保存完了: %s\n', video_filename);


function [edges, edge2elems, Sf_vec] = build_edges_and_normals(t, p, centroids)
edge_map = containers.Map('KeyType','char','ValueType','int32');
edges = []; edge2elems_map = containers.Map; edge_count = 0;
Ne = size(t,2);

for e = 1:Ne
    nodes = t(:,e);
    local_edges = [ ...
        [nodes(1); nodes(2)], ...
        [nodes(2); nodes(3)], ...
        [nodes(3); nodes(1)] ...
        ];

    for k = 1:3
        edge = sort(local_edges(:,k));
        key = sprintf('%d-%d', edge(1), edge(2));
        if ~isKey(edge_map, key)
            edge_count = edge_count + 1;
            edge_map(key) = edge_count;
            edges(edge_count,:) = edge;
            edge2elems_map(key) = [e, 0];
        else
            pair = edge2elems_map(key); pair(2) = e;
            edge2elems_map(key) = pair;
        end
    end
end

Nedges = size(edges,1);
edge2elems = zeros(Nedges,2); Sf_vec = zeros(2,Nedges);
keys = edge_map.keys;
for i = 1:Nedges
    key = keys{i}; edge = sscanf(key, '%d-%d');
    idx = edge_map(key);
    elems = edge2elems_map(key);
    edge2elems(idx,:) = elems;
    v = p(:,edge(2)) - p(:,edge(1));
    normal = [0 -1; 1 0] * v / norm(v);
    Sf = norm(v) * normal;
    if elems(2) ~= 0
        c1 = centroids(:,elems(1)); c2 = centroids(:,elems(2));
        if dot(Sf, c2 - c1) < 0, Sf = -Sf; end
    end
    Sf_vec(:,idx) = Sf;
end
end