Which vertex is where?

Hi,

I have a 2d domain (actually a cut of phasespace that resolves one spatial coordinate “z” and one velocity coordinate “vz”) and need to set the electric field Ez at all grid points with the result of a separate calculation. Physically the electric field Ez is a function of z only, i.e. the value should be the same at all vz for a given z. I tried to set Ez as in the following code, but the results are not as expected:

#!/usr/bin/env python3
  
from fenics import *
from petsc4py import PETSc
import argparse
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
from scipy.special import lambertw


parser = argparse.ArgumentParser(description='1d1v twofluid with electric field using FEM')
parser.add_argument('--degree', dest='degree', action='store', default=2, type=int, help='degree of base elements (default: 2)')
parser.add_argument('--Nv', dest='Nv', action='store', default=67, type=int, help='number of grid cells (default: 67)')
parser.add_argument('--Nz', dest='Nz', action='store', default=37, type=int, help='number of grid cells (default: 37)')
parser.add_argument('--Lz', dest='Lz', action='store', default=3.0, type=float, help='domain length (default: 3.0)')
parser.add_argument('--vmax', dest='vmax', action='store', default=3.0, type=float, help='maximum velocity (default: 3.0)')
parser.add_argument('--vthe', dest='vthe', action='store', default=0.3, type=float, help='electron thermal velocity (default: 0.3)')

args = parser.parse_args()
globals().update(args.__dict__)

# convert mesh to triangles in a format that matplotlib understands
def mesh2triang(mesh):
    xy = mesh.coordinates()
    return tri.Triangulation(xy[:, 0], xy[:, 1], mesh.cells())


# 2d phasespace
phasespace = RectangleMesh(Point(0., -vmax*vthe), Point(Lz, vmax*vthe), Nz-1, Nv-1, "left/right")

# define function space for phasespace density
P2 = FiniteElement('P', 'triangle', degree)
Ve2d = FunctionSpace(phasespace, P2)

# lower order for plotting
Ve2d_deg1 = FunctionSpace(phasespace, FiniteElement('P', 'triangle', 1))

# Electric field controlled from the outside
Ez = project(Constant(0.), Ve2d_deg1)
Ezarray = Ez.vector().vec().array.reshape( Nz, Nv )
for i in range(Nz):
    for j in range(Nv):
        Ezarray[i,j] = i
Ez_e = interpolate(Ez, Ve2d)

Cz = Ez_e.compute_vertex_values(phasespace)

# Plots stuff
Z = np.linspace(Lz, 0., Nz)
pz = plt.tripcolor(mesh2triang(phasespace), Cz, shading='gouraud')
cb3 = plt.colorbar(pz, label=r"$E_z$")

plt.tight_layout()
plt.show()
plt.close()

In the above example I simply use i as the electric field at each z, in a more realistic example this would be f(i). But either way the result in the plot is NOT vertical stripes (i.e. values independent of vz), but a complicated pattern.

So the question is: how do I figure out which element in Ezarray corresponds to which z,vz? Or alternatively how do I put a numpy array Ezvec[z] onto the 2d phase space?