3D六边形坐标扩展：从8个坐标生成27个坐标及代码调试

我现有8个3D坐标点，需要在它们之间插值生成新坐标，最终扩展为27个坐标点（对应3×3×3的六边形网格结构）。

我尝试编写了一段Python代码，但运行效果不理想：

import numpy as np

xi = [-1.0, 0.0, 1.0]
# ele_arr 为给定的8个坐标组成的数组

cord1 = np.zeros([(N+1)*(N+1)*(NPZ+1),4],dtype = 'float64')

cord1[0,:] = ele_arr[-1,:]
cord1[N,:] = ele_arr[-2,:]
cord1[N*(N+1),:] = ele_arr[3,:]
cord1[((N+1)*(N+1))-1,:] = ele_arr[2,:]
cord1[-1,:] = ele_arr[0,:]
cord1[-1-(N),:] = ele_arr[1,:]
cord1[-1-(N*(N+1)),:] = ele_arr[-4,:]
cord1[-((N+1)*(N+1)),:] = ele_arr[-3,:]


for i in range(1,NPZ):
    gap = (N+1) * (N+1)
    cord1[i*gap,1:4] = np.array([cord1[0,1] + ((cord1[-((N+1)*(N+1)),1] - cord1[0,1]) * ((xi[i] - xi[0])/2)),
                                 cord1[0,2] + ((cord1[-((N+1)*(N+1)),2] - cord1[0,2]) * ((xi[i] - xi[0])/2)),
                                 cord1[0,3] + ((cord1[-((N+1)*(N+1)),3] - cord1[0,3]) * ((xi[i] - xi[0])/2))])

    cord1[i*gap + N,1:4] = np.array([cord1[N,1] + ((cord1[-1-(N*(N+1)),1] - cord1[N,1]) * ((xi[i] - xi[0])/2)),
                                     cord1[N,2] + ((cord1[-1-(N*(N+1)),2] - cord1[N,2]) * ((xi[i] - xi[0])/2)),
                                     cord1[N,3] + ((cord1[-1-(N*(N+1)),3] - cord1[N,3]) * ((xi[i] - xi[0])/2))])

    cord1[i*gap + N + (N+1)*N,1:4] = np.array([cord1[((N+1)*(N+1))-1,1] + ((cord1[-1,1] - cord1[((N+1)*(N+1))-1,1]) * ((xi[i] - xi[0])/2)),
                                           cord1[((N+1)*(N+1))-1,2] +  ((cord1[-1,2] - cord1[((N+1)*(N+1))-1,2]) * ((xi[i] - xi[0])/2)),
                                           cord1[((N+1)*(N+1))-1,3] + ((cord1[-1,3] - cord1[((N+1)*(N+1))-1,3]) * ((xi[i] - xi[0])/2))])

    cord1[i*gap + (N+1)*N,1:4] = np.array([ cord1[N*(N+1),1] + ((cord1[-1-N,1] - cord1[N*(N+1),1]) * ((xi[i] - xi[0])/2)),
                                               cord1[N*(N+1),2] + ((cord1[-1-N,2] - cord1[N*(N+1),2]) * ((xi[i] - xi[0])/2)),
                                               cord1[N*(N+1),3] + ((cord1[-1-N,3] - cord1[N*(N+1),3]) * ((xi[i] - xi[0])/2))])


for i in range(NPZ+1):
    for j in range(1,N):
        cord1[i*(N+1)*(N+1)+j,1:4] = [cord1[i*(N+1)*(N+1),1] + ((cord1[i*(N+1)*(N+1)+N,1] - cord1[i*(N+1)*(N+1),1]) * (((xi[j] - xi[0])/2))),
                                    cord1[i*(N+1)*(N+1),2] + ((cord1[i*(N+1)*(N+1)+N,2] - cord1[i*(N+1)*(N+1),2]) * (((xi[j] - xi[0])/2))),
                                    cord1[i*(N+1)*(N+1),3] + ((cord1[i*(N+1)*(N+1)+N,3] - cord1[i*(N+1)*(N+1),3]) * (((xi[j] - xi[0])/2)))]

        cord1[i*(N+1)*(N+1)+(N*(N+1))+j,1:4] = [cord1[i*(N+1)*(N+1)+(N*(N+1)),1] + ((cord1[i*(N+1)*(N+1)+(N*(N+1))+N,1] - cord1[i*(N+1)*(N+1)+(N*(N+1)),1]) * (((xi[j] - xi[0])/2))),
                                                cord1[i*(N+1)*(N+1)+(N*(N+1)),2] + ((cord1[i*(N+1)*(N+1)+(N*(N+1))+N,2] - cord1[i*(N+1)*(N+1)+(N*(N+1)),2]) * (((xi[j] - xi[0])/2))),
                                                cord1[i*(N+1)*(N+1)+(N*(N+1)),3] + ((cord1[i*(N+1)*(N+1)+(N*(N+1))+N,3] - cord1[i*(N+1)*(N+1)+(N*(N+1)),3]) * (((xi[j] - xi[0])/2)))]

for i in range(NPZ+1):
    for j in range(N+1):
        for k in range(1,N):
            cord1[(i*gap)+(k*(N+1))+j,1:4] = [cord1[(i*gap)+j,1] + ((cord1[(i*gap)+j+(N*(N+1)),1] - cord1[(i*gap)+j,1]) * (((xi[k] - xi[0])/2))),
                                              cord1[(i*gap)+j,2] + ((cord1[(i*gap)+j+(N*(N+1)),2] - cord1[(i*gap)+j,2]) * (((xi[k] - xi[0])/2))),
                                              cord1[(i*gap)+j,3] + ((cord1[(i*gap)+j+(N*(N+1)),3] - cord1[(i*gap)+j,3]) * (((xi[k] - xi[0])/2)))]

请帮忙优化这段代码，实现从8个坐标生成27个3D六边形网格坐标的需求。

内容的提问来源于stack exchange，提问作者SAURABH AGARWAL