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.