def save_only( conf, save_names, problem = None ):
    """Save information available prior to setting equations and
    solving them."""
    if problem is None:
        problem = ProblemDefinition.from_conf( conf, init_variables = False )

    if save_names.regions is not None:
        problem.save_regions( save_names.regions )

    if save_names.field_meshes is not None:
        problem.save_field_meshes( save_names.field_meshes )

    if save_names.region_field_meshes is not None:
        problem.save_region_field_meshes( save_names.region_field_meshes )

    if save_names.ebc is not None:
        if not hasattr( problem, 'variables' ):
            problem.set_variables( conf.variables )
        try:
            ts = TimeStepper.from_conf( conf.ts )
            ts.set_step( 0 )
        except:
            ts = None
        try:
            problem.variables.equation_mapping( conf.ebcs, conf.epbcs,
                                                problem.domain.regions, ts,
                                                conf.funmod )
        except Exception, e:
            output( 'cannot make equation mapping!' )
            output( 'reason: %s' % e )
        else:
            problem.save_ebc( save_names.ebc )

def vary_omega1_size( problem ):
    """Vary size of \Omega1. Saves also the regions into options['output_dir'].

    Input:
      problem: ProblemDefinition instance
    Return:
      a generator object:
      1. creates new (modified) problem
      2. yields the new (modified) problem and output container
      3. use the output container for some logging
      4. yields None (to signal next iteration to Application)
    """
    from sfepy.fem import ProblemDefinition
    from sfepy.solvers.ts import get_print_info
    
    output.prefix = 'vary_omega1_size:'

    diameters = nm.linspace( 0.1, 0.6, 7 ) + 0.001
    ofn_trunk, output_format = problem.ofn_trunk, problem.output_format
    output_dir = problem.output_dir
    join = os.path.join

    conf = problem.conf
    cf = conf.get_raw( 'functions' )
    n_digit, aux, d_format = get_print_info( len( diameters ) + 1 )
    for ii, diameter in enumerate( diameters ):
        output( 'iteration %d: diameter %3.2f' % (ii, diameter) )

        cf['select_circ'] = (lambda coors, domain=None: 
                             select_circ(coors[:,0], coors[:,1], 0, diameter),)
        conf.edit('functions', cf)
        problem = ProblemDefinition.from_conf( conf )

        problem.save_regions( join( output_dir, ('regions_' + d_format) % ii ),
                              ['Omega_1'] )
        region = problem.domain.regions['Omega_1']
        if not region.has_cells_if_can():
            print region
            raise ValueError( 'region %s has no cells!' % region.name )

        ofn_trunk = ofn_trunk + '_' + (d_format % ii)
        problem.setup_output(output_filename_trunk=ofn_trunk,
                             output_dir=output_dir,
                             output_format=output_format)

        out = []
        yield problem, out

        out_problem, state = out[-1]

        filename = join( output_dir,
                         ('log_%s.txt' % d_format) % ii )
        fd = open( filename, 'w' )
        log_item = '$r(\Omega_1)$: %f\n' % diameter
        fd.write( log_item )
        fd.write( 'solution:\n' )
        nm.savetxt(fd, state())
        fd.close()

        yield None

def main():
    from sfepy.base.base import output
    from sfepy.base.conf import ProblemConf, get_standard_keywords
    from sfepy.fem import ProblemDefinition
    from sfepy.applications import solve_evolutionary

    output.prefix = 'therel:'

    required, other = get_standard_keywords()
    conf = ProblemConf.from_file(__file__, required, other)

    problem = ProblemDefinition.from_conf(conf, init_equations=False)

    # Setup output directory according to options above.
    problem.setup_default_output()

    # First solve the stationary electric conduction problem.
    problem.set_equations({'eq' : conf.equations['1']})
    problem.time_update()
    state_el = problem.solve()
    problem.save_state(problem.get_output_name(suffix = 'el'), state_el)

    # Then solve the evolutionary heat conduction problem, using state_el.
    problem.set_equations({'eq' : conf.equations['2']})
    phi_var = problem.get_variables()['phi_known']
    phi_var.data_from_any(state_el())
    solve_evolutionary(problem)

    output('results saved in %s' % problem.get_output_name(suffix = '*'))

    def from_conf( conf, options ):
        from sfepy.fem import ProblemDefinition

        problem = ProblemDefinition.from_conf(conf, init_equations=False)
        test = Test( problem = problem,
                     conf = conf, options = options )
        return test

def create_problem(filename):
    from sfepy.fem import ProblemDefinition

    problem = ProblemDefinition.from_conf_file(filename,
                                               init_equations=False,
                                               init_solvers=False)
    return problem

def main():
    from sfepy.base.conf import ProblemConf, get_standard_keywords
    from sfepy.fem import ProblemDefinition
    from sfepy.base.plotutils import plt

    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-n', '--no-plot',
                      action="store_true", dest='no_plot',
                      default=False, help=helps['no_plot'])
    options, args = parser.parse_args()

    required, other = get_standard_keywords()
    # Use this file as the input file.
    conf = ProblemConf.from_file( __file__, required, other )

    # Create problem instance, but do not set equations.
    problem = ProblemDefinition.from_conf( conf,
                                           init_equations = False )

    # Solve the problem. Output is ignored, results stored by using the
    # step_hook.
    u_t = solve_branch(problem, linear_tension)
    u_c = solve_branch(problem, linear_compression)

    # Get pressure load by calling linear_*() for each time step.
    ts = problem.get_timestepper()
    load_t = nm.array([linear_tension(ts, nm.array([[0.0]]), 'qp')['val']
                       for aux in ts.iter_from( 0 )],
                      dtype=nm.float64).squeeze()
    load_c = nm.array([linear_compression(ts, nm.array([[0.0]]), 'qp')['val']
                       for aux in ts.iter_from( 0 )],
                      dtype=nm.float64).squeeze()

    # Join the branches.
    displacements = {}
    for key in u_t.keys():
        displacements[key] = nm.r_[u_c[key][::-1], u_t[key]]
    load = nm.r_[load_c[::-1], load_t]


    if plt is None:
        print 'matplotlib cannot be imported, printing raw data!'
        print displacements
        print load
    else:
        legend = []
        for key, val in displacements.iteritems():
            plt.plot( load, val )
            legend.append( key )

        plt.legend( legend, loc = 2 )
        plt.xlabel( 'tension [kPa]' )
        plt.ylabel( 'displacement [mm]' )
        plt.grid( True )

        plt.gcf().savefig( 'pressure_displacement.png' )

        if not options.no_plot:
            plt.show()

    def from_conf(conf, options):
        from sfepy.fem import ProblemDefinition

        problem = ProblemDefinition.from_conf(conf)
        problem.time_update()

        test = Test(problem=problem, conf=conf, options=options)
        return test

def solve_direct(conf, options, problem=None, step_hook=None,
                 post_process_hook=None, post_process_hook_final=None,
                 pre_process_hook=None, nls_status=None):
    """Generic (simple) problem solver."""

    if problem is None:
        is_eqs = not options.solve_not
        problem = ProblemDefinition.from_conf(conf, init_equations=is_eqs)

        problem.setup_default_output(conf, options)

    if pre_process_hook is not None: # User pre_processing.
        pre_process_hook(problem)

    ofn_trunk = problem.ofn_trunk

    save_names = Struct( ebc = None, regions = None,
                         regions_as_groups = None, field_meshes = None,
                         region_field_meshes = None )
    if options.save_ebc:
        save_names.ebc = ofn_trunk + '_ebc.vtk'
    if options.save_regions:
        save_names.regions = ofn_trunk + '_region'
    if options.save_regions_as_groups:
        save_names.regions_as_groups = ofn_trunk + '_regions'
    if options.save_field_meshes:
        save_names.field_meshes = ofn_trunk + '_field'

    is_extra_save = False
    for name, val in save_names.to_dict().iteritems():
        if val is not None:
            is_extra_save = True
            break
    if is_extra_save:
        save_only( conf, save_names, problem=problem )

    if options.solve_not:
        return None, None, None

    if hasattr( conf.options, 'ts' ):
        ##
        # Time-dependent problem.
        state = solve_evolutionary_op(problem, options,
                                      step_hook=step_hook,
                                      post_process_hook=post_process_hook,
                                      nls_status=nls_status)
    else:
        ##
        # Stationary problem.
        state = solve_stationary_op(problem, options,
                                    post_process_hook=post_process_hook,
                                    nls_status=nls_status)

    if post_process_hook_final is not None: # User postprocessing.
       post_process_hook_final(problem, state)

    return problem, state

def main():
    from sfepy import data_dir

    parser = OptionParser(usage=usage, version="%prog")
    parser.add_option("-s", "--show", action="store_true", dest="show", default=False, help=help["show"])
    options, args = parser.parse_args()

    mesh = Mesh.from_file(data_dir + "/meshes/2d/rectangle_tri.mesh")
    domain = Domain("domain", mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
    eps = 1e-8 * (max_x - min_x)
    omega = domain.create_region("Omega", "all")
    gamma1 = domain.create_region("Gamma1", "nodes in x < %.10f" % (min_x + eps))
    gamma2 = domain.create_region("Gamma2", "nodes in x > %.10f" % (max_x - eps))

    field = Field("fu", nm.float64, "vector", omega, space="H1", poly_space_base="lagrange", approx_order=2)

    u = FieldVariable("u", "unknown", field, mesh.dim)
    v = FieldVariable("v", "test", field, mesh.dim, primary_var_name="u")

    m = Material("m", lam=1.0, mu=1.0)
    f = Material("f", val=[[0.02], [0.01]])

    integral = Integral("i", order=3)

    t1 = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v, u)", integral, omega, m=m, v=v, u=u)
    t2 = Term.new("dw_volume_lvf(f.val, v)", integral, omega, f=f, v=v)
    eq = Equation("balance", t1 + t2)
    eqs = Equations([eq])

    fix_u = EssentialBC("fix_u", gamma1, {"u.all": 0.0})

    bc_fun = Function("shift_u_fun", shift_u_fun, extra_args={"shift": 0.01})
    shift_u = EssentialBC("shift_u", gamma2, {"u.0": bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = ProblemDefinition("elasticity", equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups("regions")

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print nls_status

    pb.save_state("linear_elasticity.vtk", vec)

    if options.show:
        view = Viewer("linear_elasticity.vtk")
        view(vector_mode="warp_norm", rel_scaling=2, is_scalar_bar=True, is_wireframe=True)

def make_h1_projection_data(target, eval_data):
    """
    Project scalar data given by a material-like `eval_data()` function to a
    scalar `target` field variable using the :math:`H^1` dot product.
    """
    order = target.field.approx_order * 2
    integral = Integral('i', order=order)

    un = target.name
    v = FieldVariable('v', 'test', target.field, 1, primary_var_name=un)
    lhs1 = Term.new('dw_volume_dot(v, %s)' % un, integral,
                    target.field.region, v=v, **{un : target})
    lhs2 = Term.new('dw_laplace(v, %s)' % un, integral,
                    target.field.region, v=v, **{un : target})

    def _eval_data(ts, coors, mode, **kwargs):
        if mode == 'qp':
            val = eval_data(ts, coors, mode, 'val', **kwargs)
            gval = eval_data(ts, coors, mode, 'grad', **kwargs)
            return {'val' : val, 'gval' : gval}

    m = Material('m', function=_eval_data)
    rhs1 = Term.new('dw_volume_lvf(m.val, v)', integral, target.field.region,
                    m=m, v=v)
    rhs2 = Term.new('dw_diffusion_r(m.gval, v)', integral, target.field.region,
                    m=m, v=v)

    eq = Equation('projection', lhs1 + lhs2 - rhs1 - rhs2)
    eqs = Equations([eq])

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = ProblemDefinition('aux', equations=eqs, nls=nls, ls=ls)

    pb.time_update()

    # This sets the target variable with the projection solution.
    pb.solve()

    if nls_status.condition != 0:
        output('H1 projection: solver did not converge!')

def make_l2_projection(target, source):
    """
    Project `source` field variable to `target` field variable using
    :math:`L^2` dot product.
    """
    order = target.field.get_true_order()**2
    integral = Integral('i', order=order)

    un = target.name
    v = FieldVariable('v', 'test', target.field, 1, primary_var_name=un)
    lhs = Term.new('dw_mass_scalar(v, %s)' % un, integral,
                   target.field.region, v=v, **{un : target})

    def eval_variable(ts, coors, mode, **kwargs):
        if mode == 'qp':
            val = source.evaluate_at(coors)
            val.shape = val.shape + (1,)
            out = {'val' : val}
            return out

    m = Material('m', function=eval_variable)
    rhs = Term.new('dw_volume_lvf(m.val, v)', integral, target.field.region,
                   m=m, v=v)

    eq = Equation('projection', lhs - rhs)
    eqs = Equations([eq])

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = ProblemDefinition('aux', equations=eqs, nls=nls, ls=ls)

    pb.time_update()

    # This sets the target variable with the projection solution.
    pb.solve()

    if nls_status.condition != 0:
        output('L2 projection: solver did not converge!')

def solve_stationary(conf, save_names=None, nls_status=None):

    problem = ProblemDefinition.from_conf( conf )

    problem.time_update( None )

    if save_names is not None:
        save_only( conf, save_names, problem = problem )

    state = problem.solve( nls_status = nls_status )

    return problem, state

def main():
    from sfepy.base.conf import ProblemConf, get_standard_keywords
    from sfepy.fem import ProblemDefinition
    from sfepy.base.plotutils import pylab

    required, other = get_standard_keywords()
    # Use this file as the input file.
    conf = ProblemConf.from_file( __file__, required, other )

    # Create problem instance, but do not set equations.
    problem = ProblemDefinition.from_conf( conf,
                                           init_equations = False )

    options = Struct( output_filename_trunk = None )
    
    # Solve the problem. Output is ignored, results stored by using the
    # step_hook.
    u_t = solve_branch( problem, options, linear_tension )
    u_c = solve_branch( problem, options, linear_compression )

    # Get pressure load by calling linear_*() for each time step.
    ts = problem.get_timestepper()
    load_t = nm.array( [linear_tension( ts, nm.array( [[0.0]] ) )['val']
                        for aux in ts.iter_from( 0 )],
                       dtype = nm.float64 ).squeeze()
    load_c = nm.array( [linear_compression( ts, nm.array( [[0.0]] ) )['val']
                        for aux in ts.iter_from( 0 )],
                       dtype = nm.float64 ).squeeze()

    # Join the branches.
    displacements = {}
    for key in u_t.keys():
        displacements[key] = nm.r_[u_c[key][::-1], u_t[key]]
    load = nm.r_[load_c[::-1], load_t]

    if pylab is None:
        print 'pylab cannot be imported, printing raw data!'
        print displacements
        print load
    else:
        legend = []
        for key, val in displacements.iteritems():
            pylab.plot( load, val )
            legend.append( key )

        pylab.legend( legend, loc = 2 )
        pylab.xlabel( 'tension [kPa]' )
        pylab.ylabel( 'displacement [mm]' )
        pylab.grid( True )

        pylab.gcf().savefig( 'pressure_displacement.png' )
        pylab.show()

    def __init__(self, conf, options, output_prefix, init_equations=True, **kwargs):
        """`kwargs` are passed to  ProblemDefinition.from_conf()

        Command-line options have precedence over conf.options."""
        Application.__init__(self, conf, options, output_prefix)
        self.setup_options()

        is_eqs = init_equations
        if hasattr(options, "solve_not") and options.solve_not:
            is_eqs = False
        self.problem = ProblemDefinition.from_conf(conf, init_equations=is_eqs, **kwargs)

        self.setup_output_info(self.problem, self.options)

    def test_solving(self):
        from sfepy.base.base import IndexedStruct
        from sfepy.fem \
             import FieldVariable, Material, ProblemDefinition, \
                    Function, Equation, Equations, Integral
        from sfepy.fem.conditions import Conditions, EssentialBC
        from sfepy.terms import Term
        from sfepy.solvers.ls import ScipyDirect
        from sfepy.solvers.nls import Newton

        u = FieldVariable('u', 'unknown', self.field, self.dim)
        v = FieldVariable('v', 'test', self.field, self.dim,
                          primary_var_name='u')

        m = Material('m', lam=1.0, mu=1.0)
        f = Material('f', val=[[0.02], [0.01]])

        bc_fun = Function('fix_u_fun', fix_u_fun,
                          extra_args={'extra_arg' : 'hello'})

        fix_u = EssentialBC('fix_u', self.gamma1, {'u.all' : bc_fun})
        shift_u = EssentialBC('shift_u', self.gamma2, {'u.0' : 0.1})

        integral = Integral('i', order=3)

        t1 = Term.new('dw_lin_elastic_iso(m.lam, m.mu, v, u)',
                      integral, self.omega, m=m, v=v, u=u)

        t2 = Term.new('dw_volume_lvf(f.val, v)', integral, self.omega, f=f, v=v)

        eq = Equation('balance', t1 + t2)
        eqs = Equations([eq])

        ls = ScipyDirect({})

        nls_status = IndexedStruct()
        nls = Newton({}, lin_solver=ls, status=nls_status)

        pb = ProblemDefinition('elasticity', equations=eqs, nls=nls, ls=ls)
        ## pb.save_regions_as_groups('regions')

        pb.time_update(ebcs=Conditions([fix_u, shift_u]))

        state = pb.solve()

        name = op.join(self.options.out_dir, 'test_high_level_solving.vtk')
        pb.save_state(name, state)

        ok = nls_status.condition == 0
        if not ok:
            self.report('solver did not converge!')

        _ok = state.has_ebc()
        if not _ok:
            self.report('EBCs violated!')

        ok = ok and _ok

        return ok

def make_l2_projection_data(target, eval_data):
    """
    Project scalar data given by a material-like `eval_data()` function to a
    scalar `target` field variable using the :math:`L^2` dot product.
    """
    order = target.field.approx_order * 2
    integral = Integral("i", order=order)

    un = target.name
    v = FieldVariable("v", "test", target.field, primary_var_name=un)
    lhs = Term.new("dw_volume_dot(v, %s)" % un, integral, target.field.region, v=v, **{un: target})

    def _eval_data(ts, coors, mode, **kwargs):
        if mode == "qp":
            val = eval_data(ts, coors, mode, **kwargs)
            return {"val": val}

    m = Material("m", function=_eval_data)
    rhs = Term.new("dw_volume_lvf(m.val, v)", integral, target.field.region, m=m, v=v)

    eq = Equation("projection", lhs - rhs)
    eqs = Equations([eq])

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = ProblemDefinition("aux", equations=eqs, nls=nls, ls=ls)

    pb.time_update()

    # This sets the target variable with the projection solution.
    pb.solve()

    if nls_status.condition != 0:
        output("L2 projection: solver did not converge!")

def vary_y3_size( problem ):
    """Vary size of Y3 inclusion."""
    from sfepy.fem import ProblemDefinition
    from sfepy.solvers.ts import get_print_info
    
    default_printer.prefix = 'vary_y3_size:'

    y3_diameters = [0.2, 0.25, 0.3, 0.35, 0.4]
    if diameters_g is None:
        y3_diameters = nm.linspace( 0.15, 0.45, 16 )
    else:
        y3_diameters = diameters_g
#    y3_diameters = [0.45]
    ofn_trunk, output_dir = problem.ofn_trunk, problem.output_dir
    join = os.path.join

    conf = problem.conf
    cr = conf.get_raw( 'regions' )
    n_digit, aux, d_format = get_print_info( len( y3_diameters ) + 1 )
    for ii, diameter in enumerate( y3_diameters ):
        output( 'iteration %d: diameter %3.2f' % (ii, diameter) )
        opts = problem.conf.options

        cr['Y3'] = ('nodes by select_y3_circ( x, y, z, %.5f )' % diameter, {})
        conf.edit( 'regions', cr )
        problem = ProblemDefinition.from_conf( conf )

        problem.save_regions( join( output_dir, ('regions_' + d_format) % ii ),
			      ['Y2', 'Y3'] )
        for region in problem.domain.regions:
            if not region.has_cells_if_can():
                raise ValueError( 'region %s has no cells!' % region.name )

        opts.plot_options['show'] = False
        opts.fig_name = join( output_dir,
                              (('band_gaps_%s' % d_format)
                               + '_y3_%03.2f' + fig_suffix) % (ii, diameter) )
        problem.ofn_trunk = ofn_trunk + '_' + (d_format % ii)
        out = []
        yield problem, out

        evp, bg = out[-1]

        filename = join( output_dir,
                         ('band_gaps_%s.txt' % d_format) % ii )
        log_item = '$r(Y_3)$: %f\n' % diameter
        save_log( filename, bg, log_item )

        yield None

def solve_stationary( conf, data = None, save_names = None, nls_status = None ):

    if data is None:
        # Term-dependent data.
        data = {}
    problem = ProblemDefinition.from_conf( conf )

    problem.time_update( None )

    if save_names is not None:
        save_only( conf, save_names, problem = problem )

    state = problem.solve( nls_status = nls_status )

    return problem, state, data

def recover_micro_hook( micro_filename, region, macro,
                        naming_scheme = 'step_iel',
                        recovery_file_tag='' ):

    # Create a micro-problem instance.
    required, other = get_standard_keywords()
    required.remove( 'equations' )
    pb = ProblemDefinition.from_conf_file(micro_filename,
                                          required=required,
                                          other=other,
                                          init_equations=False,
                                          init_solvers=False)

    coefs_filename = pb.conf.options.get_default_attr('coefs_filename', 'coefs')
    output_dir = pb.conf.options.get_default_attr('output_dir', '.')
    coefs_filename = op.join(output_dir, coefs_filename) + '.h5'

    # Coefficients and correctors
    coefs = Coefficients.from_file_hdf5( coefs_filename )
    corrs = get_correctors_from_file( dump_names = coefs.dump_names ) 

    recovery_hook = get_default_attr( pb.conf.options,
                                      'recovery_hook', None )

    if recovery_hook is not None:
        recovery_hook = getattr( pb.conf.funmod, recovery_hook )

        aux = max(pb.domain.shape.n_gr, 2)
        format = get_print_info( aux, fill = '0' )[1] \
            + '_' + get_print_info( pb.domain.mesh.n_el, fill = '0' )[1]

        for ig, ii, iel in region.iter_cells():
            print 'ig: %d, ii: %d, iel: %d' % (ig, ii, iel)

            local_macro = {}
            for k, v in macro.iteritems():
                local_macro[k] = v[ii,0]

            out = recovery_hook( pb, corrs, local_macro )

            # save data
            suffix = format % (ig, iel)
            micro_name = pb.get_output_name(extra='recovered_'\
                                            + recovery_file_tag + suffix)
            filename = op.join(output_dir, op.basename(micro_name))
            fpv = pb.conf.options.get_default_attr('file_per_var', False)
            pb.save_state(filename, out=out,
                          file_per_var=fpv)

def solve_direct( conf, options, problem = None, step_hook = None,
                  post_process_hook = None ):
    """Generic (simple) problem solver."""
    if problem is None:
	problem = ProblemDefinition.from_conf( conf )
	if options.output_filename_trunk:
	    ofn_trunk = options.output_filename_trunk
	    problem.ofn_trunk = ofn_trunk
	if options.output_format:
	    problem.output_format = options.output_format
    ofn_trunk = problem.ofn_trunk
    
    save_names = Struct( ebc = None, regions = None, field_meshes = None,
                         region_field_meshes = None )
    if options.save_ebc:
        save_names.ebc = ofn_trunk + '_ebc.vtk'
    if options.save_regions:
        save_names.regions = ofn_trunk + '_region'
    if options.save_field_meshes:
        save_names.field_meshes = ofn_trunk + '_field'
    if options.save_region_field_meshes:
        save_names.region_field_meshes = ofn_trunk + '_region_field'

    is_extra_save = False
    for name, val in save_names.to_dict().iteritems():
        if val is not None:
            is_extra_save = True
            break
    if is_extra_save:
        save_only( conf, save_names )

    if options.solve_not:
        return None, None, None
            
    if hasattr( conf.options, 'ts' ):
        ##
        # Time-dependent problem.
        out = solve_evolutionary_op( problem, options,
				     post_process_hook = post_process_hook,
                                     step_hook = step_hook )
    else:
        ##
        # Stationary problem.
        out = solve_stationary_op( problem, options,
				   post_process_hook = post_process_hook )

    state, data = out
    return problem, state, data