122 lines
3.4 KiB
Matlab
122 lines
3.4 KiB
Matlab
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
|