clc
clear all

root = fileparts(mfilename('fullpath'));
addpath(root);
addpath(fullfile(root, '..', 'matlab 3D一阶基+散射边界条件'));
addpath(fullfile(root, '..', 'matlab 3D一阶基+散射边界条件+单周期边界'));

[phy, cfg] = case_config();

meshFile = fullfile(root, cfg.meshDatFile);
if ~isfile(meshFile)
    error(['找不到网格 %s。请先运行 getMesh.m 与 export_SBCmesh_to_dat.m。'], meshFile);
end

mesh = load_SBCmesh_dat(meshFile);
mesh.outIndex = findTri(phy.out, mesh);

solver.Ai = [];
solver.Aj = [];
solver.Av = [];
solver.b = zeros(mesh.NbrEdge, 1);

solver = assembly_equ(phy, mesh, solver);
solver = assembly_out(phy, mesh, solver);
solver = assembly_bele(phy, mesh, solver);
solver.A = sparse(solver.Ai, solver.Aj, solver.Av);
clear solver.Ai solver.Aj solver.Av

if isfield(phy, 'pec') && ~isempty(phy.pec)
    mesh.pecEdgeIds = find_pec_edges(phy.pec, mesh);
    [solver.A, solver.b] = apply_pec_bele(solver.A, solver.b, mesh, phy, mesh.pecEdgeIds);
    fprintf('已施加 PEC 边约束 %d 条\n', numel(mesh.pecEdgeIds));
end

outDir = fullfile(root, 'OutFile');
export_ab_matrix(solver.A, solver.b, outDir);

fprintf('求解中... condest(A)=%.4g\n', condest(solver.A));
solver.x = solver.A \ solver.b;

[Ex, Ey, Ez, normE] = get_ele_vertices_bele(mesh, solver, phy);
export_fields(Ex, Ey, Ez, normE, outDir);

fprintf('完成。OutFile: %s\n', outDir);
fprintf('  max(normE)=%.12g, mean(normE)=%.12g\n', max(normE), mean(normE));