Interpolate a function given as matrix on a finite element basis

Question

Hi all,

Let phi be a function given as a matrix N*N and representing the values of the function at the vertices of a unit square mesh. I want to interpolate this function on the finite element basis. The interpolated function is called psi. I tried this:

mesh = RectangleMesh(0, 0, 1, 1, N, N)
V   = FunctionSpace(mesh, "Lagrange", 1)
psi = Function(V)
phi_vector  = phi.reshape((N+1)*(N+1))
psi.vector()[:] = phi_vector

but the resulting function psi is not correct. I guess this is because the basis functions are not indexed like the mesh vertices. Any idea how to do the interpolation then? Is there a way to find the indexing of the basis functions with respect to the mesh vertices?

thanks
Antoine

Answer 1

This is the evergreen of the FEniCS questions. There are basically these approaches:

• enforce ordering of DOFs according to underlying mesh entities (here vertices) by

  parameters['reorder_dofs_serial'] = False
  

  but this will only work in serial.

• use helper member functions of V.dofmap(), in this case

  psi.vector()[:] = phi_vector[V.dofmap().vertex_to_dof_map(mesh)]
  

  should do the job if phi_vector is indexed by local vertex indices.

• interpolate or project subclassed Expression. This applies well to DG0 spaces as Expression.eval interface supplies cell index to users evaluation algorithm. Hence it is not probably useful for your case.

For further reading refer to

• http://fenicsproject.org/qa/121/python-how-to-set-expansion-coefficients-for-a-function
• http://fenicsproject.org/qa/78/assigning-expansion-coefficients-function-does-parallel
• http://fenicsproject.org/qa/207/how-to-interpolate-data-at-vertices-of-3d-cells

Comments

Following your advice I have tried this:

psi.vector()[:] = phi_vector[V.dofmap().vertex_to_dof_map(mesh)]

and it works, many thanks.