How to assemble K local matrix for a hexahedron cell in a piecewise linear fashion?

Question

import meshio
from funcs import hexaBasis 
import sympy as sp 
from sympy import diff, symbols 
from scipy.optimize import root
import numpy as np 

x, y, z = symbols('x y z')
# Define quadrature points and weights for a hexahedron (2x2x2 Gaussian quadrature)
quadrature_points = [
    [-1/np.sqrt(3), -1/np.sqrt(3), -1/np.sqrt(3)],
    [ 1/np.sqrt(3), -1/np.sqrt(3), -1/np.sqrt(3)],
    [ 1/np.sqrt(3),  1/np.sqrt(3), -1/np.sqrt(3)],
    [-1/np.sqrt(3),  1/np.sqrt(3), -1/np.sqrt(3)],
    [-1/np.sqrt(3), -1/np.sqrt(3),  1/np.sqrt(3)],
    [ 1/np.sqrt(3), -1/np.sqrt(3),  1/np.sqrt(3)],
    [ 1/np.sqrt(3),  1/np.sqrt(3),  1/np.sqrt(3)],
    [-1/np.sqrt(3),  1/np.sqrt(3),  1/np.sqrt(3)]
]
weights = [1.0] * 8  # All weights are 1 for 2x2x2 Gaussian quadrature
mesh = meshio.read("cube_of_hexahedrons_dolfinx.xdmf")
cell_block0 = mesh.cells[0]
cell0 = cell_block0.data[0]
nodes_cell0 = mesh.points[cell0]

global tab_points
tab_points = [
    [0, 0, 0],  # Node 0
    [1, 0, 0],  # Node 1
    [1, 1, 0],  # Node 2
    [0, 1, 0],  # Node 3
    [0, 0, 1],  # Node 4
    [1, 0, 1],  # Node 5
    [1, 1, 1],  # Node 6
    [0, 1, 1]   # Node 7
]

def global_to_local(x, y, z, tab_points):
    def func(xi_eta_zeta):
        xi, eta, zeta = xi_eta_zeta
        N = [
            (1 - xi) * (1 - eta) * (1 - zeta) / 8,
            (1 + xi) * (1 - eta) * (1 - zeta) / 8,
            (1 + xi) * (1 + eta) * (1 - zeta) / 8,
            (1 - xi) * (1 + eta) * (1 - zeta) / 8,
            (1 - xi) * (1 - eta) * (1 + zeta) / 8,
            (1 + xi) * (1 - eta) * (1 + zeta) / 8,
            (1 + xi) * (1 + eta) * (1 + zeta) / 8,
            (1 - xi) * (1 + eta) * (1 + zeta) / 8,
        ]
        x_calc = sum(N[i] * tab_points[i][0] for i in range(8))
        y_calc = sum(N[i] * tab_points[i][1] for i in range(8))
        z_calc = sum(N[i] * tab_points[i][2] for i in range(8))
        return [x_calc - x, y_calc - y, z_calc - z]

    xi_eta_zeta_initial_guess = [0, 0, 0]
    solution = root(func, xi_eta_zeta_initial_guess)
    return solution.x


def jacobian(Φs, nodes_local, hexahedron_nodes):
    ML = sp.zeros(3, 8)
    MR = sp.zeros(8,3)

    for i in range(8):
        ML[0,i] = diff(Φs[i], x).subs({x: nodes_local[i][0], y: nodes_local[i][1], z: nodes_local[i][2]})
        ML[1,i] = diff(Φs[i], y).subs({x: nodes_local[i][0], y: nodes_local[i][1], z: nodes_local[i][2]})
        ML[2,i] = diff(Φs[i], z).subs({x: nodes_local[i][0], y: nodes_local[i][1], z: nodes_local[i][2]})

    for i in range(3):
        MR[0,i] = hexahedron_nodes[0][i]
        MR[1,i] = hexahedron_nodes[1][i]
        MR[2,i] = hexahedron_nodes[2][i]
        MR[3,i] = hexahedron_nodes[3][i]
        MR[4,i] = hexahedron_nodes[4][i]
        MR[5,i] = hexahedron_nodes[5][i]
        MR[6,i] = hexahedron_nodes[6][i]
        MR[7,i] = hexahedron_nodes[7][i]

    return ML * MR

Φs = hexaBasis()
gradΦs = []
NBs = []



nodes_local =[]
for i in range(8):
    nodes_local.append( global_to_local(nodes_cell0[i][0], nodes_cell0[i][1], nodes_cell0[i][2], tab_points ) ) 


for i in range(8):
    gradΦ0 = diff(Φs[0],x).subs({x: nodes_local[i][0], y: nodes_local[i][1], z: nodes_local[i][2], }) 
    + diff(Φs[0],y).subs({x: nodes_local[i][0], y: nodes_local[i][1], z: nodes_local[i][2], }) 
    + diff(Φs[0],z).subs({x: nodes_local[i][0], y: nodes_local[i][1], z: nodes_local[i][2], })
    gradΦs.append(gradΦ0)


M = sp.zeros(8)
for i in range(8):
    M[i,0] = gradΦs[i] * gradΦs[0]
    M[i,1] = gradΦs[i] * gradΦs[1]
    M[i,2] = gradΦs[i] * gradΦs[2]
    M[i,3] = gradΦs[i] * gradΦs[3]
    M[i,4] = gradΦs[i] * gradΦs[4]
    M[i,5] = gradΦs[i] * gradΦs[5]
    M[i,6] = gradΦs[i] * gradΦs[6]
    M[i,7] = gradΦs[i] * gradΦs[7]


J = jacobian(Φs, nodes_local, nodes_cell0)

M = M * J.det()


print(M)
print('done...')

'''
[1.93359375000000, 1.93359375000000, 1.45019531250000, 1.45019531250000, 1.45019531250000, 1.45019531250000, 1.08764648437500, 1.08764648437500]
'''

'''
dolfinx output of the same K local from the same mesh of cell
[[ 8.33333333e-02 -3.95818496e-18  2.67916521e-18 ...  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [-5.22280515e-18  1.66666667e-01 -2.08333333e-02 ...  0.00000000e+00
   0.00000000e+00  0.00000000e+00]

'''

There seems to be something not quite right yet. The local K matrix from the same cell0 for the same mesh xdmf that was imported here is reporting different values than I have. Seems that things maybe need to be integrated some way. Ideally I need to match the dolfinx values within a small margin. So far I attempted to construct the hexahedron in a piecewise linear fashion. Why are my results so far away from what dolfinx says for the same cell?

Is there a way to integrate the gradients over the domain?

Answer 1

Φs = hexaBasis()
gradΦs = []

nodes_local = []
for i in range(8):
    nodes_local.append(global_to_local(nodes_cell0[i][0], nodes_cell0[i][1], nodes_cell0[i][2], tab_points))

for i in range(8):
    gradΦ0 = diff(Φs[0], x).subs({x: quadrature_points[i][0], y: quadrature_points[i][1], z: quadrature_points[i][2]}) 
    + diff(Φs[0], y).subs({x: quadrature_points[i][0], y: quadrature_points[i][1], z: quadrature_points[i][2]}) 
    + diff(Φs[0], z).subs({x: quadrature_points[i][0], y: quadrature_points[i][1], z: quadrature_points[i][2]})
    gradΦs.append(gradΦ0)

M_k = sp.zeros(8)
for i in range(8):
    M_k[i, 0] = gradΦs[i] * gradΦs[0]
    M_k[i, 1] = gradΦs[i] * gradΦs[1]
    M_k[i, 2] = gradΦs[i] * gradΦs[2]
    M_k[i, 3] = gradΦs[i] * gradΦs[3]
    M_k[i, 4] = gradΦs[i] * gradΦs[4]
    M_k[i, 5] = gradΦs[i] * gradΦs[5]
    M_k[i, 6] = gradΦs[i] * gradΦs[6]
    M_k[i, 7] = gradΦs[i] * gradΦs[7]

# Integrate each element of M_k symbolically using SymPy
M_k_integrated = sp.zeros(8, 8)

for i in range(8):
    for j in range(8):
        M_k_integrated[i, j] = M_k[i, j].integrate((x, x_min, x_max), (y, y_min, y_max), (z, z_min, z_max))

# Convert elements of M_k_integrated to floating-point numbers
M_k_integrated_decimal = M_k_integrated.evalf()

# Print the integrated M_k in decimal format
print("Integrated M_k:")
for row in M_k_integrated_decimal.tolist():
    print([float(elem) for elem in row])

This is so far a tentative answer to the question of how to integrate the system. It provides values similar to other FEM software however the results are not identical. So far is there is anything incorrect here I am not aware of it at this time.

The limits of integration are the limits of the element in the global reference frame.