function [P,T,ID] = diskretizujObdelnik(a,b,h)

m1 = ceil(a/h);
x1 = linspace(0,a,m1+1);
m2 = ceil(b/h);
x2 = linspace(0,b,m2+1); 

m = 2*m1*m2;
n = (m1+1)*(m2+1);

[X1,X2] = meshgrid(x1,x2);
P = [X1(:),X2(:)]';

T = zeros(3,m);
ID = [];

idxT = 1;
for i=1:m1
    idxP = (i-1)*(m2+1);
    for j=1:m2
        T(:,idxT) = idxP + [j;m2+1+j;m2+1+j+1];
        T(:,idxT+1) = idxP + [j;m2+1+j+1;j+1];
        idxT = idxT + 2;
    end
end

ID = unique([1:(m2+1),n:-1:(n-m2),1:(m2+1):n,(m2+1):(m2+1):n]);