how to evaluate higher derivative of a Function at a point

Question

Suppose u is a function in a FunctionSpace, say CG3, and I need to evaluate a partial derivative of u at a point, e.g., evaluate d^2 u/dx dy at (.5, .5). What is the correct way to do this?

I could project u.dx(0).dx(1) into DG1 and then evaluate, but this seems inefficient.

Answer 1

Hi, here's my attempt. I hope there will be more suggestions.

from dolfin import *
import numpy as np

# Have the compiler generate code for evaluating derivatives
parameters['form_compiler']['no-evaluate_basis_derivatives'] = False

mesh = UnitSquareMesh(3, 3)
V = FunctionSpace(mesh, 'CG', 3)
f = interpolate(Expression('3*x[0]*x[0]*x[1]-2*x[1]*x[1]*x[0]'), V)
el = V.element()

# Where to evaluate
x = np.array([0.33, 0.55])

# Find the cell with point
x_point = Point(*x) 
cell_id = mesh.bounding_box_tree().compute_first_entity_collision(x_point)
cell = Cell(mesh, cell_id)
coordinate_dofs = cell.get_vertex_coordinates()

# Array for values with derivatives of all basis functions. 4 * element dim
values = np.zeros(4*el.space_dimension(), dtype=float)
# Compute all 2nd order derivatives
el.evaluate_basis_derivatives_all(2, values, x, coordinate_dofs, cell.orientation())
# Reshape such that colums are [d/dxx, d/dxy, d/dyx, d/dyy]
values = values.reshape((-1, 4))

# Get expansion coefs on cell. Alternative to this is f.restrict(...)
dofs = V.dofmap().cell_dofs(cell_id)
dofs = f.vector()[dofs]

# Combine
print np.sum(values[:, 2]*dofs)