u = TrialFunction(V)
v = TestFunction(V)
M = assemble(u*v*dx)

i wonder if one can compute integral of v_iv_jv_k over the domain, where v_i,v_j and v_k are the testfunctions. I thought about exporting the testfunctions from somewhere, but i dont know where and how …

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))