def METHOD_NAME(self):
    p =self.p
    mesh = self.mesh
    space = self.space
    edge2cell = mesh.ds.edge_to_cell()
    qf = GaussLegendreQuadrature(p + 4)
    ws = qf.weights
    eh = mesh.entity_measure('edge')
    jump = self.Jump()
    gphi = self.Grad_average()
    n = mesh.edge_unit_normal()
    # Construct the flux matrix of computing element
    S = np.einsum('i, ijk, ijmp, jp, j->jkm', ws, jump, gphi, n, eh)
    # Construct the global flux matrix
    ldof = space.number_of_local_dofs(p=p, doftype='cell')
    cell2dof = space.cell_to_dof()
    dof1  = cell2dof[edge2cell[:,0]]
    dof2  = cell2dof[edge2cell[:,1]]
    dof = np.append(dof1,dof2,axis = 1)
    I = np.einsum('k, ij->ijk', np.ones(2*ldof), dof)
    J = I.swapaxes(-1, -2)
    gdof = space.number_of_global_dofs(p=p)
    S = csr_matrix((S.flat, (I.flat, J.flat)), shape=(gdof, gdof))
    return S