Efficiently Finding Maximum Function Value in a Mesh with FEniCSx

What do you mean by this? If you use a P1 function space, the dofs are located at the vertices, whatever those vertices are are not affected by V.tabulate_dof_coordinates() which I just used to illustrate to you want the coordinates are.
The true value at the vertex is the value that is in u.x.array, u.eval will be affected by floating point accuracy in the push forward.

Secondly:

is not correct, as I mentioned earlier. See for instance: What is the difference between compute_incident_entities and entities_to_geometry? - #2 by dokken
and Getting coordinates of facets in dolfinx - #2 by dokken