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

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 = zeros(12,12);  % memory allocation

M = 1/12*abs(V)*diag([1 1 1 1 1 1 1 1 1 1 1 1]);

M__  = 1/12*V*[2 0 1 0 1 0;
               0 2 0 1 0 1;
               1 0 2 0 1 0;
               0 1 0 2 0 1;
               1 0 1 0 2 0;
               0 1 0 1 0 2];


% displacement vector: [x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4]