def METHOD_NAME(model: BDF, pid: int,
                             x: np.ndarray,
                             y: np.ndarray,
                             z: np.ndarray,
                             nx: int, ny: int, nz: int,
                             eid: int=1, nid: int=1) -> tuple[int, int]:
    """
    Parameters
    ----------
    p1 / p2 / p3 / p4 : (3, ) float ndarray
        points defining the quad
    nx : int
        points in the p1-p2 direction
    ny : int
        points in the p1-p4 direction
    nid / eid : int
        node / element id offset
    Returns
    -------
    nid : int
        ???
    eid : int
        ???
    """
    nid0 = nid
    nnodes = nx * ny * nz
    nelements = (nx - 1) * (ny - 1) * (nz - 1)
    assert nelements > 0, f'nx={nx} ny={ny} nz={nz} nelements={nelements}'
    #points, elements = points_elements_from_cube_points(
        #p1, p2, p3, p4,
        #p5, p6, p7, p8,
        #x, y, z, dtype='int32')
    xv, yv, zv = np.meshgrid(x, y, z)
    for point in zip(xv.ravel(), yv.ravel(), zv.ravel()):
        model.add_grid(nid, point)
        nid += 1
    node_ids = np.arange(nnodes, dtype='int32').reshape(nx, ny, nz)
    #print(node_ids)
    p1 = node_ids[:nx-1, :ny-1, :nz-1].ravel()
    p2 = node_ids[1:,    :ny-1, :nz-1].ravel()
    p3 = node_ids[1:,       1:, :nz-1].ravel()
    p4 = node_ids[:nx-1,    1:, :nz-1].ravel()
    p5 = node_ids[:nx-1, :ny-1, 1:].ravel()
    p6 = node_ids[1:,    :ny-1, 1:].ravel()
    p7 = node_ids[1:,       1:, 1:].ravel()
    p8 = node_ids[:nx-1,    1:, 1:].ravel()
    elements = np.vstack([p1, p2, p3, p4, p5, p6, p7, p8]).T
    #print(elements)
    #print(elements.shape)
    for node_ids in (elements + nid0).tolist():
        #print(node_ids, type(node_ids))
        model.add_chexa(eid, pid, node_ids)
        eid += 1
    return nid, eid