I am trying to implement a fem code on tet10 elements.  I follow the lecture notes for tet10 implementation given inhttp://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdfI 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 0node = [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 sameend;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 referenceY=480; %Youngs Modulusnu=1/3; % Poisson ratioalfa=(5.0+3.0*sqrt(5.0))/20; %Gauss Pointsbeta =(5.0-sqrt(5.0))/20;    %Gauss Pointsweight=0.25; %weights for integration pointsGaussPoints=[alfa beta beta beta;beta alfa beta beta;beta beta alfa  beta;beta beta beta alfa];%Elasticity MatrixE=[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 MatrixK=zeros(30,30);for j=1:4xi=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);JdetendThe 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 
FEM on tet10 element: negetive determinant at the Gauss point
finite element;jacobian
For 4 node tetrahedron if you take any face and calculate normal of it has to point in the direction of opposite to this face node then a determination of the Jacobian is positive.In the case of 10 node tetrahedron is more complicated, for example, if you move mid nodes on the edges, to quarter point distance, your Jacobian become singular (this is used in fracture mechanics at the crack tip to approximate singularity). If you can beyond quarter point and then  determinant of  Jacobian will be negative.In other words, not only ordering of the nodes is important but also the position of mid nodes. You have to check this.