def METHOD_NAME(self, source=None, mu=None, kappa=1, eps=None):
    p = 2
    lspace = LagrangeFiniteElementSpace(self.mesh, p=p, spacetype='D')
    luh = lspace.function()
    qf = self.mesh.integrator(q=8)
    bcs, ws = qf.get_quadrature_points_and_weights()
    point = self.mesh.bc_to_point(bcs)
    cellmeasure = self.mesh.entity_measure('cell')
    NC = self.mesh.number_of_cells()
    NQ = len(bcs)
    val = np.zeros(shape=(NQ, NC, 2))  # 被积函数
    '''
    1. 估计单元内的后验误差
    R₁ = ▽ × (1/μ ▽ × Eₕ) - k²𝜀Eₕ - F
    1.1 计算 1/μ ▽ × Eₕ , 并将其插值到 Lagrange 空间
    '''
    # 0.0 构造一个 L2 的函数
    ib = self.mesh.multi_index_matrix(p) / p
    ps = self.mesh.bc_to_point(ib)
    # 1/μ ▽ × Eₕ
    if mu == None:
        luh[:] = self.uh.curl_value(ib).flatten()
    else:
        luh[:] = self.uh.curl_value(ib).flatten() / mu(ps).flatten()
    grad_luh = lspace.grad_value(luh, bcs)
    p_x = grad_luh[..., 0]  # (NQ, NC)
    p_y = grad_luh[..., 1]
    '''
    若 u 是一个标量, 那么 ▽ × u = (∂ᵧu, -∂ₓu)ᵀ
    '''
    # cuh 存储了 (∂ᵧu, -∂ₓu)ᵀ
    val[..., 0] = p_y
    val[..., 1] = -p_x
    k2 = kappa ** 2
    eh = self.space.value(self.uh, bcs)
    if eps == None:
        val[:] -= k2 * eh
    else:
        eps = eps(point)
        val[:] -= k2 * np.einsum('qcij, qci -> qcj', eps, eh)
    if source == None:
        val = 0 - val
    else:
        f = source(point)
        val = f - val
    '''
    下面是对 R1 在每个单元上做积分
    '''
    val = np.power(val, 2)
    intR1 = np.einsum('i, ija, j -> j', ws, val, cellmeasure)
    return intR1