mesh.pl.get_numbering(SF::NBR_REF) are again global indices. You are comparing global with local indices.

Lets say we are looking at the local vertex index 5. Then mesh.xyz[5*3+0] till mesh.xyz[5*3+2] are the vertex coordinates (of course usually we dont access them this way, but via an element_view), mesh.pl.get_numbering(SF::NBR_REF)[5] is the (global) index in the mesh on the hard disk of this vertex, and mesh.pl.get_numbering(SF::NBR_PETSC)[5] is the (global) index in the parallel, distributed linear equation system.

The whole point of storing the numberings is to not have to communicate.

Coming back to algebraic_nodes(): Lets denote Np := mesh.l_numpts on rank p. Then the implicit local index range is [0, Np-1] (this local range is not stored since trivial). algebraic_nodes() is some subset of this interval.

For example, if Np := 5, then the local range is [0 1 2 3 4] and one possible algebraic_nodes() set would be [1 3].. Hope that makes sense now.

The reason we have algebraic_nodes() is that in many algorithms one needs to work on a non-overlapping parallel distribution of the nodes, while in other algorithms one needs a non-overlapping distribution of the *elements* (thus a overlapping node distribution). By storing algebraic_nodes() we can access both overlapping and non-overlapping local nodal index ranges.