function [mapM, vmapM, mapP, vmapP, mapB, vmapB] = BuildMaps2D()

% connectivity and boundary tables in the K # of Np elements
% vmapM, vmapP and vmapB - global node numbering in the mesh
% mapM, mapP and mapB - numbering with node duplicity in mind

Globals2D;

nodeIDs = reshape(1:(K * N_p),N_p,K);
vmapM = zeros(N_fp, N_faces, K);
vmapP = zeros(N_fp, N_faces, K);
mapM = (1:(K * N_fp * N_faces))';
mapP = reshape(mapM, N_fp, N_faces, K);

for element1 = 1:K
    for face1 = 1:N_faces
        vmapM(:,face1,element1) = nodeIDs(Fmask(:,face1),element1);
    end
end

for element1 = 1:K
    for face1 = 1:N_faces
        element2 = EtoE(element1,face1);
        face2 = EtoF(element1,face1);
        
        vertex1 = E(face1,element1);
        vertex2 = E(1 + mod(face1,N_faces),element1);
        ref_length = sqrt((P(1,vertex1) - P(1,vertex2)) ^ 2 + (P(2,vertex1) - P(2,vertex2)) ^ 2);
        
        vertexIDM = vmapM(:,face1,element1);
        vertexIDP = vmapM(:,face2,element2);
        x1 = x(vertexIDM); x2 = x(vertexIDP);
        y1 = y(vertexIDM); y2 = y(vertexIDP);
        
        % vytvoreni ctvercovych matic z daneho sloupcoveho vektoru
        X1 = zeros(N_fp); X2 = X1; Y1 = X1; Y2 = X1;
        for i = 1:N_fp
            X1(:,i) = x1;
            X2(:,i) = x2;
            Y1(:,i) = y1;
            Y2(:,i) = y2;
        end
        
        D = (X1 - X2') .^ 2 + (Y1 - Y2') .^ 2;
        [IDM, IDP] = find(sqrt(D) < NODETOL * ref_length);
        vmapP(IDM,face1,element1) = vertexIDP(IDP);
        mapP(IDM,face1,element1) = IDP + (face2 - 1) * N_fp + (element2 - 1) * N_faces * N_fp; % mapP je od 1 do K * N_faces * N_fp
    end
end

vmapP = vmapP(:);
vmapM = vmapM(:);
mapP = mapP(:);

mapB = find(vmapM == vmapP);
vmapB = vmapM(mapB);

end