def METHOD_NAME(self, p1, p2, p3):
    l = self.l
    GD = self.GD
    ldof1 = self.number_of_dofs(p1)
    ldof2 = self.number_of_dofs(p2)
    ldof3 = self.number_of_dofs(p3)
    phi1 = self.basis(p1)
    phi2 = self.basis(p2)
    phi3 = self.basis(p3)
    #S = np.zeros(shape = (ldof1, ldof2, ldof3))
    S = sp.tensor.array.MutableDenseNDimArray(sp.zeros(ldof1*ldof2*ldof3)\
            ,(ldof1, ldof2 ,ldof3))
    for i in range(ldof1):
        for j in range(ldof2):
            for k in range(ldof3):
                    temp= phi1[i]*phi2[j]*phi3[k]
                    S[i, j, k] = self.integrate(temp) 
    return S