def METHOD_NAME(self):
    mesh = self.mesh
    area = self.area
    gradphi = self.gradphi 
    A = self.A
    B = self.B
    NC = mesh.number_of_cells() 
    NN = mesh.number_of_nodes() 
    node = mesh.entity('node')
    edge = mesh.entity('edge')
    cell = mesh.entity('cell')
    edge2cell = mesh.ds.edge_to_cell()
    isBdEdge = (edge2cell[:,0]==edge2cell[:,1])
    bdEdge = edge[isBdEdge]
    cellIdx = edge2cell[isBdEdge, [0]]
    # construct the unit outward normal on the boundary
    W = np.array([[0, -1], [1, 0]], dtype=np.int)
    n = (node[bdEdge[:,1],] - node[bdEdge[:,0],:])@W
    h = np.sqrt(np.sum(n**2, axis=1)) 
    n /= h.reshape(-1, 1)
    # 计算梯度的恢复矩阵
    I = np.einsum('ij, k->ijk',  cell, np.ones(3))
    J = I.swapaxes(-1, -2)
    val = np.einsum('i, ij, ik->ijk', area, gradphi[:, :, 0], gradphi[:, :, 0])
    P = csc_matrix((val.flat, (I.flat, J.flat)), shape=(NN, NN))
    val = np.einsum('i, ij, ik->ijk', area, gradphi[:, :, 0], gradphi[:, :, 1])
    Q = csc_matrix((val.flat, (I.flat, J.flat)), shape=(NN, NN))
    val = np.einsum('i, ij, ik->ijk', area, gradphi[:, :, 1], gradphi[:, :, 1])
    S = csc_matrix((val.flat, (I.flat, J.flat)), shape=(NN, NN))
    K = A.transpose()@P@A + A.transpose()@Q@B + B.transpose()@Q.transpose()@A+B.transpose()@S@B 
    # 中间的边界上的两项
    I = np.einsum('ij, k->ijk', bdEdge, np.ones(3))
    J = np.einsum('ij, k->ikj', cell[cellIdx], np.ones(2)) 
    val0 = 0.5*h.reshape(-1, 1)*n[:, [0]]*gradphi[cellIdx, :, 0]  
    val0 = np.repeat(val0, 2, axis=0).reshape(-1, 2, 2)
    P0 = csc_matrix((val0.flat, (I.flat, J.flat)), shape=(NN, NN))
    val0 = 0.5*h.reshape(-1, 1)*n[:, [0]]*gradphi[cellIdx, :, 1]  
    val0 = np.repeat(val0, 2, axis=0).reshape(-1, 2, 2)
    Q0 = csc_matrix((val0.flat, (I.flat, J.flat)), shape=(NN, NN))
    val0 = 0.5*h.reshape(-1, 1)*n[:, [1]]*gradphi[cellIdx, :, 0]  
    val0 = np.repeat(val0, 2, axis=0).reshape(-1, 2, 2)
    P1 = csc_matrix((val0.flat, (I.flat, J.flat)), shape=(NN, NN))
    val0 = 0.5*h.reshape(-1, 1)*n[:, [1]]*gradphi[cellIdx, :, 1]  
    val0 = np.repeat(val0, 2, axis=0).reshape(-1, 2, 2)
    Q1 = csc_matrix((val0.flat, (I.flat, J.flat)), shape=(NN, NN))
    M = A.transpose()@P0@A + A.transpose()@Q0@B + B.transpose()@P1@A + B.transpose()@Q1@B
    K -= (M + M.transpose())
    K *= self.pde.epsilon**2
    # 边界上两个方向导数相乘的积分
    I = np.einsum('ij, k->ijk', bdEdge, np.ones(2))
    J = I.swapaxes(-1, -2)
    val = np.array([(1/3, 1/6), (1/6, 1/3)])
    val0 = np.einsum('i, jk->ijk', n[:, 0]*n[:, 0]/h, val)        
    P = csc_matrix((val0.flat, (I.flat, J.flat)), shape=(NN, NN))
    val0 = np.einsum('i, jk->ijk', n[:, 0]*n[:, 1]/h, val)
    Q = csc_matrix((val0.flat, (I.flat, J.flat)), shape=(NN, NN))
    val0 = np.einsum('i, jk->ijk', n[:, 1]*n[:, 1]/h, val)
    S = csc_matrix((val0.flat, (I.flat, J.flat)), shape=(NN, NN))
    K +=self.pde.epsilon**2*(A.transpose()@P@A + A.transpose()@Q@B + B.transpose()@Q@A + B.transpose()@S@B)
    return K 