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 compute_eigenmomenta(em_equation, var_name, problem, eig_vectors,
                         transform=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)

    for ii in xrange(n_eigs):

        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.set_data(vec_phi.copy())

            val = eval_equations(equations, variables)

            eigenmomenta[ii, :] = val

    return eigenmomenta

    def get_coef(self, row, col, volume, problem, data):
        problem = get_default(problem, self.problem)
        term_mode = self.term_mode
        equations, variables = problem.create_evaluable(self.expression,
                                                        term_mode=term_mode)

        coef = nm.zeros((len(row), len(col)), dtype=self.dtype)

        for ir, (irr, icr) in enumerate(row):
            if isinstance(self.set_variables, list):
                self.set_variables_default(variables, irr, icr, 'row',
                                           self.set_variables, data)
            else:
                self.set_variables(variables, irr, icr, 'row', **data)

            for ic, (irc, icc) in enumerate(col):
                if isinstance(self.set_variables, list):
                    self.set_variables_default(variables, irc, icc, 'col',
                                               self.set_variables, data)
                else:
                    self.set_variables(variables, irc, icc, 'col', **data)

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

        coef /= self._get_volume(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:`Problem.create_evaluable` and
        :func:`eval_equations() <sfepy.discrete.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:`Problem.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, state_dp):
        """
        Objective function evaluation for given direct problem state.
        """
        var_data = state_dp.get_parts()
        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 = []

        output('computing sensitivity of %d variables...' % idsgs)

        shape = (n_mesh_nod, dim)
        for ii, nu in enumerate(self.generate_mesh_velocity(shape, idsgs)):
            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 )
        output('...done')

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

    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()