% PML

sigma_x_PML = zeros(N_p,K);
sigma_y_PML = zeros(N_p,K);

D_PML = 10;
sigma_m = 2;
beta = 2;

indexes_R = find(x > x_u + NODETOL);
indexes_L = find(x < x_l - NODETOL);
% toto jeste doresit, jaky tam ma byt nasobek pro sigma_x
sigma_x_PML(indexes_R) = sigma_m * (1 - u_0 ^ 2) * (abs(x(indexes_R) - x_u) / D_PML) .^ beta;
sigma_x_PML(indexes_L) = sigma_m * (1 - u_0 ^ 2) * (abs(x(indexes_L) - x_l) / D_PML) .^ beta;

indexes_T = find(y > y_u + NODETOL);
indexes_B = find(y < y_l - NODETOL);
sigma_y_PML(indexes_T) = sigma_m * (abs(y(indexes_T) - y_u) / D_PML) .^ beta;
sigma_y_PML(indexes_B) = sigma_m * (abs(y(indexes_B) - y_l) / D_PML) .^ beta;

sigma_x_PML_plus_sigma_y_PML = sigma_x_PML + sigma_y_PML;
sigma_x_PML_times_sigma_y_PML = sigma_x_PML .* sigma_y_PML;


% PML access information
[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);

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);
vmapPML_M = vmapM(mapPML_M);
vmapPML_P = vmapM(mapPML_P);