%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-04-14 17:43:46
 % @LastEditTime: 2022-04-25 10:42:27
%--------------------------
function Nh = init_Nh_triangle(nodes, faces)
    x = nodes(:, 1); % [n, 1]
    x = x(faces);    % [m, 3]
    y = nodes(:, 2); % [n, 1]
    y = y(faces);    % [m, 3]
    [a, b, c, area] = paras(x, y); % [n,3]

    Nh = zeros(size(x,1), 3, 3); % [m,3,3], [ele_num, nodes, 3 for 6]
    for i = 1:3
        % Nhj = zeros(size(x,1),1,3);
        for j = 1:3
            Nh(:,i,j) = (a(:,j) + b(:,j) .* x(:,i) + c(:,j) .* y(:,i)) / 2 ./ area;
        end
        % Nh(:,i,:) = kron(Nhj, eye(2)); % [2,6]
    end
    % Nh1 = kron(Nh, eye(2)); % [m,6]
end

function [a,b,c, area] = paras(x, y)
    % a = [a1, a2, a3], n*3 , ...
    a1 = x(:,2) .* y(:,3) - x(:,3) .* y(:,2);
    a2 = x(:,3) .* y(:,1) - x(:,1) .* y(:,3);
    a3 = x(:,1) .* y(:,2) - x(:,2) .* y(:,1);
    a = [a1, a2, a3];

    b1 = y(:,2) - y(:,3);
    b2 = y(:,3) - y(:,1);
    b3 = y(:,1) - y(:,2);
    b = [b1, b2, b3];

    c1 = -x(:,2) + x(:,3);
    c2 = -x(:,3) + x(:,1);
    c3 = -x(:,1) + x(:,2);
    c = [c1, c2, c3];

    area = 0.5 * ((x(:,1).* y(:,2) - x(:,2).*y(:,1)) + ...
                  (x(:,2).* y(:,3) - x(:,3).*y(:,2)) + ...
                  (x(:,3).* y(:,1) - x(:,1).*y(:,3)));
    area = abs(area); % important!
end