function edgeIds = find_pec_edges(pecDomains, mesh)
%FIND_PEC_EDGES  收集 PEC 边界三角面上的边自由度（1-based edge id）

triIdx = findTri(pecDomains, mesh);
edgeIds = zeros(numel(triIdx) * 3, 1);

for n = 1:numel(triIdx)
    numTet = mesh.ConnOfTri(triIdx(n), 1);
    numFace = mesh.ConnOfTri(triIdx(n), 2);
    if numFace == 1
        e = mesh.EdgeOfTet(numTet, [1, 2, 4]);
    elseif numFace == 2
        e = mesh.EdgeOfTet(numTet, [1, 3, 5]);
    elseif numFace == 3
        e = mesh.EdgeOfTet(numTet, [2, 3, 6]);
    else
        e = mesh.EdgeOfTet(numTet, [4, 5, 6]);
    end
    edgeIds((n - 1) * 3 + (1:3)) = e;
end

edgeIds = unique(edgeIds);

end