function mesh = load_SBCmesh_dat(datFile)
%LOAD_SBCMESH_DAT  读取 OpticsFEM/COMSOL 导出的 SBCmesh.dat（1-based 索引）。
%  与 3D opticsfem-master Mesh_3D 读入格式一致，供 main.m 与 C++ 共用同一网格。

fid = fopen(datFile, 'r');
if fid < 0
    error('无法打开网格文件: %s', datFile);
end
cleanup = onCleanup(@() fclose(fid));

tag = read_tag(fid);
assert(strcmp(tag, 'NbrVertex'), 'Expected NbrVertex');
mesh.NbrVertex = read_int(fid);

tag = read_tag(fid);
assert(strcmp(tag, 'Vertex'), 'Expected Vertex');
mesh.Vertex = reshape(fscanf(fid, '%f', mesh.NbrVertex * 3), 3, []).';

tag = read_tag(fid);
assert(strcmp(tag, 'NbrTet'), 'Expected NbrTet');
mesh.NbrTet = read_int(fid);

tag = read_tag(fid);
assert(strcmp(tag, 'Tet'), 'Expected Tet');
mesh.Tet = reshape(fscanf(fid, '%d', mesh.NbrTet * 4), 4, []).';

tag = read_tag(fid);
assert(strcmp(tag, 'DomainOfTet'), 'Expected DomainOfTet');
mesh.DomainOfTet = fscanf(fid, '%d', mesh.NbrTet);

tag = read_tag(fid);
assert(strcmp(tag, 'NbrEdge'), 'Expected NbrEdge');
mesh.NbrEdge = read_int(fid);

tag = read_tag(fid);
assert(strcmp(tag, 'Edge'), 'Expected Edge');
mesh.Edge = reshape(fscanf(fid, '%d', mesh.NbrEdge * 2), 2, []).';

tag = read_tag(fid);
assert(strcmp(tag, 'EdgeOfTet'), 'Expected EdgeOfTet');
mesh.EdgeOfTet = reshape(fscanf(fid, '%d', mesh.NbrTet * 6), 6, []).';

tag = read_tag(fid);
assert(strcmp(tag, 'NbrTri'), 'Expected NbrTri');
mesh.NbrTri = read_int(fid);

tag = read_tag(fid);
assert(strcmp(tag, 'Tri'), 'Expected Tri');
mesh.Tri = reshape(fscanf(fid, '%d', mesh.NbrTri * 3), 3, []).';

tag = read_tag(fid);
assert(strcmp(tag, 'DomainOfTri'), 'Expected DomainOfTri');
mesh.DomainOfTri = fscanf(fid, '%d', mesh.NbrTri);

tag = read_tag(fid);
assert(strcmp(tag, 'ConnOfTri'), 'Expected ConnOfTri');
mesh.ConnOfTri = reshape(fscanf(fid, '%d', mesh.NbrTri * 2), 2, []).';

tag = read_tag(fid);
assert(strcmp(tag, 'NormOfFace'), 'Expected NormOfFace');
nbrNorm = read_int(fid);
maxDom = 14;
mesh.NormOfFace = zeros(maxDom, 3);
for k = 1:nbrNorm
    row = fscanf(fid, '%f', 4);
    if numel(row) < 4
        error('NormOfFace 第 %d 条记录不完整（文件可能被截断）', k);
    end
    domainId = round(row(1));
    if domainId > maxDom
        extra = zeros(domainId, 3);
        extra(1:maxDom, :) = mesh.NormOfFace;
        mesh.NormOfFace = extra;
        maxDom = domainId;
    end
    mesh.NormOfFace(domainId, :) = row(2:4);
end
% 仅对仍为 0 的 SBC 模板域补默认（PBC 网格 1..9 已在 dat 中给出）
defaults = {
    1, [-1, 0, 0];
    2, [0, -1, 0];
    3, [0, 0, -1];
    4, [0, 0, 1];
    5, [0, 1, 0];
    8, [0, 0, 1];
    14, [1, 0, 0];
};
for k = 1:size(defaults, 1)
    d = defaults{k, 1};
    if norm(mesh.NormOfFace(d, :)) < 1e-30
        mesh.NormOfFace(d, :) = defaults{k, 2};
    end
end

if min(mesh.EdgeOfTet(:)) < 1
    error(['SBCmesh.dat 中 EdgeOfTet 必须为 1-based（最小值=%d）。' ...
        '请使用 export_SBCmesh_to_dat 或修正后再读。'], min(mesh.EdgeOfTet(:)));
end

fprintf('已加载 %s: NbrTet=%d, NbrEdge=%d, NbrTri=%d\n', ...
    datFile, mesh.NbrTet, mesh.NbrEdge, mesh.NbrTri);
end

function tag = read_tag(fid)
tag = '';
while true
    line = fgetl(fid);
    if ~ischar(line)
        error('Unexpected end of mesh file.');
    end
    line = strtrim(line);
    if ~isempty(line)
        tag = line;
        return;
    end
end
end

function value = read_int(fid)
value = fscanf(fid, '%d', 1);
end