对于方形 20 节点六边形元素中的不同高斯积分点，我在雅可比行列式上得到不同的符号，这是否正确？移动电源

如何解决对于方形 20 节点六边形元素中的不同高斯积分点，我在雅可比行列式上得到不同的符号，这是否正确？移动电源

我正在尝试使用 20 个节点六面体元素计算 3D 案例的刚度和质量矩阵贡献。

虽然我得到了非常奇怪的结果。对于不同的高斯点（我使用 3×3×3 积分，所以是 27 点），雅可比矩阵到处都是。它可以是正数和超过 1e3，负数低于 1.5e3 或其他东西，有时是 0。我对雅可比矩阵的解释是它将元素转换为参考元素以便于积分，行列式是一种查看多少的方法原始元素偏离“最佳”参考形状。有时行列式与元素的体积一致，但我认为这并不总是正确的，或者？

如果 Jacobian 为负或非常接近于零，则意味着元素失真或节点编号错误。但对我来说，元素和节点的顺序都不是错误的。我对节点编号进行了两次和三次检查，但我的逻辑中仍然可能存在一些小错误或错误。

我的测试用例是一个 20×20×20 的立方体。形状函数来自 Carlos Felippa 的书“高级有限元方法。我遵循 Abaqus 中的节点排序。

形状函数如下，对于下面的参考元素。我使用 Python 中的符号工具箱集成了形状函数，请参阅下面的 MWE。

下图只是我的参考元素和全局坐标中元素的小图。它显示元素轮廓，蓝色点代表 20 个节点，黑色穿过 27 个积分点。

下面的代码是一个最小的工作示例，它为 27 个高斯点中的每一个打印雅可比行列式。有人能发现问题吗？虽然阅读别人的代码总是很困难，但我认为如果您了解 FEM，这个案例就很简单了。

import sympy as s
import numpy as np


xhi = s.Symbol('xhi')
eta = s.Symbol('eta')
my  = s.Symbol('my')
# node 1-20 isoparametric values
#         1   2   3   4   5   6  7   8   9  10   11   12  13   14   15   16    17  18   19    20 
xhi_i = np.array((-1,1,-1,-1)) 
eta_i = np.array((-1,1))
my_i  = np.array((-1,0))   

# Coordinates of my element in the global xyz coordinate system
coord = np.array(((20*xhi_i),(20*eta_i),(20*my_i))).transpose()

v = 0.3
E = 210e3
D = np.array(((1-v,v,0),(   v,1-v,(   0,(1-2*v)/2,(1-2*v)/2)))
D *= E/((1+v)*(1-2*v))



# Define and derive shape functions,and save a list
N = []
# list with shape function derivates in each direction:
dNdXhi,dNdEta,dNdMy = [],[],[]
for i in range(0,20):
    if i+1 in [1,2,3,4,5,6,7,8]:
        N_i = (1/8)*(1+xhi*xhi_i[i])*(1+eta*eta_i[i])*(1+my*my_i[i])*((my*my_i[i]) + (eta*eta_i[i]) + (my*my_i[i]) - 2)
    if i+1 in [9,11,13,15]:
        N_i = (1/4)*(1-(xhi*xhi))*(1+(eta*eta_i[i]))*(1+(my*my_i[i]))        
    if i+1 in [10,12,14,16]:
        N_i = (1/4)*(1-(eta*eta))*(1+(xhi*xhi_i[i])*(1+(my*my_i[i])))
    if i+1 in [17,18,19,20]:
        N_i = (1/4)*(1-(my*my))*(1+(xhi*xhi_i[i])*(1+(eta*eta_i[i])))
    # derive using symbolic toolBox
    dNdXhi.append(s.diff(N_i,xhi))
    dNdEta.append(s.diff(N_i,eta))
    dNdMy.append(s.diff(N_i,my))

shape_functions_derivates = [dNdXhi,dNdMy]

# loop over 27 Gauss points:
# i,j & k are the points where we evaluate the stiffness contribution
for i in [-0.774596669241483,0.0,0.774596669241483]:
    for j in [-0.774596669241483,0.774596669241483]:
        for k in [-0.774596669241483,0.774596669241483]:
            if i==0.0 or j==0.0 or k==0.0:
                w = 0.888888888888889
            else:
                w= 0.555555555555556
            # evaluate shape function derivates
            dNdXhidEtadmy = np.zeros((3,20))
            for node in range(0,19):
                dNdXhidEtadmy[0,node] = dNdXhi[node].evalf(subs={xhi:i,eta:j,my:k})
                dNdXhidEtadmy[1,node] = dNdEta[node].evalf(subs={xhi:i,my:k})
                dNdXhidEtadmy[2,node] = dNdMy[node].evalf(subs={xhi:i,my:k})
            #  @ means matrix multiplication
            Jacobian = dNdXhidEtadmy@coord
            detJ = np.linalg.det(Jacobian)
            print('detJ='+str(detJ)+',for xhi='+str(i)+',eta='+str(j)+',my='+str(k))
            # Find stiffness- and mass matrix later by forming B matrix and some multiplications...

输出如下：

detJ=-17855.99999999991,for xhi=-0.77459,eta=-0.77459,my=-0.77459
detJ=-3.197442310920449e-13,my=0.0
detJ=-17855.99999999994,my=0.77459
detJ=8063.999999999968,eta=0.0,my=-0.77459
detJ=8.758115402030098e-47,my=0.0
detJ=8063.999999999968,my=0.77459
detJ=-17855.99999999994,eta=0.77459,my=-0.77459
detJ=3.197442310920437e-13,my=0.77459
detJ=-4607.999999999974,for xhi=0.0,my=-0.77459
detJ=1.2621774483536082e-29,my=0.0
detJ=-4607.999999999974,my=0.77459
detJ=5759.999999999983,my=-0.77459
detJ=0.0,my=0.0
detJ=5759.999999999983,my=-0.77459
detJ=1.2621774483536262e-29,for xhi=0.77459,my=-0.77459
detJ=1.385558334732187e-12,my=0.77459