%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-11-01 20:07:04
 % @LastEditTime: 2021-12-07 20:10:03
%--------------------------
function [sK, area] = intKE_triangle(nodes, faces, E, nu)
    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);

    if any(area < 1e-3)
        error("area tri-mesh is 0!");
    end

    E0 = E / 4 / (1-nu^2) ./ area;
    
    Ke = zeros(6, 6, size(faces,1));

    for r = 1:3
        for s = 1:3
            [k1, k2, k3, k4] = krs(r, s, nu, b, c); % [m,1]

            Ke((r-1)*2+1, (s-1)*2+1, :) = k1;
            Ke((r-1)*2+1, s*2, :)       = k3;
            Ke(r*2, (s-1)*2+1, :)       = k2;
            Ke(r*2, s*2, :)             = k4;
        end
    end
    sK = reshape(Ke, 36, []) .* E0';
end

function [k1, k2, k3, k4] = krs(r, s, nu, b, c)
    k1 = b(:,r) .* b(:,s) + (1-nu)/2 * c(:,r) .* c(:,s);
    k2 = nu * c(:,r) .* b(:,s) + (1-nu)/2 * b(:,r) .* c(:,s);
    k3 = nu * b(:,r) .* c(:,s) + (1-nu)/2 * c(:,r) .* b(:,s);
    k4 = c(:,r) .* c(:,s) + (1-nu)/2 * b(:,r) .* b(:,s);
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