def __call__(self, volume, problem=None, data=None):
        problem = get_default(problem, self.problem)

        coef = nm.zeros((self.dim, self.dim), dtype=self.dtype)

        term_mode = self.term_mode
        equations, variables = problem.create_evaluable(self.expression,
                                                        term_mode=term_mode)

        if isinstance(self.set_variables, list):
            for ir in range(self.dim):
                self.set_variables_default(variables, ir, None, 'row',
                                           self.set_variables, data)
                for ic in range(self.dim):
                    self.set_variables_default(variables, None, ic, 'col',
                                               self.set_variables, data)
                    val = eval_equations(equations, variables,
                                         term_mode=term_mode)
                    coef[ir,ic] = val
        else:
            for ir in range(self.dim):
                self.set_variables(variables, ir, None, 'row', **data)
                for ic in range(self.dim):
                    self.set_variables(variables, None, ic, 'col', **data)
                    val = eval_equations(equations, variables,
                                         term_mode=term_mode)
                    coef[ir,ic] = val

        coef /= self._get_volume(volume)

        return coef

    def check_custom_sensitivity(self, term_desc, idsg, delta, dp_var_data, vec_ap):
        pb = self.apb

        domain = pb.domain
        regions = domain.regions
        materials = pb.materials

        variables = self.ofg_equations.variables
        aux = self.dpb.create_evaluable(term_desc, try_equations=False, var_dict=variables, extra_args={"mode": 0})
        check0_equations, check0_variables = aux

        aux = self.dpb.create_evaluable(term_desc, try_equations=False, var_dict=variables, extra_args={"mode": 1})
        check1_equations, check1_variables = aux

        var_data = pb.equations.get_state_parts(vec_ap)
        var_data.update(dp_var_data)

        check0_equations.set_data(var_data, ignore_unknown=True)
        check1_equations.set_data(var_data, ignore_unknown=True)

        dim = self.sp_boxes.dim
        n_mesh_nod = domain.shape.n_nod

        a_grad = []
        d_grad = []

        coors0 = domain.mesh.coors

        for nu in self.generate_mesh_velocity((n_mesh_nod, dim), [idsg]):
            check1_variables["Nu"].data_from_any(nu.ravel())

            aux = eval_equations(check1_equations, check1_variables)
            a_grad.append(aux)

            coorsp = coors0 + delta * nu
            pb.set_mesh_coors(coorsp, update_state=True)
            valp = eval_equations(check0_equations, check0_variables)

            coorsm = coors0 - delta * nu
            pb.set_mesh_coors(coorsm, update_state=True)
            valm = eval_equations(check0_equations, check0_variables)

            d_grad.append(0.5 * (valp - valm) / delta)

        pb.set_mesh_coors(coors0, update_state=True)

        a_grad = nm.array(a_grad, nm.float64)
        d_grad = nm.array(d_grad, nm.float64)

        output(term_desc + ":")
        output("       a: %.8e" % a_grad)
        output("       d: %.8e" % d_grad)
        output("-> ratio:", a_grad / d_grad)
        pause()

    def __call__( self, volume, problem = None, data = None ):
        problem = get_default( problem, self.problem )

        dim, sym = problem.get_dim( get_sym = True )

        aux = self.get_filename( data, 0, 0 )
        ts = TimeStepper( *HDF5MeshIO( aux ).read_time_stepper() )

        coef = nm.zeros( (ts.n_step, sym), dtype = nm.float64 )

        equations, variables = problem.create_evaluable(self.expression)

        self.set_variables(variables, None, None, 'col', **data)

        for ii, (ir, ic) in enumerate( iter_sym( dim ) ):
            io = HDF5MeshIO( self.get_filename( data, ir, ic ) )
            for step, time in ts:
                self.set_variables(variables, io, step, 'row', **data)

                val = eval_equations(equations, variables)

                coef[step,ii] = val

        coef /= volume

        return coef

    def evaluate(self, expression, try_equations=True, auto_init=False,
                 preserve_caches=False, copy_materials=True, integrals=None,
                 ebcs=None, epbcs=None, lcbcs=None,
                 ts=None, functions=None,
                 mode='eval', dw_mode='vector', term_mode=None,
                 var_dict=None, strip_variables=True, ret_variables=False,
                 verbose=True, extra_args=None, **kwargs):
        """
        Evaluate an expression, convenience wrapper of
        :func:`ProblemDefinition.create_evaluable` and
        :func:`eval_equations() <sfepy.fem.evaluate.eval_equations>`.

        Parameters
        ----------
        dw_mode : 'vector' or 'matrix'
            The assembling mode for 'weak' evaluation mode.
        term_mode : str
            The term call mode - some terms support different call modes
            and depending on the call mode different values are
            returned.
        ret_variables : bool
            If True, return the variables that were created to evaluate
            the expression.
        other : arguments
            See docstrings of :func:`ProblemDefinition.create_evaluable`.

        Returns
        -------
        out : array
            The result of the evaluation.
        variables : Variables instance
            The variables that were created to evaluate
            the expression. Only provided if `ret_variables` is True.
        """
        aux = self.create_evaluable(expression,
                                    try_equations=try_equations,
                                    auto_init=auto_init,
                                    preserve_caches=preserve_caches,
                                    copy_materials=copy_materials,
                                    integrals=integrals,
                                    ebcs=ebcs, epbcs=epbcs, lcbcs=lcbcs,
                                    ts=ts, functions=functions,
                                    mode=mode, var_dict=var_dict,
                                    strip_variables=strip_variables,
                                    extra_args=extra_args,
                                    verbose=verbose, **kwargs)
        equations, variables = aux

        out = eval_equations(equations, variables,
                             preserve_caches=preserve_caches,
                             mode=mode, dw_mode=dw_mode, term_mode=term_mode)

        if ret_variables:
            out = (out, variables)

        return out

    def obj_fun(self, vec_dp):
        """
        Objective function evaluation for given direct problem state.
        """
        var_data = self.dpb.equations.get_state_parts(vec_dp)
        var_data = remap_dict(var_data, self.var_map)

        self.of_equations.set_data(var_data, ignore_unknown=True)

        val = eval_equations(self.of_equations, self.of_variables)

        return nm.squeeze(val)

    def sensitivity(self, dp_var_data, state_ap, select=None):
        """
        Sensitivity of objective function evaluation for given direct
        and adjoint problem states.
        """
        apb = self.apb

        var_data = state_ap.get_parts()
        var_data.update(dp_var_data)

        self.ofg_equations.set_data(var_data, ignore_unknown=True)

        dim = self.sp_boxes.dim
        n_mesh_nod = apb.domain.shape.n_nod

        if select is None:
            idsgs = nm.arange( self.dsg_vars.n_dsg, dtype = nm.int32 )
        else:
            idsgs = select

        sa = []

        pbar = MyBar('sensitivity:')
        pbar.init(len(idsgs))

        shape = (n_mesh_nod, dim)
        for ii, nu in enumerate(self.generate_mesh_velocity(shape, idsgs)):
            pbar.update(ii)
            self.ofg_variables['Nu'].set_data(nu.ravel())

            ## from sfepy.base.ioutils import write_vtk
            ## cc = nla.norm( vec_nu )
            ## nun = nu / cc
            ## out = {'v' : Struct( mode = 'vertex', data = nun,
            ##                      ap_name = 'nic', dof_types = (0,1,2) )}
            ## fd = open( 'anim/pert_%03d.pvtk' % (ii+1), 'w' )
            ## write_vtk( fd, domain.mesh, out )
            ## fd.close()
            ## print ii

            val = eval_equations(self.ofg_equations, self.ofg_variables,
                                 term_mode=1, preserve_caches=True)

            sa.append( val )

        vec_sa = nm.array( sa, nm.float64 )
        return vec_sa

def compute_eigenmomenta(em_equation, var_name, problem, eig_vectors,
                         transform=None, progress_bar=None):
    """
    Compute the eigenmomenta corresponding to given eigenvectors.
    """
    n_dof, n_eigs = eig_vectors.shape

    equations, variables = problem.create_evaluable(em_equation)
    var = variables[var_name]

    n_c = var.n_components
    eigenmomenta = nm.empty((n_eigs, n_c), dtype=nm.float64)

    if progress_bar is not None:
        progress_bar.init(n_eigs - 1)

    for ii in xrange(n_eigs):
        if progress_bar is not None:
            progress_bar.update(ii)

        else:
            if (ii % 100) == 0:
                output('%d of %d (%f%%)' % (ii, n_eigs,
                                            100. * ii / (n_eigs - 1)))

        if transform is None:
            vec_phi, is_zero = eig_vectors[:,ii], False

        else:
            vec_phi, is_zero = transform(eig_vectors[:,ii], (n_dof / n_c, n_c))

        if is_zero:
            eigenmomenta[ii, :] = 0.0

        else:
            var.data_from_any(vec_phi.copy())

            val = eval_equations(equations, variables)

            eigenmomenta[ii, :] = val

    return eigenmomenta

    def eval_equations(self, names=None, preserve_caches=False,
                   mode='eval', dw_mode='vector', term_mode=None,
                   verbose=True):
        """
        Evaluate (some of) the problem's equations, convenience wrapper of
        :func:`eval_equations() <sfepy.fem.evaluate.eval_equations>`.

        Parameters
        ----------
        names : str or sequence of str, optional
            Evaluate only equations of the given name(s).
        preserve_caches : bool
            If True, do not invalidate evaluate caches of variables.
        mode : one of 'eval', 'el_avg', 'qp', 'weak'
            The evaluation mode - 'weak' means the finite element
            assembling, 'qp' requests the values in quadrature points,
            'el_avg' element averages and 'eval' means integration over
            each term region.
        dw_mode : 'vector' or 'matrix'
            The assembling mode for 'weak' evaluation mode.
        term_mode : str
            The term call mode - some terms support different call modes
            and depending on the call mode different values are
            returned.
        verbose : bool
            If False, reduce verbosity.

        Returns
        -------
        out : dict or result
            The evaluation result. In 'weak' mode it is the vector or sparse
            matrix, depending on `dw_mode`. Otherwise, it is a dict of results
            with equation names as keys or a single result for a single
            equation.
        """
        return eval_equations(self.equations, self.equations.variables,
                              names=names, preserve_caches=preserve_caches,
                              mode=mode, dw_mode=dw_mode, term_mode=term_mode,
                              verbose=verbose)

def compute_eigenmomenta( em_equation, u_name, pb, eig_vectors,
                          transform = None, pbar = None ):
    dim = pb.domain.mesh.dim
    n_dof, n_eigs = eig_vectors.shape

    eigenmomenta = nm.empty( (n_eigs, dim), dtype = nm.float64 )

    equations, variables = pb.create_evaluable(em_equation)

    var = variables[u_name]

    if pbar is not None:
        pbar.init( n_eigs - 1 )

    for ii in xrange( n_eigs ):
        if pbar is not None:
            pbar.update( ii )
        else:
            if (ii % 100) == 0:
                output( '%d of %d (%f%%)' % (ii, n_eigs,
                                             100. * ii / (n_eigs - 1)) )
            
        if transform is None:
            vec_phi, is_zero = eig_vectors[:,ii], False
        else:
            vec_phi, is_zero = transform(eig_vectors[:,ii], (n_dof / dim, dim))

        if is_zero:
            eigenmomenta[ii,:] = 0.0

        else:
            var.data_from_any(vec_phi.copy())

            val = eval_equations(equations, variables)

            eigenmomenta[ii,:] = val

    return eigenmomenta

    def __call__(self, volume, problem=None, data=None):
        problem = get_default(problem, self.problem)

        coef = nm.zeros((self.dim,), dtype=self.dtype)
        term_mode = self.term_mode

        for ir in range(self.dim):
            expression = self.expression % self.expr_pars[ir]
            equations, variables = \
              problem.create_evaluable(expression, term_mode=term_mode)

            if isinstance(self.set_variables, list):
                self.set_variables_default(variables, ir, self.set_variables,
                                           data)
            else:
                self.set_variables(variables, ir, **data)

            val = eval_equations(equations, variables,
                                 term_mode=term_mode)
            coef[ir] = val

        coef /= self._get_volume(volume)

        return coef

    def iterate(self, v_hxc_qp, eig_solver,
                mtx_a_equations, mtx_a_variables, mtx_b, log, file_output,
                n_electron=None):
        from sfepy.physics import dft

        self.itercount += 1

        pb = self.problem
        opts = self.app_options

        n_electron = get_default( n_electron, opts.n_electron )

        sh = self.qp_shape

        v_hxc_qp = nm.array(v_hxc_qp, dtype=nm.float64)
        v_hxc_qp.shape = (sh[0] * sh[1],) + sh[2:]

        mat_v = Materials(mtx_a_equations.collect_materials())['mat_v']
        mat_v.set_extra_args(vhxc=v_hxc_qp)
        mat_v.time_update(None, pb.domain, mtx_a_equations)

        v_hxc_qp.shape = sh

        v_ion_qp = mat_v.get_data(('Omega', 'i1'), 0, 'V_ion')

        output( 'assembling lhs...' )
        tt = time.clock()
        mtx_a = eval_equations(mtx_a_equations, mtx_a_variables,
                               mode='weak', dw_mode='matrix')
        output( '...done in %.2f s' % (time.clock() - tt) )

        assert_( nm.alltrue( nm.isfinite( mtx_a.data ) ) )

        output( 'computing the Ax=Blx Kohn-Sham problem...' )
        tt = time.clock()
        eigs, mtx_s_phi = eig_solver( mtx_a, mtx_b,
                                      opts.n_eigs, eigenvectors = True )
        output( '...done in %.2f s' % (time.clock() - tt) )
        n_eigs_ok = len(eigs)

        output( 'setting-up smearing...' )
        e_f, smear_tuned = setup_smearing( eigs, n_electron )
        output( 'Fermi energy:', e_f )
        if e_f is None:
            raise Exception("cannot find Fermi energy - exiting.")
        weights = smear_tuned(eigs)
        output( '...done' )

        if (weights[-1] > 1e-12):
            output("last smearing weight is nonzero (%s eigs ok)!" % n_eigs_ok)

        output( "saving solutions, iter=%d..." % self.itercount )
        out = {}
        var_name = mtx_a_variables.get_names(kind='state')[0]
        for ii in xrange( n_eigs_ok ):
            vec_phi = mtx_a_variables.make_full_vec(mtx_s_phi[:,ii])
            update_state_to_output( out, pb, vec_phi, var_name+'%03d' % ii )
        name = op.join( opts.output_dir, "iter%d" % self.itercount )
        pb.save_state('.'.join((name, opts.output_format)), out=out)
        output( "...solutions saved" )

        output('computing total charge...')
        tt = time.clock()
        aux = pb.create_evaluable('dq_state_in_volume_qp.i1.Omega(Psi)')
        psi_equations, psi_variables = aux
        var = psi_variables['Psi']

        n_qp = nm.zeros_like(v_hxc_qp)
        for ii in xrange( n_eigs_ok ):
            vec_phi = mtx_a_variables.make_full_vec(mtx_s_phi[:,ii])
            var.data_from_any(vec_phi)

            phi_qp = eval_equations(psi_equations, psi_variables)
            n_qp += weights[ii] * (phi_qp ** 2)
        output('...done in %.2f s' % (time.clock() - tt))

        ap, vg = var.get_approximation(('i1', 'Omega', 0), 'Volume')

        det = vg.variable(1)
        charge = (det * n_qp).sum()
        ## Same as above.
        ## out = nm.zeros((n_qp.shape[0], 1, 1, 1), dtype=nm.float64)
        ## vg.integrate(out, n_qp)
        ## charge = out.sum()

        vec_n = self._interp_to_nodes(n_qp)

        var.data_from_any(vec_n)
        charge_n = pb.evaluate('di_volume_integrate.i1.Omega(Psi)', Psi=var)

        ##
        # V_xc in quadrature points.
        v_xc_qp = nm.zeros((nm.prod(self.qp_shape),), dtype=nm.float64)
        for ii, val in enumerate(n_qp.flat):
            ## print ii, val
            v_xc_qp[ii] = dft.getvxc(val, 0)
        assert_(nm.isfinite(v_xc_qp).all())
        v_xc_qp.shape = self.qp_shape

        mat_key = mat_v.datas.keys()[0]
#.........这里部分代码省略.........

    def check_custom_sensitivity(self, term_desc, idsg, delta,
                                 dp_var_data, state_ap):
        pb = self.apb

        domain = pb.domain

        possible_mat_names = get_expression_arg_names(term_desc)
        materials = self.dpb.create_materials(possible_mat_names).as_dict()

        variables = self.ofg_equations.variables
        aux = self.dpb.create_evaluable(term_desc,
                                        try_equations=False,
                                        var_dict=variables,
                                        verbose=False,
                                        **materials)
        check0_equations, check0_variables = aux

        aux = self.dpb.create_evaluable(term_desc,
                                        try_equations=False,
                                        var_dict=variables,
                                        verbose=False,
                                        **materials)
        check1_equations, check1_variables = aux

        var_data = state_ap.get_parts()
        var_data.update(dp_var_data)

        check0_equations.set_data(var_data, ignore_unknown=True)
        check1_equations.set_data(var_data, ignore_unknown=True)

        dim = self.sp_boxes.dim
        n_mesh_nod = domain.shape.n_nod

        a_grad = []
        d_grad = []

        coors0 = domain.mesh.coors

        for nu in self.generate_mesh_velocity( (n_mesh_nod, dim), [idsg] ):
            check1_variables['Nu'].set_data(nu.ravel())

            aux = eval_equations(check1_equations, check1_variables,
                                 term_mode=1)
            a_grad.append( aux )

            coorsp = coors0 + delta * nu
            pb.set_mesh_coors( coorsp, update_fields=True )
            valp = eval_equations(check0_equations, check0_variables,
                                  term_mode=0)

            coorsm = coors0 - delta * nu
            pb.set_mesh_coors( coorsm, update_fields=True )
            valm = eval_equations(check0_equations, check0_variables,
                                  term_mode=0)

            d_grad.append( 0.5 * (valp - valm) / delta )

        pb.set_mesh_coors( coors0, update_fields=True )

        a_grad = nm.array( a_grad, nm.float64 )
        d_grad = nm.array( d_grad, nm.float64 )

        output( term_desc + ':' )
        output( '       a: %.8e' % a_grad )
        output( '       d: %.8e' % d_grad )
        output( '-> ratio:', a_grad / d_grad )
        pause()