I'm trying to calculate the stiffness- and mass matrix contributions for a 3D case, using 20 node hexahedron elements.
I get very strange results though. The Jacobian is all over the place for the different Gauss points (I use 3-by-3-by-3 integration, so 27 points). It can be positive and over 1e3, and negative below 1.5e3 or something, and sometimes 0. My interpretation of the Jacobian is that it transforms the element to the reference element for easier integration, and that the determinant is a way to see how much the original element stray from the "optimal" reference shape. Sometimes the determinant aligns with the volume of the element, but I don't think that's always true, or?
If the Jacobian is negative or very close to zero, it should mean the element is distorted or that the node numbering is faulty. But to me neither the element nor the node ordering is wrong. I have double and triple checked the node numbering, but there could still be some small bug or error in my logic.
My test case is a 20-by-20-by-20 cube. The shape functions are from Carlos Felippa's book "Advanced Finite element methods. I follow the node ordering from Abaqus.
The shape functions are as follows, for the reference element below. I integrate the shape functions using the symbolic toolbox in Python, see MWE below.
The figure below is just a small graph of my reference element and element in global coordinates. It shows the element outlines, the blue dots the 20 nodes, and the black crosses the 27 integrations points.
The code below is a minimal working example that prints the determinant of the Jacobian for each of the 27 Gauss points. Could anyone spot something wrong? While it's always hard to read someone else's code I think this case is pretty straight forward if you know 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, -1, 1, 1, -1, 0, 1, 0, -1, 0, 1, 0, -1, -1, 1, 1, -1))
eta_i = np.array((-1, -1, 1, 1, -1, -1, 1, 1, -1, 0, 1, 0, -1, 0, 1, 0 , -1, -1, 1, 1))
my_i = np.array((-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1 , 0, 0, 0, 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, v, 0, 0, 0),
( v, 1-v, v, 0, 0, 0),
( v, v, 1-v, 0, 0, 0),
( 0, 0, 0, (1-2*v)/2, 0, 0),
( 0, 0, 0, 0, (1-2*v)/2, 0),
( 0, 0, 0, 0, 0, (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,dNdEta,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.0, 0.774596669241483]:
for k in [-0.774596669241483, 0.0, 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,eta:j,my:k})
dNdXhidEtadmy[2,node] = dNdMy[node].evalf(subs={xhi:i,eta:j,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...
The output is as follows:
detJ=-17855.99999999991, for xhi=-0.77459, eta=-0.77459, my=-0.77459
detJ=-3.197442310920449e-13, for xhi=-0.77459, eta=-0.77459, my=0.0
detJ=-17855.99999999994, for xhi=-0.77459, eta=-0.77459, my=0.77459
detJ=8063.999999999968, for xhi=-0.77459, eta=0.0, my=-0.77459
detJ=8.758115402030098e-47, for xhi=-0.77459, eta=0.0, my=0.0
detJ=8063.999999999968, for xhi=-0.77459, eta=0.0, my=0.77459
detJ=-17855.99999999994, for xhi=-0.77459, eta=0.77459, my=-0.77459
detJ=3.197442310920437e-13, for xhi=-0.77459, eta=0.77459, my=0.0
detJ=-17855.99999999994, for xhi=-0.77459, eta=0.77459, my=0.77459
detJ=-4607.999999999974, for xhi=0.0, eta=-0.77459, my=-0.77459
detJ=1.2621774483536082e-29, for xhi=0.0, eta=-0.77459, my=0.0
detJ=-4607.999999999974, for xhi=0.0, eta=-0.77459, my=0.77459
detJ=5759.999999999983, for xhi=0.0, eta=0.0, my=-0.77459
detJ=0.0, for xhi=0.0, eta=0.0, my=0.0
detJ=5759.999999999983, for xhi=0.0, eta=0.0, my=0.77459
detJ=-4607.999999999974, for xhi=0.0, eta=0.77459, my=-0.77459
detJ=1.2621774483536262e-29, for xhi=0.0, eta=0.77459, my=0.0
detJ=-4607.999999999974, for xhi=0.0, eta=0.77459, my=0.77459
detJ=-17855.99999999994, for xhi=0.77459, eta=-0.77459, my=-0.77459
detJ=1.385558334732187e-12, for xhi=0.77459, eta=-0.77459, my=0.0
detJ=-17855.99999999994, for xhi=0.77459, eta=-0.77459, my=0.77459
detJ=8063.999999999968, for xhi=0.77459, eta=0.0, my=-0.77459
detJ=0.0, for xhi=0.77459, eta=0.0, my=0.0
detJ=8063.999999999968, for xhi=0.77459, eta=0.0, my=0.77459
detJ=-17855.99999999994, for xhi=0.77459, eta=0.77459, my=-0.77459
detJ=3.197442310920437e-13, for xhi=0.77459, eta=0.77459, my=0.0
detJ=-17855.99999999994, for xhi=0.77459, eta=0.77459, my=0.77459
