1
$\begingroup$

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

http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ How to check this? The grid generator just places the points at the mid point of the line joining the two nodes. Does the positioning matter even for a constant metric element? $\endgroup$ – kaush May 23 '17 at 9:49
  • $\begingroup$ If the all mid nodes are in the middle of the edges, then all should be all ok, and Jacobian should be equal to 4 node tet. This is how you can verify this. $\endgroup$ – likask May 23 '17 at 13:39
  • $\begingroup$ Yes that is indeed the case for the data given above. Except that there may be some approximation error $\endgroup$ – kaush May 23 '17 at 13:50
  • $\begingroup$ You can do it like that, calculate normal of the face, and check if it pointing into volume side. If not this is problem with ordering. For problematic element, simply, take two edges on the face, calculate cross product and check direction. If direction is ok, then you have problem somewhere with your code or equations. $\endgroup$ – likask May 23 '17 at 14:08
  • $\begingroup$ A bit more elaborate. Volumes and face normals are calculated just considering the corner nodes and point outside the face. The examples taken hand constructs the nodes so the doubt. I will verify it once again giving a different tet coordinate. $\endgroup$ – kaush May 23 '17 at 15:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.