function [noveUzly,noveHrany,novePrvky,noveOkrHrany] = refine2d(uzly,hrany,prvky,okrHrany)

nuzlu = size(uzly,2);
nhran = size(hrany,2);
nprvku = size(prvky,2);
nokrHran = length(okrHrany);

noveUzly = zeros(2,nuzlu+nhran);
noveHrany = zeros(2,2*nhran+3*nprvku);
novePrvky = zeros(3,4*nprvku);
noveOkrHrany = zeros(1,2*nokrHran);


noveUzly(:,1:nuzlu) = uzly;

for i=1:nhran
    hranai = hrany(:,i);
    
    novyUzel = 0.5*sum(uzly(:,hranai),2);
    noveUzly(:,nuzlu+i) = novyUzel;
    
    noveHrany(:,2*i-1:2*i) = [ [hranai(1);nuzlu+i] , [nuzlu+i;hranai(2)] ];
end

idxNoveHrany = 2*nhran+1;
for i=1:nprvku
    prveki = prvky(:,i);
    j = prveki(1);
    k = prveki(2);
    l = prveki(3);
    
    % Takto vyhledevat hrany je nevhodne, vypocetni slozitost je
    % O(nprvku*nhran), tj. O(n^2), coz by zrusilo slozitost O(n)
    % multigridu. Spravne je pamatovat si take pole prvek->hrany.
    for h=1:nhran
        if (hrany(1,h)==j && hrany(2,h)==k) || (hrany(1,h)==k && hrany(2,h)==j)
            break;
        end
    end
    p = h;
    for h=1:nhran
        if (hrany(1,h)==k && hrany(2,h)==l) || (hrany(1,h)==l && hrany(2,h)==k)
            break;
        end
    end
    q = h;
    for h=1:nhran
        if (hrany(1,h)==l && hrany(2,h)==j) || (hrany(1,h)==j && hrany(2,h)==l)
            break;
        end
    end
    r = h;
    
    idx1 = nuzlu+p;
    idx2 = nuzlu+q;
    idx3 = nuzlu+r;
    
    novePrvky(:,4*(i-1)+[1:4]) = [ [j;idx1;idx3] , [idx3;idx2;l] , [idx1;k;idx2] , [idx1;idx2;idx3] ];
    noveHrany(:,idxNoveHrany:idxNoveHrany+2) = [ [idx1;idx2] , [idx2;idx3] , [idx3;idx1] ];
    
    idxNoveHrany = idxNoveHrany+3;
end

for i=1:nokrHran
    noveOkrHrany(2*i-1:2*i) = [2*okrHrany(i)-1,2*okrHrany(i)];
end