37 lines
1.5 KiB
Matlab
37 lines
1.5 KiB
Matlab
[phy, cfg] = case_config();
|
|
load(cfg.meshMatFile);
|
|
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);
|
|
solver = assembly_rbc(mesh, solver);
|
|
solver.A = solver.P' * solver.A * solver.P;
|
|
solver.b = solver.P' * solver.b;
|
|
rr = symamd(solver.A);
|
|
solver.A = solver.A(rr, rr);
|
|
solver.b = solver.b(rr);
|
|
solver.x = solver.A \ solver.b;
|
|
q = zeros(length(rr), 1);
|
|
for i = 1:length(rr), q(rr(i)) = i; end
|
|
solver.x = solver.x(q);
|
|
solver.x = solver.P * solver.x;
|
|
triIndex = findTri([2,5], mesh);
|
|
[x,y,z,normE] = get_ele(triIndex, mesh, solver);
|
|
fid = fopen('matlab_norme_stats.txt','w');
|
|
fprintf(fid,'dof=%d\n', solver.dof);
|
|
fprintf(fid,'reduced=%d\n', size(solver.P,2));
|
|
fprintf(fid,'nbrPBC=%d\n', size(mesh.PBCIndex,1));
|
|
fprintf(fid,'max_normE_face25=%.6g\n', max(normE(:)));
|
|
fprintf(fid,'mean_normE_face25=%.6g\n', mean(normE(:)));
|
|
fprintf(fid,'abs_x_max=%.6g\n', max(abs(solver.x)));
|
|
fclose(fid);
|
|
dlmwrite('matlab_x_real.txt', real(solver.x), 'precision', 16);
|
|
dlmwrite('matlab_x_imag.txt', imag(solver.x), 'precision', 16);
|
|
fprintf('OK max|E|=%.6g reduced=%d nbrPBC=%d\n', max(normE(:)), size(solver.P,2), size(mesh.PBCIndex,1)); |