function [Uzly,Trojuhelniky,HranicniUzlyMinus,HranicniUzlyPlus,Oblasti] = diskretizace3_elektrostatika_MKP(h)

% [Uzly,Trojuhelniky,HranicniUzly,Oblasti] = diskretizace3_elektrostatika_MKP(h)
%       uvazujeme oblasti Omega := (-5,5)x(-5,5), Omegar = (-1,1)x(-1,1),
%       h znaci krok diskretizace, tj. delku hrany trojuhelniku.

sirkaOmega = 10; %28;
vyskaOmega = 10; %20;

m1 = sirkaOmega/h;
m2 = vyskaOmega/h;
m = 2*m1*m2;        % pocet trojuhelniku

n1 = m1+1;
n2 = m2+1;
n = n1*n2;          % pocet uzlu

sirkaOmegar = 2;
vyskaOmegar = 2;

Uzly = zeros(2,n);
Trojuhelniky = zeros(3,m);
HranicniUzlyMinus = zeros(1,n2);
HranicniUzlyPlus = zeros(1,n2);
Oblasti = zeros(1,m);

x1 = -sirkaOmega/2:h:sirkaOmega/2;
x2 = -vyskaOmega/2:h:vyskaOmega/2;

idxUzly = 1;
idxTrojuhelniky = 1;
idxHraniceMinus = 1;
idxHranicePlus = 1;

for i1=1:n1
    for i2=1:n2
        
        Uzly(:,idxUzly) = [x1(i1);x2(i2)];
        
        if i1<n1 && i2<n2
            Trojuhelniky(:,idxTrojuhelniky) = [idxUzly;idxUzly+n2;idxUzly+n2+1];
            Trojuhelniky(:,idxTrojuhelniky+1) = [idxUzly;idxUzly+n2+1;idxUzly+1];
            
            if x1(i1)>=-sirkaOmegar/2 && x1(i1)<sirkaOmegar/2 && x2(i2)>=-vyskaOmegar/2 && x2(i2)<vyskaOmegar/2
                Oblasti(idxTrojuhelniky) = 1;
                Oblasti(idxTrojuhelniky+1) = 1;
            else
                Oblasti(idxTrojuhelniky) = 0;
                Oblasti(idxTrojuhelniky+1) = 0;
            end
            
            idxTrojuhelniky = idxTrojuhelniky+2;
        end
        
        if i1==1
            HranicniUzlyMinus(idxHraniceMinus) = idxUzly;
            idxHraniceMinus = idxHraniceMinus+1;
        elseif i1==n1
            HranicniUzlyPlus(idxHranicePlus) = idxUzly;
            idxHranicePlus = idxHranicePlus+1;
        end
        
        idxUzly = idxUzly+1;
        
    end
end
        
        