47 lines
1.4 KiB
Matlab
47 lines
1.4 KiB
Matlab
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));
|