function [Ex, Ey, Ez, normE] = get_ele_vertices(mesh, solver)
%GET_ELE_VERTICES  Recover E at all mesh vertices (1st-order Nedelec).

Ex = zeros(mesh.NbrVertex, 1);
Ey = zeros(mesh.NbrVertex, 1);
Ez = zeros(mesh.NbrVertex, 1);
nEx = zeros(mesh.NbrVertex, 1);
nEy = zeros(mesh.NbrVertex, 1);
nEz = zeros(mesh.NbrVertex, 1);

u = [1, 0, 0, 0];
v = [0, 1, 0, 0];
w = [0, 0, 1, 0];

for n = 1:mesh.NbrTet
    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);

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

    for j = 1:4
        vi = mesh.Tet(n, j);
        for i = 1:6
            coef = solver.x(mesh.EdgeOfTet(n, i));
            Ex(vi) = Ex(vi) + E(1, i, j) * coef;
            Ey(vi) = Ey(vi) + E(2, i, j) * coef;
            Ez(vi) = Ez(vi) + E(3, i, j) * coef;
        end
        nEx(vi) = nEx(vi) + 1;
        nEy(vi) = nEy(vi) + 1;
        nEz(vi) = nEz(vi) + 1;
    end
end

Ex = Ex ./ nEx;
Ey = Ey ./ nEy;
Ez = Ez ./ nEz;
normE = sqrt(abs(Ex .* conj(Ex) + Ey .* conj(Ey) + Ez .* conj(Ez)));

end