Could you please verify if I understand your question as the following;

Let us assume that the function space has N degrees of freedom and i,j,k\in [0, N-1].

A function u in this function space is defined as u(x)= \sum_{l=0}^n c_l \phi_l(x) where \phi_l is the $l$th basis function.

You would like to assemble \int_{\Omega} \phi_i \phi_j\phi_k ~\mathrm{d} x?

Then you could do something like the following:

```
from dolfin import *
mesh = UnitSquareMesh(3, 3)
V = FunctionSpace(mesh, "CG", 1)
print(V.tabulate_dof_coordinates())
i, j, k = 0, 1, 2
ui = Function(V)
ui.vector()[i] = 1
uj = Function(V)
uj.vector()[j] = 1
uk = Function(V)
uk.vector()[k] = 1
print(assemble(ui*uj*uk*dx))
```