Get node ordering of an element

[derne058]

Hello everyone,

I am writing a code where I want to calculate the deformations and stresses at the Gauss points of an element and compare them with the results calculated by APDL.
I am using the isoparametric formulation for hexahedral linear elements and have implemented the corresponding shape functions.
I would need to extrapolate for each element its nodes, i.e. a vector in which in the first position is the node ID of node i, in the second position the node ID of node J, and so on up to index P.
So far I haven't found anything that works for me.
Does anyone have any ideas or has ever done anything similar?

[mikerife]

Hi @derne058 Fields 10-30 have the element nodal connectivity: https://mapdl.docs.pyansys.com/version/stable/api/_autosummary/ansys.mapdl.core.mesh_grpc.MeshGrpc.elem.html

If that is too much we can *get (or get_value in PyMAPDL) the specific node for an element. Like return the 'i' node, or the 'j' node. Where the nodes traditionally are numbered/lettered i, j, k, l, m, n,......
Mike

[derne058]

Hi @mikerife, thanks a lot!

[derne058]

Hi @mikerife, I think I'm able now to get the right node ordering.
After calculating the displacement at integration points (I used the command mapdl.eresx("NO")) I tried to compute the deformation and displacement at those points using the strain-displacement matrix that I built:
{deformation} = [B] {u}
{stress} = [D] [B] {u}.
I found a huge differences between what I compute and what ANSYS compute, but honestly I don't know where I'm doing wrong.
Below the class I implemented to calculate the [B] matrix.

`class ElementFormulation:
def init(self, E, ni, coords):
"""
Initialize the Hexahedral Element with material properties and nodal coordinates
"""
self.E = E # Young's modulus [Pa]
self.ni = ni # Poisson's ratio
self.coords = np.array(coords) # Nodal coordinates (8 nodes)
self.D = self.constitutive_matrix() # Elasticity matrix

    # Gauss quadrature points and weights
    self.gauss_pts = np.array([-1 / np.sqrt(3), 1 / np.sqrt(3)])
    self.weights = np.array([1, 1])

def constitutive_matrix(self):
    """
    Compute the elastic stiffness (constitutive) matrix for an isotropic material
    """
    G = self.E / (2 * (1 + self.ni))    # Shear modulus [Pa]
    D1 = self.E * (1 - self.ni) / ((1 + self.ni) * (1 - 2 * self.ni))    # First term
    D2 = self.ni * self.E / ((1 + self.ni) * (1 - 2 * self.ni))    # Second term
    D3 = G    # Shear modulus terms

    D = np.array([[D1, D2, D2, 0, 0, 0],
                  [D2, D1, D2, 0, 0, 0],
                  [D2, D2, D1, 0, 0, 0],
                  [0, 0, 0, D3, 0, 0],
                  [0, 0, 0, 0, D3, 0],
                  [0, 0, 0, 0, 0, D3]])
    return D

def shape_functions(self, xi, eta, zeta):    # this function can be removed, it doesn't do anything for computing [B] and [Ke]
    """
    Compute the shape functions for the 8-node hexahedral element
    Following the ANSYS convention:
    xi   --> s
    eta  --> t
    zeta --> r
    """
    N = np.zeros(8)
    N[0] = 0.125 * (1 - xi) * (1 - eta) * (1 - zeta)
    N[1] = 0.125 * (1 + xi) * (1 - eta) * (1 - zeta)
    N[2] = 0.125 * (1 + xi) * (1 + eta) * (1 - zeta)
    N[3] = 0.125 * (1 - xi) * (1 + eta) * (1 - zeta)
    N[4] = 0.125 * (1 - xi) * (1 - eta) * (1 + zeta)
    N[5] = 0.125 * (1 + xi) * (1 - eta) * (1 + zeta)
    N[6] = 0.125 * (1 + xi) * (1 + eta) * (1 + zeta)
    N[7] = 0.125 * (1 - xi) * (1 + eta) * (1 + zeta)
    return N

def shape_function_derivatives(self, xi, eta, zeta):
    """
    Compute the derivatives of the shape functions with respect to xi, eta, and zeta (they are scalar values)
    """
    dN_dxi = np.array([-(1 - eta) * (1 - zeta), (1 - eta) * (1 - zeta), (1 + eta) * (1 - zeta), -(1 + eta) * (1 - zeta),
                       -(1 - eta) * (1 + zeta), (1 - eta) * (1 + zeta), (1 + eta) * (1 + zeta), -(1 + eta) * (1 + zeta)]) / 8
    dN_deta = np.array([-(1 - xi) * (1 - zeta), -(1 + xi) * (1 - zeta), (1 + xi) * (1 - zeta), (1 - xi) * (1 - zeta),
                        -(1 - xi) * (1 + zeta), -(1 + xi) * (1 + zeta), (1 + xi) * (1 + zeta), (1 - xi) * (1 + zeta)]) / 8
    dN_dzeta = np.array([-(1 - xi) * (1 - eta), -(1 + xi) * (1 - eta), -(1 + xi) * (1 + eta), -(1 - xi) * (1 + eta),
                         (1 - xi) * (1 - eta), (1 + xi) * (1 - eta), (1 + xi) * (1 + eta), (1 - xi) * (1 + eta)]) / 8
    grad_N = np.vstack((dN_dxi, dN_deta, dN_dzeta))

    return grad_N

def jacobian_matrix(self, grad_N):
    """
    Compute the Jacobian matrix for the transformation from the reference element to the real element
    """
    J = grad_N @ self.coords
    detJ = np.linalg.det(J)
    invJ = np.linalg.inv(J)
    return J, invJ, detJ

def compute_B_matrix(self, invJ, grad_N):
    """
    Compute the strain-displacement matrix B using the inverse Jacobian 
    """
    dN_dxyz = invJ @ grad_N

    B = np.zeros((6, 24))
    for i in range(8):
        B[0, i*3] = dN_dxyz[0, i]      # ∂N_i/∂x (normal strain in x)
        B[1, i*3 + 1] = dN_dxyz[1, i]  # ∂N_i/∂y (normal strain in y)
        B[2, i*3 + 2] = dN_dxyz[2, i]  # ∂N_i/∂z (normal strain in z)

        B[3, i*3] = dN_dxyz[1, i]      # ∂N_i/∂y (shear strain in xy)
        B[3, i*3 + 1] = dN_dxyz[0, i]  # ∂N_i/∂x (shear strain in xy)

        B[4, i*3 + 1] = dN_dxyz[2, i]  # ∂N_i/∂z (shear strain in yz)
        B[4, i*3 + 2] = dN_dxyz[1, i]  # ∂N_i/∂y (shear strain in yz)

        B[5, i*3] = dN_dxyz[2, i]      # ∂N_i/∂z (shear strain in zx)
        B[5, i*3 + 2] = dN_dxyz[0, i]  # ∂N_i/∂x (shear strain in zx)
        
    return B`
    
    
  Then I compute for a specific Gauss point (giving the coordinate xhi, eta, my) the elasticity matrix [D], the shape function derivatives (dN), the Jacobian (J) and so on:

`` element = ElementFormulation(E, ni, coord) D = element.constitutive_matrix() xi = -1/(3)**0.5 eta = -1/(3)**0.5 zeta = -1/(3)**0.5 dN = element.shape_function_derivatives(xi, eta, zeta) J, invJ, detJ = element.jacobian_matrix(dN) B = element.compute_B_matrix(invJ, dN)

I don't know what to check and where the error is.