function [A,f] = assemble(ey,ec,x,e,ID,IN)

% [A,f] = assemble(ey,ec,x,e,ID,IN)
%          It assembles a system (permeability) matrix
%          and a right-hand side vector both arising
%          from a finite-element approximation
%          of the following 2d magnetostatic problem:
%           Find uH from Vh : a(uH,v) = f(v) for all v from Vh,
%          where
%           a(v,u) = \int_{Omega} (1/mu(x))*grad(uH(x))*grad(v(x)) dx,
%           f(v) = \int_{Omega} J(x)*v(x) dx,
%           Vh is given by Lagrange shape functions
%           with the homogeneous Dirichlet boundary condition
%           on Gamma_D, specified by the indices ID,
%           and with the homogeneous Neumann b.c. on Gamma_N,
%           specified by the indices IN.

%sprintf('%s','assemble')

% number of nodes (= number of degrees of freedom)
n = max(size(x));

% number of elements
ne = max(size(e));

A = sparse(n,n);
f = zeros(n,1);

% permeability of the air
mu0 = 4*pi*1e-7;
% permeability of the ferromagnetic
mur = mu0*5100;

% current density
J = 10e6;

for i=1:ne
%    if mod(i,100)==0
%        i
%    end
    ei = e(i,:);
    
    if in(ey,i)
        mui = mur;
    else
        mui = mu0;
    end
    
    Aei = aloc_mu(mui,x,ei);
    
    if in(ec,i)
        Ji = J;
    else
        Ji = 0;
    end

    bei = bloc_J(Ji,x,ei);
    
    nei = max(size(Aei));
    
    for j=1:nei
        for k=1:nei          
            if j==k | (~in(ID,ei(j)) & ~in(ID,ei(k)))
                A(ei(j),ei(k)) = A(ei(j),ei(k)) + Aei(j,k);
            end
        end
        if ~in(ID,ei(j))
            f(ei(j)) = f(ei(j)) + bei(j);
        end
    end
end
