D3_TETRAH (x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, Em, miu)

function [K,M,B,E] = D3_TETRAH (x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, Em, miu)

% returns a 3D tetrahedron (CST) element K & M matrix
% in global coordinates
% INPUT:
%       x1, y1, z1 - coordinates of joint 1 of the triangle element
%       x2, y2, z2 - coordinates of joint 2 of the triangle element
%       x3, y3, z3 - coordinates of joint 3 of the triangle element
%       x4, y4, z4 - coordinates of joint 3 of the triangle element
%       Em     - elastic modulus for isotropic material
%       miu    - Poisson's ratio
% ---------------------------------------------------------------------------

J = [1 1 1 1; x1 x2 x3 x4; y1 y2 y3 y4; z1 z2 z3 z4];
detJ = det(J);
V = 1/6 * detJ;
Q = inv(J)*detJ;

a1 = Q(1,2); b1 = Q(1,3); c1 = Q(1,4);
a2 = Q(2,2); b2 = Q(2,3); c2 = Q(2,4);
a3 = Q(3,2); b3 = Q(3,3); c3 = Q(3,4);
a4 = Q(4,2); b4 = Q(4,3); c4 = Q(4,4);

s=(1-miu);
E  = Em*(1-miu)/((1+miu)*(1-2*miu)).*...
         [   1   miu/s  miu/s       0                0                0;
           miu/s   1    miu/s       0                0                0;
           miu/s miu/s    1         0                0                0;
             0     0      0   (1-2*miu)/(2*s)        0                0;
             0     0      0         0         (1-2*miu)/(2*s)         0;
             0     0      0         0                0        (1-2*miu)/(2*s)];

B = 1/(6*V)*[ a1    0   0  a2   0   0   a3   0   0   a4   0   0;
               0   b1   0  0   b2   0   0   b3   0   0   b4   0;
               0    0  c1  0    0  c2   0    0  c3   0    0  c4;
               0   c1  b1  0   c2  b2   0   c3  b3   0   c4  b4;
              c1    0  a1 c2    0  a2  c3    0  a3  c4    0  a4;
              b1   a1   0 b2   a2   0  b3   a3   0  b4   a4   0];

K = V.*B'*E*B ;


M  = V/20*  [2 0 0 1 0 0 1 0 0 1 0 0;
               0 2 0 0 1 0 0 1 0 0 1 0;
               0 0 2 0 0 1 0 0 1 0 0 1;
               1 0 0 2 0 0 1 0 0 1 0 0;
               0 1 0 0 2 0 0 1 0 0 1 0;
               0 0 1 0 0 2 0 0 1 0 0 1;
               1 0 0 1 0 0 2 0 0 1 0 0;
               0 1 0 0 1 0 0 2 0 0 1 0;
               0 0 1 0 0 1 0 0 2 0 0 1;
               1 0 0 1 0 0 1 0 0 2 0 0;
               0 1 0 0 1 0 0 1 0 0 2 0;
               0 0 1 0 0 1 0 0 1 0 0 2];

% displacement vector: [x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4]
% node 4 is in the back
%     2 o
%      .  .
%     .    .
%    .   o4 .
%   o........o
%  1          3
% ------------------end