title: I get different signs on the Jacobian's determinant for different Gauss integration points in a square 20 node hexa element, can that be correct? MWE
body: <p>I'm trying to calculate the stiffness- and mass matrix contributions for a 3D case, using 20 node hexahedron elements.</p>
  <p>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 &quot;optimal&quot; reference shape. Sometimes the determinant aligns with the volume of the element, but I don't think that's always true, or?</p>
  <p>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.</p>
  <p>My test case is a 20-by-20-by-20 cube. The shape functions are from Carlos Felippa's book &quot;Advanced Finite element methods. I follow the node ordering from Abaqus.</p>
  <p>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.</p>
  <p><a href="https://i.stack.imgur.com/pAlmP.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/pAlmP.png" alt="shape functions for reference element" /></a></p>
  <p>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.</p>
  <p><a href="https://i.stack.imgur.com/co9Zk.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/co9Zk.png" alt="elements in local and global coordinates" /></a></p>
  <p>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.</p>
  <pre><code>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 &amp; 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...
  </code></pre>
  <p>The output is as follows:</p>
  <pre><code>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
  </code></pre>