# FEM on tet10 element: negetive determinant at the Gauss point

I am trying to implement a fem code on tet10 elements. I follow the lecture notes for tet10 implementation given in

I have validated my stiffness matrix with the example provided. The numerical values of the Jacobian at the Gauss points add up to the volume of the element when integrated against the unit function.

However when I take a different set of nodes with exactly the same ordering of the nodes, I get a negative determinant. I am now not sure of how to decide my code is correct or not.

I am providing the matlab code snippet for the same both with the original data given in the paper as well as the new data.

clc;
if 0
node = [19.39935000986642   0.1580775331976262  9.374695831233881
19.85814434000072   0.6598506371289292  0.0841144105427571
19.82724989196882   0.4316427461960032  9.89092611758984
0.1049331859852762  0.5417518219571752  9.70657837113129
19.62874717493357   0.4089640851632777  4.729405120888319
19.84269711598477   0.5457466916624663  4.987520264066299
19.61329995091762   0.2948601396968147  9.632810974411861
9.752141597925849   0.3499146775774007  9.540637101182586
9.981538762992999   0.6008012295430523  4.895346390837023
9.966091538977048   0.4866972840765892  9.798752244360564
];% I get wrong results with same ordering uncomment to run the same
end;
node = [2 3 4
6 3 2
2 5 1
4 3 6
4 3 3
4 4 1.5
2 4 2.5
3 3 5
5 3 4
3 4 3.5]; % I get right results as in the reference
Y=480; %Youngs Modulus
nu=1/3; % Poisson ratio
alfa=(5.0+3.0*sqrt(5.0))/20; %Gauss Points
beta =(5.0-sqrt(5.0))/20;    %Gauss Points
weight=0.25; %weights for integration points
GaussPoints=[alfa beta beta beta;beta alfa beta beta;beta beta alfa  beta;beta beta beta alfa];
%Elasticity Matrix
E=[1-nu nu nu 0 0 0;nu 1-nu nu 0 0 0;nu nu 1-nu 0 0 0;0 0 0 0.5-nu 0 0;0 0 0 0 0.5-nu 0;0 0 0 0 0 0.5-nu];
E=Y/((1+nu)*(1-2*nu))*E;
%Stiffness Matrix
K=zeros(30,30);

for j=1:4
xi=GaussPoints(:,j);
jx1=4.0*(node(1,1)*(xi(1)-0.25)+node(5,1)*xi(2)+node(7,1)*xi(3)+node(8,1)*xi(4));
jy1=4.0*(node(1,2)*(xi(1)-0.25)+node(5,2)*xi(2)+node(7,2)*xi(3)+node(8,2)*xi(4));
jz1=4.0*(node(1,3)*(xi(1)-0.25)+node(5,3)*xi(2)+node(7,3)*xi(3)+node(8,3)*xi(4));

jx2=4.0*(node(5,1)*xi(1)+node(2,1)*(xi(2)-0.25)+node(6,1)*xi(3)+node(9,1)*xi(4));
jy2=4.0*(node(5,2)*xi(1)+node(2,2)*(xi(2)-0.25)+node(6,2)*xi(3)+node(9,2)*xi(4));
jz2=4.0*(node(5,3)*xi(1)+node(2,3)*(xi(2)-0.25)+node(6,3)*xi(3)+node(9,3)*xi(4));

jx3=4.0*(node(7,1)*xi(1)+node(6,1)*xi(2)+node(3,1)*(xi(3)-0.25)+node(10,1)*xi(4));
jy3=4.0*(node(7,2)*xi(1)+node(6,2)*xi(2)+node(3,2)*(xi(3)-0.25)+node(10,2)*xi(4));
jz3=4.0*(node(7,3)*xi(1)+node(6,3)*xi(2)+node(3,3)*(xi(3)-0.25)+node(10,3)*xi(4));

jx4=4.0*(node(8,1)*xi(1)+node(9,1)*xi(2)+node(10,1)*xi(3)+node(4,1)*(xi(4)-0.25));
jy4=4.0*(node(8,2)*xi(1)+node(9,2)*xi(2)+node(10,2)*xi(3)+node(4,2)*(xi(4)-0.25));
jz4=4.0*(node(8,3)*xi(1)+node(9,3)*xi(2)+node(10,3)*xi(3)+node(4,3)*(xi(4)-0.25));

J=[1 1 1 1;jx1 jx2 jx3 jx4;jy1 jy2 jy3 jy4;jz1 jz2 jz3 jz4];

Jdet=det(J);
Jinv=inv(J);
Iaug=[0 0 0;1 0 0;0 1 0;0 0 1];
P=Jinv*Iaug;

a1=P(1,1);
a2=P(2,1);
a3=P(3,1);
a4=P(4,1);

b1=P(1,2);
b2=P(2,2);
b3=P(3,2);
b4=P(4,2);

c1=P(1,3);
c2=P(2,3);
c3=P(3,3);
c4=P(4,3);

Nfx=[(4.0*xi(1)-1)*a1 (4.0*xi(2)-1)*a2  (4.0*xi(3)-1)*a3 (4.0*xi(4)-1)*a4 ...
4.0*(a1*xi(2)+a2*xi(1)) 4.0*(a2*xi(3)+a3*xi(2)) 4.0*(a1*xi(3)+a3*xi(1)) ...
4.0*(a1*xi(4)+a4*xi(1)) 4.0*(a2*xi(4)+a4*xi(2)) 4.0*(a3*xi(4)+a4*xi(3))];

Nfy=[(4.0*xi(1)-1)*b1 (4.0*xi(2)-1)*b2 (4.0*xi(3)-1)*b3 (4.0*xi(4)-1)*b4 ...
4.0*(b1*xi(2)+b2*xi(1)) 4.0*(b2*xi(3)+b3*xi(2)) 4.0*(b1*xi(3)+b3*xi(1)) ...
4.0*(b1*xi(4)+b4*xi(1)) 4.0*(b2*xi(4)+b4*xi(2)) 4.0*(b3*xi(4)+b4*xi(3))];

Nfz=[(4.0*xi(1)-1)*c1 (4.0*xi(2)-1)*c2 (4.0*xi(3)-1)*c3 (4.0*xi(4)-1)*c4 ...
4.0*(c1*xi(2)+c2*xi(1)) 4.0*(c2*xi(3)+c3*xi(2)) 4.0*(c1*xi(3)+c3*xi(1)) ...
4.0*(c1*xi(4)+c4*xi(1)) 4.0*(c2*xi(4)+c4*xi(2)) 4.0*(c3*xi(4)+c4*xi(3))];

B=[Nfx(1) 0 0 Nfx(2) 0 0 Nfx(3) 0 0 Nfx(4) 0 0 Nfx(5) 0 0 Nfx(6) 0 0 Nfx(7) 0 0 Nfx(8) 0 0 Nfx(9) 0 0 Nfx(10) 0 0;
0 Nfy(1) 0 0 Nfy(2) 0 0 Nfy(3) 0 0 Nfy(4) 0 0 Nfy(5) 0 0 Nfy(6) 0 0 Nfy(7) 0 0 Nfy(8) 0 0 Nfy(9) 0 0 Nfy(10) 0;
0 0 Nfz(1) 0 0 Nfz(2) 0 0 Nfz(3) 0 0 Nfz(4) 0 0 Nfz(5) 0 0 Nfz(6) 0 0 Nfz(7) 0 0 Nfz(8) 0 0 Nfz(9) 0 0 Nfz(10);
Nfy(1) Nfx(1) 0 Nfy(2) Nfx(2) 0 Nfy(3) Nfx(3) 0 Nfy(4) Nfx(4) 0 Nfy(5) Nfx(5) 0 Nfy(6) Nfx(6) 0 Nfy(7) Nfx(7) 0 Nfy(8) Nfx(8) 0 Nfy(9) Nfx(9) 0 Nfy(10) Nfx(10) 0;
0 Nfz(1) Nfy(1) 0 Nfz(2) Nfy(2) 0 Nfz(3) Nfy(3) 0 Nfz(4) Nfy(4) 0 Nfz(5) Nfy(5) 0 Nfz(6) Nfy(6) 0 Nfz(7) Nfy(7) 0 Nfz(8) Nfy(8) 0 Nfz(9) Nfy(9) 0 Nfz(10) Nfy(10);
Nfz(1) 0 Nfx(1) Nfz(2) 0 Nfx(2) Nfz(3) 0 Nfx(3) Nfz(4) 0 Nfx(4) Nfz(5) 0 Nfx(5) Nfz(6) 0 Nfx(6) Nfz(7) 0 Nfx(7) Nfz(8) 0 Nfx(8) Nfz(9) 0 Nfx(9) Nfz(10) 0 Nfx(10)];

K=K+weight*(B'*E*B*Jdet/6.0);
Jdet

end


The node ordering is as follows 1 2 3 4 (1+2) (2+3) (1+3) (1+4) (2+4) (3+4) where 1 2 3 4 denote the corner nodes and terms in the bracket (a+b) denote the node which forms the midpoint of the line node a and node b. The ordering of the nodes in the both the set of data is exactly the same