Hi all,

I have a code in dolfinx that solves a system of PDEs in a rather complex geometry that I mesh using gmsh. The result is a vector field \mathbf m. I would now need to calculate some quantities for each cell of the mesh, for which I need to identify all the vertices that belong to them. Graphically, it would be something like this:

```
0 -------- 1 -------- 4
| | |
| cell 0 | cell 1 |
| | |
2 -------- 3 -------- 5
```

I would need to know which vertices (in this case, 0, 1, 2, 3) belong to cell 0, calculate some stuff with the value of \mathbf m at them, then do the same for cell 1 (with vertices 1, 4, 5, 3), etc.

Can anyone give me a hint on how to achieve this?

Thanks in advance!

Here’s a MWE on how I create the mesh:

```
import gmsh
model = gmsh.model; geom = gmsh.model.geo
from mpi4py import MPI
mesh_comm = MPI.COMM_WORLD
ms = 0.1
Rbase = 1
gdim = 2
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
model.add("model_acorn")
p = [
geom.addPoint(Rbase, np.pi, 0, ms),
geom.addPoint(Rbase, np.pi/2, 0, ms),
geom.addPoint(0, np.pi/2, 0, ms),
geom.addPoint(0, np.pi, 0, ms)
]
l = [
geom.addLine(p[0], p[1]),
geom.addLine(p[1], p[2]),
geom.addLine(p[2], p[3]),
geom.addLine(p[3], p[0])
]
cl = geom.addCurveLoop([l[j] for j in range(len(l))])
s = geom.addPlaneSurface([cl])
acorn = model.addPhysicalGroup(gdim, [s])
model.setPhysicalName(gdim, acorn, "Flat acorn")
base = model.addPhysicalGroup(gdim-1, [l[0]])
model.setPhysicalName(gdim-1, base, "Base")
cap = model.addPhysicalGroup(gdim-1, [l[1]])
model.setPhysicalName(gdim-1, cap, "Cap")
geom.synchronize()
model.mesh.generate(2)
domain, cell_markers, facet_markers = gmshio.model_to_mesh(model, mesh_comm, 0, gdim=gdim)
```