I looked into both topics, but I’m still uncertain about how to proceed. I’ll try to explain my issue using a simpler example.

Let’s consider a 2D mesh with 16 vertices and a vector-valued finite element function space of first order

```
from mpi4py import MPI
import dolfinx
import basix.ufl
# mesh
quadrilateral = dolfinx.mesh.CellType.quadrilateral
msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 3, 3, cell_type=quadrilateral)
num_vertices = msh.topology.index_map(0).size_local
num_cells = msh.topology.index_map(msh.topology.dim).size_local
# function space
FE = basix.ufl.element("Lagrange", msh.basix_cell(), degree=1, shape=(msh.geometry.dim,))
V = dolfinx.fem.functionspace(msh, FE)
```

And then I have the external data for each vertex in the following format

```
# x y z f_x f_y f_z
0.33 0.33 0 2.2 1 0
0 0.66 0 3.1 1 0
0.66 0 0 1.3 1 0
1 0.33 0 2.4 1 0
1 1 0 4.4 1 0
0.66 0.33 0 2.3 1 0
0 0 0 1.1 1 0
0.66 0.66 0 3.3 1 0
1 0.66 0 3.4 1 0
0 0.33 0 2.1 1 0
0.33 0 0 1.2 1 0
0 1 0 4.1 1 0
0.66 1 0 4.3 1 0
1 0 0 1.4 1 0
0.33 0.66 0 3.2 1 0
0.33 1 0 4.2 1 0
```

which I read as

```
import numpy as np
data = np.loadtxt("file.dat")[:, 3:]
```

The values for each vertex are in random order and my goal is to assign them to

```
source = fem.Function(V)
source.x.array[:] = data[?, :mesh.geometry.dim].flatten()
```

accordingly, so the coordinates of external values match the coordinates of DOFs.

So from application-of-point-forces-mapping-vertex-indices-to-corresponding-dofs I tried

```
def vertex_to_dof_mapping(mesh, V):
vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
c2v = msh.topology.connectivity(msh.topology.dim, 0)
for cell in range(num_cells):
vertices = c2v.links(cell)
dofs = V.dofmap.cell_dofs(cell)
for i, vertex in enumerate(vertices):
vertex_to_dof_map[vertex] = dofs[V.dofmap.dof_layout.entity_dofs(0, i)]
return vertex_to_dof_map
vertex_to_dof_map = vertex_to_dof_mapping(msh, V)
```

This i guess should answer the dof->vertex mapping, however in my case this returns an array from 0 to 15, without any changes in the ordering?

And in case of vertex->input geometry mapping, thats where I just cycle through the external file and find the corresponding indices, where the coordinates match with *mesh.geometry.x*, right?