function [A,b] = MKP_magnetostatika(Uzly,Trojuhelniky,Hranice,Oblasti,mur,J)

mu0 = 4*pi*1e-7;

% pocet uzlu
n = size(Uzly,2);

% pocet trojuhelniku
m = size(Trojuhelniky,2);

% pocet uzlu na hranici
nHranice = length(Hranice);

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

for k=1:m
    
    Trojuhelnikk = Trojuhelniky(:,k);
    Uzlyk = Uzly(:,Trojuhelnikk);
    
    if Oblasti(k)==1 % feromagnetikum
        nuk = 1/mur;
    else
        nuk = 1;
    end
    
    Jk = 0;
    if Oblasti(k)==2     % leva cast civky
        Jk = -mu0*J;
    elseif Oblasti(k)==3 % prava cast civky
        Jk = mu0*J;
    end
    
    Ak = lokalni_matice_magnetostatika(Uzlyk,nuk);
    bk = lokalni_vektor_magnetostatika(Uzlyk,Jk);
    
    A(Trojuhelnikk,Trojuhelnikk) = A(Trojuhelnikk,Trojuhelnikk) + Ak;
    b(Trojuhelnikk) = b(Trojuhelnikk) + bk;
    
end

% hranicni podminky: ui = 0 pro xi z Gamma
A(Hranice,:) = 0;
A(:,Hranice) = 0;
A(Hranice,Hranice) = speye(nHranice);
b(Hranice) = 0;
