function exportAbMatrixFem4(outDir)
%EXPORTABMATRIXFEM4  Export reduced A/b (after PEC, before symrcm) for Fem4 port case.
%
%   exportAbMatrixFem4()                    % -> OutFile_fem4_ab/
%   exportAbMatrixFem4('OutFile_fem4_ab')
%
% Output format matches C++ Test_OutputMatrix:
%   Ai.txt, Aj.txt (0-based indices), Av.txt (re,im), Bv_real.txt, Bv_imag.txt
%
% Run from this folder:
%   cd('三维matlab代码/2023-2-端口激励问题(四面体网格)')
%   exportAbMatrixFem4

rootDir = fileparts(mfilename('fullpath'));
if nargin < 1 || isempty(outDir)
    outDir = fullfile(rootDir, 'OutFile_fem4_ab');
end
if ~exist(outDir, 'dir')
    mkdir(outDir);
end

exportDir = fullfile(rootDir, '..', 'matlab 3D一阶散射问题');
addpath(rootDir, exportDir);

meshFile = fullfile(rootDir, 'MeshData2x.mat');
S = load(meshFile);
Mesh = S.Mesh;
if isfield(S, 'Physic')
    Physic = S.Physic;
    Material = S.Material;
    [~, ~, Port1Solver, Port2Solver] = configPhysicMaterialFem4();
else
    [Physic, Material, Port1Solver, Port2Solver] = configPhysicMaterialFem4();
end
if ~isfield(Mesh, 'NbrNodes')
    Mesh.NbrNodes = length(Mesh.Nodes);
    Mesh.NbrElements = length(Mesh.Elements);
end

fprintf('Port eigenmode solve (same as MainFem2x)...\n');
Port1Mesh = GetBMesh(Mesh.Nodes, Mesh.Faces, Mesh.FacesIndex, Physic.Port1);
[Port1Solver, Port1Mesh] = BoundaryEigenMode(Port1Mesh, Physic, Material, Port1Solver);
Port1Solver = GetPowerCoef(Port1Solver, Port1Mesh, Material.murF, 1);
Port1Mesh.normal = [0; 0; -1];

Port2Mesh = GetBMesh(Mesh.Nodes, Mesh.Faces, Mesh.FacesIndex, Physic.Port2);
[Port2Solver, Port2Mesh] = BoundaryEigenMode(Port2Mesh, Physic, Material, Port2Solver);
Port2Mesh.normal = [0; 0; 1];

fprintf('Global matrix assembly (wave + ports + PEC reduce)...\n');
epsilonr = Material.epsilonr - 1i * Material.sigma * Material.temp;
[A, Mesh.Edges, Mesh.EdgesOfElements, Dof] = PhysicMatrixAssembly( ...
    Mesh.Nodes, Mesh.Elements, Mesh.Domains, Physic.lam0, epsilonr, Material.mur);
b = complex(zeros(Dof, 1));

el2no = Mesh.Elements';
f1 = el2no([1 1 1 2], :);
f2 = el2no([2 2 3 3], :);
f3 = el2no([3 4 4 4], :);
f_array = [f1(:) f2(:) f3(:)];
[Faces, ~, FacesOfElements] = unique(f_array, 'rows');
Mesh.AllFaces = Faces;
Mesh.FacesOfElements = FacesOfElements;
addDof = 0;

Port1Mesh.FacesConn = GetBFacesConn(Port1Mesh.OldFaces, Faces, FacesOfElements);
[A, b, addDof] = PortBCMatrixAssembly(A, b, Mesh.Nodes, Mesh.Elements, Mesh.EdgesOfElements, ...
    Mesh.Domains, Material.mur, Port1Mesh, Port1Solver, Dof, addDof);

Port2Mesh.FacesConn = GetBFacesConn(Port2Mesh.OldFaces, Faces, FacesOfElements);
[A, b, addDof] = PortBCMatrixAssembly(A, b, Mesh.Nodes, Mesh.Elements, Mesh.EdgesOfElements, ...
    Mesh.Domains, Material.mur, Port2Mesh, Port2Solver, Dof, addDof);
addDof(1) = [];

[PECFaces, ~] = GetBFaces(Mesh.Faces, Mesh.FacesIndex, Physic.PEC);
fl2no = PECFaces';
n1 = fl2no([1 1 2], :);
n2 = fl2no([2 3 3], :);
fl_fd2no_array = [n1(:) n2(:)];
PECEdges = unique(fl_fd2no_array, 'rows');
PECInd = find(ismember(Mesh.Edges, PECEdges, 'rows'));
Free_ind = 1:(Dof + sum(addDof));
Free_ind(PECInd) = [];
A = A(Free_ind, Free_ind);
b = b(Free_ind);

fprintf('Reduced system: size=%d, nnz=%d, |b|=%g\n', ...
    size(A, 1), nnz(A), norm(b));

export_ab_matrix(A, b, outDir);

metaFile = fullfile(outDir, 'meta.txt');
fid = fopen(metaFile, 'w');
fprintf(fid, 'source=MATLAB exportAbMatrixFem4\n');
fprintf(fid, 'mesh=MeshData2x.mat\n');
fprintf(fid, 'stage=after_PEC_before_symrcm\n');
fprintf(fid, 'n=%d\n', size(A, 1));
fprintf(fid, 'nnz=%d\n', nnz(A));
fprintf(fid, 'nPEC=%d\n', numel(PECInd));
fprintf(fid, 'Dof=%d\n', Dof);
fprintf(fid, 'addDof=%d\n', sum(addDof));
fprintf(fid, 'powerCoef=%.15g\n', Port1Solver.powerCoef);
fclose(fid);

fprintf('Done. Compare with C++ OutFile/Ai.txt, Aj.txt, Av.txt, Bv_*.txt\n');
fprintf('  MATLAB: %s\n', outDir);
fprintf('  C++:    port/Release/OutFile/\n');

end