def METHOD_NAME(self, p1, p2, p3=None):
    ldof1 = self.number_of_dofs(p1)
    ldof2 = self.number_of_dofs(p2)
    phi1 = self.basis(p1)
    phi2 = self.basis(p2)
    M = sp.tensor.array.MutableDenseNDimArray(sp.zeros(ldof1*ldof2),(1,ldof1,ldof2))
    for i in range(ldof1):
        for j in range(ldof2):
            M[0,i, j] = self.integrate(phi1[i]*phi2[j])
    return M