计算科学 - tet10 元素上的 FEM：高斯点的负行列式 - 吾爱随笔录

我正在尝试在 tet10 元素上实现 fem 代码。我遵循 tet10 实施的讲义

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

我已经使用提供的示例验证了我的刚度矩阵。当与单位函数积分时，高斯点处雅可比的数值加起来就是单元的体积。

但是，当我采用节点顺序完全相同的一组不同节点时，我得到一个负行列式。我现在不确定如何确定我的代码是否正确。

我提供了与论文中给出的原始数据以及新数据相同的 matlab 代码片段。

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

节点排序如下 1 2 3 4 (1+2) (2+3) (1+3) (1+4) (2+4) (3+4) 其中 1 2 3 4 表示角节点，括号 (a+b) 中的术语表示形成线节点 a 和节点 b 的中点的节点。两组数据中节点的顺序完全相同