[rows1,cols1] = find(x > x_u + NODETOL);
[rows2,cols2] = find(x < x_l - NODETOL);
[rows3,cols3] = find(y > y_u + NODETOL);
[rows4,cols4] = find(y < y_l - NODETOL);
elements_PML = unique([cols1; cols2; cols3; cols4]);

elements_Euler = setdiff((1:K)',elements_PML);

% postup jako v BuildBCMaps2D od Hesthavena
% verze cislo 1
%elements_PML_with_neighbours = unique(EtoE(elements_PML,:));
%Faces_PML = zeros(3,K);
%Faces_PML(:,elements_PML_with_neighbours) = ones(3,length(elements_PML_with_neighbours));
%Faces_PML = Faces_PML(:);
%Nodes_PML = ones(N_fp,1) * Faces_PML';
%Nodes_PML = Nodes_PML(:);
%mapPML = find(Nodes_PML == 1);
%vmapPML = vmapM(mapPML);

% verze cislo 2
Faces_PML = zeros(3,K);
Faces_PML(:,elements_PML) = ones(3,length(elements_PML));
Faces_PML = Faces_PML(:);
Nodes_PML = ones(N_fp,1) * Faces_PML';
Nodes_PML = Nodes_PML(:);
mapPML = find(Nodes_PML == 1); % faces nalezici elementum v PML
mapPML_M = unique([mapPML;mapP(mapPML)]); % body na faces na vsech elementech lezicich v PML nebo sousedicich s elementem v PML
mapPML_P = mapP(mapPML_M); % prilehajici body na body z mapPML_P
%vmapPML = vmapM(mapPML);

figure;
trimesh(E',P(1,:),P(2,:),'Color','black'); hold on;
%trimesh(E(:,elements_PML_with_neighbours)',P(1,:),P(2,:),'Color','red','LineWidth',1.5); hold on;
trimesh(E(:,elements_PML)',P(1,:),P(2,:),'Color','red','LineWidth',1.5); hold on;
plot([x_l x_u x_u x_l x_l],[y_l y_l y_u y_u y_l],'Color','blue','LineWidth',2); hold on;
X = [(x_l - D_PML) (x_u + D_PML) (x_u + D_PML) (x_l - D_PML) (x_l - D_PML)];
Y = [(y_l - D_PML) (y_l - D_PML) (y_u + D_PML) (y_u + D_PML) (y_l - D_PML)];
plot(X,Y,'Color','blue','LineWidth',2);
axis equal;
xlim([(x_l - D_PML) (x_u + D_PML)]); ylim([(y_l - D_PML) (y_u + D_PML)]);

figure;
trimesh(E',P(1,:),P(2,:),'Color','black'); hold on;
plot([x_l x_u x_u x_l x_l],[y_l y_l y_u y_u y_l],'Color','red','LineWidth',2); hold on;
plot(X,Y,'Color','red','LineWidth',2); hold on;
%plot(x(vmapPML),y(vmapPML),'*');
plot(Fx(mapPML_M),Fy(mapPML_M),'*');
axis equal;
xlim([(x_l - D_PML) (x_u + D_PML)]); ylim([(y_l - D_PML) (y_u + D_PML)]);