Evaluation on curvilinear boundary

I have a domain with a curved boundary (in fact, it is just a disc, for now). I would like to evaluate my function (pointwise), along the boundary, but I am struggling with the issue that, in fenics, the boundary is actually polygonal. Consequently, simply evaluating at
(r\cos(\theta), r\sin(\theta)) can generate errors if it’s outside the domain. Is there a straightforward way to do this? It would even be sufficient if I could just extract the nodal values.

The easiest approach would be to set the functions to allow extrapolation:

from dolfin import *
from mshr import *
import math

mesh = generate_mesh(Circle(Point(0,0),1),100)
x = SpatialCoordinate(mesh)
f = project(x[0]**2,FunctionSpace(mesh,"CG",1))
f.set_allow_extrapolation(True)

N = 10
thetaInc = 2.0*math.pi/N
for i in range(0,N):
    thetai = i*thetaInc
    xi = math.cos(thetai)
    yi = math.sin(thetai)
    print(thetai , xi**2, f(xi,yi))

One follow up question: Is there a straightforward way to either extract function elements and/or evaluate the function at the boundary nodes?

You can dig some mappings up from the API, e.g.,

from dolfin import *
from mshr import *
mesh = generate_mesh(Circle(Point(0,0),1),10)
bmesh = BoundaryMesh(mesh,"exterior")
V = FunctionSpace(mesh,"CG",1)
v2d = vertex_to_dof_map(V) # Map vertices to DoFs
entity_dim = 0 # For vertices
b2v = bmesh.entity_map(entity_dim).array() # Map bmesh verts to mesh verts
bdofs = v2d[b2v] # Compose mappings for indices of boundary DoFs in V

# Test: 
x = SpatialCoordinate(mesh)
r = project(sqrt(dot(x,x)),V)
# Should all be ~1 if mappings are correct:
print(r.vector()[bdofs])

although I’m not sure what the interpretation of vertex_to_dof_map() is for vector or mixed function spaces.