def METHOD_NAME():
    parent_mesh = UnitSquareMesh(10, 10)
    V = FunctionSpace(parent_mesh, "CG", 2)
    # Create a function f on the parent mesh to point evaluate
    x, y = SpatialCoordinate(parent_mesh)
    f = Function(V).interpolate(x**2 + y**2)
    # 3 points (i.e. vertices) at which to point evaluate f
    points = [[0.1, 0.1], [0.2, 0.2], [0.3, 0.3]]
    vom = VertexOnlyMesh(parent_mesh, points)
    # P0DG is the only function space you can make on a vertex-only mesh
    P0DG = FunctionSpace(vom, "DG", 0)
    # Interpolation performs point evaluation
    f_at_points = interpolate(f, P0DG)
    print(f_at_points.dat.data_ro)
    # Make a P0DG function on the input ordering vertex-only mesh again
    P0DG_input_ordering = FunctionSpace(vom.input_ordering, "DG", 0)
    f_at_input_points = Function(P0DG_input_ordering)
    # We interpolate the other way this time
    f_at_input_points.interpolate(f_at_points)
    print(f_at_input_points.dat.data_ro)  # will print the values at the input points
    import numpy as np
    f_at_input_points.dat.data_wo[:] = np.nan
    f_at_input_points.interpolate(f_at_points)
    print(f_at_input_points.dat.data_ro)  # any points not found will be NaN