% 将 SBCmesh.mat 导出为 OpticsFEM3D 可读文本网格 SBCmesh.dat
% 用法：先运行 getMesh.m 生成 SBCmesh.mat，再运行本脚本

clc;
if ~isfile('MPDmesh.mat')
    error('请先运行 getMesh.m 生成 MPDmesh.mat');
end
load('MPDmesh.mat', 'mesh');

fid = fopen('MPDmesh.dat', 'w');
fprintf(fid, 'NbrVertex\n%d\nVertex\n', mesh.NbrVertex);
% 兼容 Vertex 为 3×N 或 N×3
V = mesh.Vertex;
if size(V,1) == 3 && size(V,2) ~= 3
    V = V';
end
for i = 1:mesh.NbrVertex
    fprintf(fid, '%.16g %.16g %.16g\n', V(i,1), V(i,2), V(i,3));
end

T = mesh.Tet;
if size(T,1) == 4 && size(T,2) ~= 4
    T = T';
end
fprintf(fid, 'NbrTet\n%d\nTet\n', mesh.NbrTet);
for i = 1:mesh.NbrTet
    fprintf(fid, '%d %d %d %d\n', T(i,1), T(i,2), T(i,3), T(i,4));
end

fprintf(fid, 'DomainOfTet\n');
for i = 1:mesh.NbrTet
    fprintf(fid, '%d\n', mesh.DomainOfTet(i));
end

fprintf(fid, 'NbrEdge\n%d\nEdge\n', mesh.NbrEdge);
for i = 1:mesh.NbrEdge
    fprintf(fid, '%d %d\n', mesh.Edge(i,1), mesh.Edge(i,2));
end

EOT = mesh.EdgeOfTet;
if size(EOT,1) == 6
    EOT = EOT';
end
fprintf(fid, 'EdgeOfTet\n');
for i = 1:mesh.NbrTet
    fprintf(fid, '%d %d %d %d %d %d\n', EOT(i,1), EOT(i,2), EOT(i,3), EOT(i,4), EOT(i,5), EOT(i,6));
end

Tr = mesh.Tri;
if size(Tr,1) == 3
    Tr = Tr';
end
fprintf(fid, 'NbrTri\n%d\nTri\n', mesh.NbrTri);
for i = 1:mesh.NbrTri
    fprintf(fid, '%d %d %d\n', Tr(i,1), Tr(i,2), Tr(i,3));
end

fprintf(fid, 'DomainOfTri\n');
for i = 1:mesh.NbrTri
    fprintf(fid, '%d\n', mesh.DomainOfTri(i));
end

fprintf(fid, 'ConnOfTri\n');
for i = 1:mesh.NbrTri
    fprintf(fid, '%d %d\n', mesh.ConnOfTri(i,1), mesh.ConnOfTri(i,2));
end

% 面法向（与 getMesh.m 中 NormOfFace 一致）
normRows = [1;2;3;4;5;14];
fprintf(fid, 'NormOfFace\n%d\n', numel(normRows));
for k = 1:length(normRows)
    d = normRows(k);
    n = mesh.NormOfFace(d,:);
    fprintf(fid, '%d %.16g %.16g %.16g\n', d, n(1), n(2), n(3));
end

fclose(fid);
fprintf('已写出 MPDmesh.dat\n');