function main_export() %MAIN_EXPORT 二阶单周期散射:导出 Ab 矩阵与 normE(格式同 3D 一阶散射 OutFile)。 % % 输出目录:./OutFile/ % Ai.txt, Aj.txt, Av.txt, Bv_real.txt, Bv_imag.txt — 投影后线性系统 P'*A*P, P'*b % Ex, Ey, Ez, normE — 求解后全场(每顶点一行) % % 另导出 ./OutFile_asm/:PBC 投影前的装配矩阵(与一阶散射 main.m 时机一致)。 % % 依赖:case_config.m 中 cfg.meshMatFile(默认 c6.mat)须在本目录。 clc; [phy, cfg] = case_config(); root = fileparts(mfilename('fullpath')); outDir = fullfile(root, 'OutFile'); outAsm = fullfile(root, 'OutFile_asm'); meshMat = fullfile(root, cfg.meshMatFile); if ~isfile(meshMat) error('main_export:NoMesh', '找不到网格 %s,请从 COMSOL 导出 c6.mat 到本目录。', meshMat); end load(meshMat); mesh.incIndex = findTri(phy.inc, mesh); mesh.outIndex = findTri(phy.out, mesh); [mesh.PBCIndex, mesh.origIndex] = findPBCIndex(phy.src, phy.dst, cfg.pbcAngle, mesh); mesh.PBCphi = cfg.pbcPhi; solver.Ai = []; solver.Aj = []; solver.Av = []; solver.dof = mesh.NbrEdge * 2 + mesh.NbrFace * 2; solver.b = zeros(solver.dof, 1); solver = assembly_equ(phy, mesh, solver); solver = assembly_out(phy, mesh, solver); solver = assembly_inc(phy, mesh, solver); solver.A = sparse(solver.Ai, solver.Aj, solver.Av); % 装配矩阵(无 PBC 投影),与一阶散射 export_ab_matrix 时机一致 export_ab_matrix(solver, outAsm); % Bloch 周期边界 + 投影 solver = assembly_rbc(mesh, solver); solver.A = solver.P' * solver.A * solver.P; solver.b = solver.P' * solver.b; % 投影后线性系统(实际求解的 Ab) export_ab_matrix(solver, outDir); % 求解(与 main_RBC.m 相同) r = symamd(solver.A); solver.A = solver.A(r, r); solver.b = solver.b(r); solver.x = solver.A \ solver.b; q = zeros(length(r), 1); for i = 1:length(r) q(r(i)) = i; end solver.x = solver.x(q); solver.x = solver.P * solver.x; [Ex, Ey, Ez, normE] = get_ele_vertices(mesh, solver); export_fields(Ex, Ey, Ez, normE, outDir); fprintf('\n完成。投影后系统已写入 %s\n', outDir); fprintf('装配矩阵(无 PBC)已写入 %s\n', outAsm); end