function [n_x,n_y,surfJ] = Normals2D()

Globals2D;

x_r = D_r * x; x_s = D_s * x;
y_r = D_r * y; y_s = D_s * y;

facex_r = x_r(Fmask(:),:); facex_s = x_s(Fmask(:),:);
facey_r = y_r(Fmask(:),:); facey_s = y_s(Fmask(:),:);

n_x = zeros(N_faces * N_fp,K); n_y = n_x;
face1IDs = (1:N_fp)'; face2IDs = ((N_fp + 1):(2 * N_fp))'; face3IDs = ((2 * N_fp + 1):(3 * N_fp))';

% je zrejme, ze u face 1 a face 2 staci prohodit x a y, ale navolit znamenka, aby se jednalo o vnejsi normalu?

% face 1 - 1. and 2. triangle vertex
n_x(face1IDs,:) = facey_r(face1IDs,:);
n_y(face1IDs,:) = -facex_r(face1IDs,:);

% face 2 - 2. and 3. triangle vertex
n_x(face2IDs,:) = facey_s(face2IDs,:) - facey_r(face2IDs,:);
n_y(face2IDs,:) = -facex_s(face2IDs,:) + facex_r(face2IDs,:);

% face 3 - 3. and 1. triangle vertex
n_x(face3IDs,:) = -facey_s(face3IDs,:);
n_y(face3IDs,:) = facex_s(face3IDs,:);

surfJ = sqrt(n_x .* n_x + n_y .* n_y);
n_x = n_x ./ surfJ; n_y = n_y ./ surfJ;

end