%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-26 16:15:45
 % @LastEditTime: 2022-01-26 19:14:58
%--------------------------
function [sK, v] = intKE_tet(node, tet, E, NU)
    
    x = node(:,1); % n,1
    y = node(:,2);
    z = node(:,3);
    x = x(tet); % m,4
    y = y(tet);
    z = z(tet);
    
    [a,b,c,v] = abc(x,y,z);
 
    D = (E/((1+NU)*(1-2*NU)))* ...
            [   1-NU NU NU 0 0 0 ; 
                NU 1-NU NU 0 0 0 ; 
                NU NU 1-NU 0 0 0 ;
                0 0 0 (1-2*NU)/2 0 0 ; 
                0 0 0 0 (1-2*NU)/2 0 ; 
                0 0 0 0 0 (1-2*NU)/2];

    sK = intKE(a,b,c,v, D);
    v = v/6;
end

function sK = intKE(a,b,c,v, D)
 
    ntet = size(a, 1);
    sK = zeros(12*12, ntet);

    for i = 1 : ntet
        Be = 1 / v(i) * ...
                [a(i,1) 0 0 a(i,2) 0 0 a(i,3) 0 0 a(i,4) 0 0;
                 0 b(i,1) 0 0 b(i,2) 0 0 b(i,3) 0 0 b(i,4) 0;
                 0 0 c(i,1) 0 0 c(i,2) 0 0 c(i,3) 0 0 c(i,4);
                 b(i,1) a(i,1) 0 b(i,2) a(i,2) 0 b(i,3) a(i,3) 0 b(i,4) a(i,4) 0;
                 0 c(i,1) b(i,1) 0 c(i,2) b(i,2) 0 c(i,3) b(i,3) 0 c(i,4) b(i,4);
                 c(i,1) 0 a(i,1) c(i,2) 0 a(i,2) c(i,3) 0 a(i,3) c(i,4) 0 a(i,4)];
        ke = Be' * D * Be * v(i)/6;
        sK(:, i) = ke(:);
    end

end

function [a,b,c,v] = abc(x,y,z)
    
    x1 = x(:,1); x2 = x(:,2); x3 = x(:,3); x4 = x(:,4); 
    y1 = y(:,1); y2 = y(:,2); y3 = y(:,3); y4 = y(:,4); 
    z1 = z(:,1); z2 = z(:,2); z3 = z(:,3); z4 = z(:,4); 

    a(:,1) =   y2 .* (z4 - z3) - y3 .* (z4 - z2) + y4 .* (z3 - z2);
    a(:,2) = - y1 .* (z4 - z3) + y3 .* (z4 - z1) - y4 .* (z3 - z1);
    a(:,3) =   y1 .* (z4 - z2) - y2 .* (z4 - z1) + y4 .* (z2 - z1);
    a(:,4) = - y1 .* (z3 - z2) + y2 .* (z3 - z1) - y3 .* (z2 - z1);
    b(:,1) = - x2 .* (z4 - z3) + x3 .* (z4 - z2) - x4 .* (z3 - z2);
    b(:,2) =   x1 .* (z4 - z3) - x3 .* (z4 - z1) + x4 .* (z3 - z1);
    b(:,3) = - x1 .* (z4 - z2) + x2 .* (z4 - z1) - x4 .* (z2 - z1);
    b(:,4) =   x1 .* (z3 - z2) - x2 .* (z3 - z1) + x3 .* (z2 - z1);
    c(:,1) =   x2 .* (y4 - y3) - x3 .* (y4 - y2) + x4 .* (y3 - y2);
    c(:,2) = - x1 .* (y4 - y3) + x3 .* (y4 - y1) - x4 .* (y3 - y1);
    c(:,3) =   x1 .* (y4 - y2) - x2 .* (y4 - y1) + x4 .* (y2 - y1);
    c(:,4) = - x1 .* (y3 - y2) + x2 .* (y3 - y1) - x3 .* (y2 - y1);

    % v is 6V, 6 * vol
    v = (x2-x1) .* ((y3-y1) .* (z4-z1) - (y4-y1) .* (z3-z1)) + ...
        (y2-y1) .* ((x4-x1) .* (z3-z1) - (x3-x1) .* (z4-z1)) + ...
        (z2-z1) .* ((x3-x1) .* (y4-y1) - (x4-x1) .* (y3-y1));
end