def METHOD_NAME(self, p1, p2):
    l = self.l
    GD = self.GD
    ldof1 = self.number_of_dofs(p1)
    ldof2 = self.number_of_dofs(p2)
    phi1 = self.basis(p1)
    phi2 = self.basis(p2)
    S =sp.tensor.array.MutableDenseNDimArray(sp.zeros(ldof1*ldof2*(GD+1))*(GD+1)\
            , (ldof1,ldof2,GD+1,GD+1))
    for i in range(ldof1):
        for j in range(ldof2):
            for m in range(GD + 1):
                for n in range(GD + 1):
                    temp= sp.diff(phi1[i],l[m])*sp.diff(phi2[j],l[n])
                    S[i,j,m,n] = self.integrate(temp) 
    return S