function solver = assembly_bele(phy, mesh, solver)
%ASSEMBLY_BELE  BELE 背景场源项：方法 A 直接体积分
%
%  b -= ∫ J_b·N dV，J_b = curl curl E_b − k₀²ε E_b  （弱形式 Tt − St）
%  不拆 Tt−St，不在 SBC 外表面做 IBP / Be(E_b)

k0 = 2 * pi / phy.lda0;

[u, v, w, weight, nbrGP] = getGaussPoints(1);

for n = 1:mesh.NbrTet
    domain = mesh.DomainOfTet(n);
    if ~ismember(domain, phy.beleDomains)
        continue;
    end

    x = mesh.Vertex(mesh.Tet(n, :), 1);
    y = mesh.Vertex(mesh.Tet(n, :), 2);
    z = mesh.Vertex(mesh.Tet(n, :), 3);

    Jac = zeros(3, 3);
    Jac(1, 1) = x(1) - x(4); Jac(1, 2) = y(1) - y(4); Jac(1, 3) = z(1) - z(4);
    Jac(2, 1) = x(2) - x(4); Jac(2, 2) = y(2) - y(4); Jac(2, 3) = z(2) - z(4);
    Jac(3, 1) = x(3) - x(4); Jac(3, 2) = y(3) - y(4); Jac(3, 3) = z(3) - z(4);
    DetJac = abs(det(Jac));

    E = zeros(3, 6, nbrGP);
    for gp = 1:nbrGP
        for j = 1:6
            E(:, j, gp) = getBF(1, j, u(gp), v(gp), w(gp));
            E(:, j, gp) = Jac \ E(:, j, gp);
        end
    end

    eps = phy.eps(domain);

    for i = 1:6
        contrib = 0;
        for gp = 1:nbrGP
            px = x(4) + Jac(1, 1) * u(gp) + Jac(2, 1) * v(gp) + Jac(3, 1) * w(gp);
            py = y(4) + Jac(1, 2) * u(gp) + Jac(2, 2) * v(gp) + Jac(3, 2) * w(gp);
            pz = z(4) + Jac(1, 3) * u(gp) + Jac(2, 3) * v(gp) + Jac(3, 3) * w(gp);
            bE = phy.bele.Eb(px, py, pz);
            bCurlCurlE = phy.bele.curlcurlEb(px, py, pz);
            Jb = bCurlCurlE - k0^2 * eps * bE;
            contrib = contrib + weight(gp) * DetJac * dot(E(:, i, gp), Jb);
        end
        edgeId = mesh.EdgeOfTet(n, i);
        % Weak-form BELE: b += Tt - St (same as 2D/C++ after IPP with A = St - Tt)
        solver.b(edgeId) = solver.b(edgeId) - contrib;
    end
end

end