Note that V.tabulate_dof_coordinates() gives you all dof coordinates on your process. You can in turn use V.dofmap().cell_dofs(i)` to Get the indices for the ith cell.
You have already asked questions regarding this at: How to get Gaussian integration point directly on a mesh element/cell? - #2 by dokken
In legacy dolfin this is not straightforward, as you would have to evaluate the basis functions and then combine them to be the function at that point. See Properly calling evaluate basis derivatives all() function in python - #3 by Ryan_Vogt2
In DOLFINx, you have a larger variety of quadrature rules, as well as the possibility to evaluate function derivatives at any point in the cell (including interpolation points). See for instance Solving a time-dependent problem — FEniCS 22 tutorial