Ordering of basis functions and cell dofs

General
1

I have small subroutine to assemble an interpolation matrix for a mixed Lagrange element with P3/P2 such that,

W0 = W.sub(0)
V, _ = W0.collapse()

# Assemble interpolation matrix.
N_dofs = W.dofmap.index_map.size_global #  Number of velocity + pressure dofs.
N_dofs_v = V.dofmap.index_map.size_global*3 # Number of velocity dofs.
basis_functions_v = W0.element.basix_element.tabulate(0, xi_obs_inside)[0,:,:,0] # Basis functions for velocity evaluated at observation points.
N_basis_functions_v = np.shape(basis_functions_v)[1] # Number of distinct velocity basis functions.
M = PETSc.Mat().createAIJ(size = (N_obs_inside*3, N_dofs), nnz=N_basis_functions_v) 
M.setUp()

for i in range(N_obs_inside):
    cell = obs_inside_2_cell[i]
    block_dofs = W0.dofmap.cell_dofs(cell)
    for j in range(3):
        M.setValues(3*i+j, block_dofs[j::3], basis_functions_v[i,:], PETSc.InsertMode.INSERT_VALUES)

M.assemble()

I am worried that the order of block_dofs[j::3] and basis_functions_v[i,:] do not match. I received False and True for print(Vh.element.basix_element.dof_transformations_are_identity) and
print(Vh.element.basix_element.dof_transformations_are_permutations), respectively. Consequently, I believe they match upto a tranformation. I am lost on how to integrate this. Thank you.

2

Please make the code reproducible.
You reference Vh, which is not defined above.
Use a unit square or cube as the domain.