function konvergenceMKP2d()

% -laplace(u(x)) = -2*(x1*(x1-1)+x2*(x2-1)) x z (0,1)^2
% u(x) = 0
% analyt. reseni: u(x) = x1*(x1-1)*x2*(x2-1),
% grad(u) = [(2*x1-1)*x2*(x2-1);(2*x2-1)*x1*(x1-1)]

load GaussQuadratureForTriangles
w = weights{6};
ksi = points{6};
e = w; e(:) = 1;

hcka = 0.5*2.^[0:-1:-3];

for i=1:length(hcka)
    h = hcka(i);
    [P,T,ID,IN] = diskretizujObdelnik(1,1,h);
    m = size(T,2);
    mN = 0; %size(IN,2);
    nvlnka = size(P,2);
    ID = find(abs(P(1,:))<1e-8 | abs(P(1,:)-1)<1e-8 | abs(P(2,:))<1e-8 | abs(P(2,:)-1)<1e-8);
    IN = [];
    nD = length(ID);

    % data ulohy
    k = ones(1,m);
    Pmid = (P(:,T(1,:))+P(:,T(2,:))+P(:,T(3,:)))/3;
    f = -2*(Pmid(1,:).*(Pmid(1,:)-1)+Pmid(2,:).*(Pmid(2,:)-1));
    uDirich = zeros(nD,1);
    idxNonDirich = setdiff(1:nvlnka,ID);
    tauNeum = ones(1,mN);

    % MKP
    [A,b] = mkp(P,T,k,f,IN,tauNeum,ID,uDirich);
    bvlnka = zeros(nvlnka,1);
    for j=1:m
        Tj = T(:,j);
        Rj = [P(:,Tj(2))-P(:,Tj(1)),P(:,Tj(3))-P(:,Tj(1))];
        xj = P(:,Tj(1))*e+Rj*ksi;
        fj = -2*(xj(1,:).*(xj(1,:)-1)+xj(2,:).*(xj(2,:)-1));
        bj = abs(det(Rj)) * [w*(fj.*(1-ksi(1,:)-ksi(2,:)))';w*(fj.*ksi(1,:))';w*(fj.*ksi(2,:))'];
        bvlnka(Tj) = bvlnka(Tj) + bj;
    end
        
    uvlnka = zeros(nvlnka,1);
    uH = A\b;
    uvlnka(idxNonDirich) = uH;
    uvlnka(ID) = uDirich;
    
%     trisurf(T',P(1,:),P(2,:),uvlnka')
%     pause

    gradu = zeros(2,m);
    err(i) = 0;
    for j=1:m
        Tj = T(:,j);
        Rj = [P(:,Tj(2))-P(:,Tj(1)),P(:,Tj(3))-P(:,Tj(1))];
        uj = uvlnka(Tj);
        graduj = (Rj')\([[-1;-1] [1;0] [0;1]] * uj);
        xj = P(:,Tj(1))*e+Rj*ksi;
        errj = graduj*e - [(2*xj(1,:)-1).*xj(2,:).*(xj(2,:)-1);(2*xj(2,:)-1).*xj(1,:).*(xj(1,:)-1)];
        errj = w * (errj(1,:).^2+errj(2,:).^2)' * abs(det(Rj));
        err(i) = err(i) + errj;
    end
    err(i) = sqrt(err(i));
end

plot(log(hcka),log(err),'b',log(hcka),log(hcka),'r--')
xlabel('log(h)')
ylabel('log(H1error)')