function [A,b] = MKP_elektrostatika(Uzly,Trojuhelniky,Hranice,Oblasti,epsilonr,q)

epsilon0 = 8.854e-12;

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

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

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

A = sparse(n,n); % nevhodne, lepe pouzit sparse(n,n)
b = zeros(n,1);

for k=1:m
    
    Trojuhelnikk = Trojuhelniky(:,k);
    Uzlyk = Uzly(:,Trojuhelnikk);
    
    if Oblasti(k)==1 % dielektrikum
        epsilonk = epsilon0*epsilonr;
    else
        epsilonk = epsilon0;
    end
    
    qk = 0;
    if Oblasti(k)==2
        qk = -q;
    elseif Oblasti(k)==3
        qk = q;
    end
    
    Ak = lokalni_matice_elektrostatika(Uzlyk,epsilonk);
    bk = lokalni_vektor_elektrostatika(Uzlyk,qk);
    
    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,Hranice) = speye(nHranice);
b(Hranice) = 0;
