I want to share what is now working much faster for me than my old approach via the coordinates.

In preprocessing (serial), read in your mesh with the standard XDMF routine, set values to a function, and write a text file using

`igi = msh.geometry.input_global_indices`

with the following structure

```
igi1 fx1 fy1 fz1
igi2 fx2 fy2 fz2
...
ign fxn fyn fzn
```

so the input node index followed by the x,y,z values of your vector field (or the one value of a scalar field).

Then, in parallel, use

```
def readfunction(f, V, datafile, msh):
# block size of vector
bs = f.vector.getBlockSize()
# load data and input node indices
data = np.loadtxt(datafile,usecols=(np.arange(1,bs+1)),ndmin=2)
ind_file = np.loadtxt(datafile,usecols=(0),dtype=int)
# index map and input indices
im = np.asarray(V.dofmap.index_map.local_to_global(np.arange(V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts, dtype=np.int32)), dtype=PETSc.IntType)
igi = msh.geometry.input_global_indices
# since in parallel, the ordering of the dof ids might change, so we have to find the
# mapping between original and new id via the coordinates
ci = 0
for i in im:
ind = np.where(ind_file == igi[ci])[0]
# only write if we've found the index
if len(ind):
for j in range(bs):
f.vector[bs*i+j] = data[ind[0],j]
ci+=1
f.vector.assemble()
f.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
```