def main():
    parser = ArgumentParser()
    parser.add_argument("--version", action="version",
                        version="%(prog)s " + sfepy.__version__)
    parser.add_argument('--debug',
                        action='store_true', dest='debug',
                        default=False, help=helps['debug'])
    parser.add_argument("-o", metavar='filename', action="store",
                        dest="output_filename_trunk",
                        default=None, help=helps['filename'])
    parser.add_argument('filename_in')
    options = parser.parse_args()

    if options.debug:
        from sfepy.base.base import debug_on_error; debug_on_error()

    filename_in = options.filename_in

    required, other = get_standard_keywords()
    required.remove('equations')

    conf = ProblemConf.from_file(filename_in, required, other)

    app = HomogenizationApp(conf, options, 'homogen:')
    opts = conf.options
    if hasattr(opts, 'parametric_hook'):  # Parametric study.
        parametric_hook = conf.get_function(opts.parametric_hook)
        app.parametrize(parametric_hook)
    app()

def main():
    parser = OptionParser(usage=usage, version="%prog " + sfepy.__version__)
    parser.add_option(
        "-s", "--server", action="store_true", dest="server_mode", default=False, help=help["server_mode"]
    )
    parser.add_option("-a", "--adjoint", action="store_true", dest="adjoint", default=False, help=help["adjoint"])
    parser.add_option("-d", "--direct", action="store_true", dest="direct", default=False, help=help["direct"])
    parser.add_option(
        "-t", "--test", type=int, metavar="idsg", action="store", dest="test", default=None, help=help["test"]
    )
    parser.add_option(
        "", "--dump", metavar="filename", action="store", dest="dump_filename", default=None, help=help["dump"]
    )
    parser.add_option(
        "",
        "--pert-mesh",
        metavar="filename",
        action="store",
        dest="pert_mesh_filename",
        default=None,
        help=help["pert"],
    )
    parser.add_option("-f", "--full", action="store_true", dest="optimize", default=False, help=help["optimize"])

    options, args = parser.parse_args()

    if options.test is not None:
        options.adjoint = options.direct = True

    if options.optimize:
        options.adjoint = options.direct = False

    if (len(args) == 1) and (options.direct or options.adjoint or options.optimize):
        filename_in = args[0]
    else:
        parser.print_help(),
        return

    required, other = get_standard_keywords()
    required.remove("equations")
    if options.adjoint:
        required += ["equations_adjoint_.*", "filename_vp", "equations_direct_.*"]
        options.direct = True
    elif options.direct:
        required += ["equations_direct_.*"]
    elif options.optimize:
        required += ["equations_direct_.*", "equations_adjoint_.*", "equations_sensitivity_.*", "filename_vp"]

    conf = ProblemConf.from_file(filename_in, required, other)

    if options.direct:
        dpb, state_dp, data = solve_direct(conf, options)
    else:
        dpb, state_dp, data = None, None, None

    if options.adjoint:
        solve_adjoint(conf, options, dpb, state_dp, data)

    if options.optimize:
        solve_optimize(conf, options)

def pde_solve(conf_filename, options=None, **app_options):
    required, other = get_standard_keywords()
    conf = ProblemConf.from_file(conf_filename, required, other)

    opts = conf.options = dict_to_struct(app_options, flag=(1,)) + conf.options

    output_prefix = opts.get_default_attr('output_prefix', None)
    if output_prefix is None:
        output_prefix = output.prefix 

    if options is None:
        options = Struct(output_filename_trunk = None,
                         save_ebc = False,
                         save_regions = False,
                         save_field_meshes = False,
                         save_region_field_meshes = False,
                         save_regions_as_groups = False,
                         solve_not = False)
        
    app = SimpleApp(conf, options, output_prefix)
    if hasattr( opts, 'parametric_hook' ): # Parametric study.
        parametric_hook = getattr(conf, opts.parametric_hook)
        app.parametrize(parametric_hook)

    return app()

def main():
    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-c', '--counts',
                      action='store_true', dest='counts',
                      default=False, help=helps['counts'])
    parser.add_option('-u', '--unused',
                      action='store_true', dest='unused',
                      default=False, help=helps['unused'])
    options, args = parser.parse_args()

    if len(args) > 0:
        pdf_dir = os.path.realpath(args[0])

    else:
        parser.print_help(),
        return

    required, other = get_standard_keywords()

    terms_use = dict_from_keys_init(term_table.keys(), set)

    for filename in locate_files('*.py', pdf_dir):
        base = filename.replace(pdf_dir, '').lstrip(os.path.sep)
        output('trying "%s"...' % base)

        try:
            conf = ProblemConf.from_file(filename, required, other,
                                         verbose=False)

        except:
            output('...failed')
            continue

        use = conf.options.get('use_equations', 'equations')
        eqs_conf = getattr(conf, use)
        for key, eq_conf in eqs_conf.iteritems():
            term_descs = parse_definition(eq_conf)
            for td in term_descs:
                terms_use[td.name].add(base)

        output('...ok')
    output('...done')

    if options.unused:
        output('unused terms:')

        unused = [name for name in terms_use.keys()
                  if len(terms_use[name]) == 0]
        for name in sorted(unused):
            output('  ' + name)

        output('total: %d' % len(unused))

    else:
        output('terms use:')
        for name, ex_names in ordered_iteritems(terms_use):
            output('%s: %d' % (name, len(ex_names)))
            if not options.counts:
                for ex_name in sorted(ex_names):
                    output('  ' + ex_name)

    def test_solution(self):

        from sfepy.base.base import Struct
        from sfepy.base.conf import ProblemConf, get_standard_keywords
        from sfepy.applications import solve_pde, assign_standard_hooks
        import numpy as nm
        import os.path as op

        solutions = {}
        ok = True

        for hp, pb_filename in input_names.iteritems():

            required, other = get_standard_keywords()
            input_name = op.join(op.dirname(__file__), pb_filename)
            test_conf = ProblemConf.from_file(input_name, required, other)

            name = output_name_trunk + hp
            solver_options = Struct(output_filename_trunk=name,
                                    output_format='vtk',
                                    save_ebc=False, save_ebc_nodes=False,
                                    save_regions=False,
                                    save_regions_as_groups=False,
                                    save_field_meshes=False,
                                    solve_not=False)
            assign_standard_hooks(self, test_conf.options.get, test_conf)

            self.report( 'hyperelastic formulation: %s' % (hp, ) )

            status = NLSStatus(conditions=[])

            pb, state = solve_pde(test_conf,
                                  solver_options,
                                  nls_status=status,
                                  output_dir=self.options.out_dir,
                                  step_hook=self.step_hook,
                                  post_process_hook=self.post_process_hook,
                                  post_process_hook_final=self.post_process_hook_final)

            converged = status.condition == 0
            ok = ok and converged

            solutions[hp] = state.get_parts()['u']
            self.report('%s solved' % input_name)

        rerr = 1.0e-3
        aerr = nm.linalg.norm(solutions['TL'], ord=None) * rerr

        self.report('allowed error: rel = %e, abs = %e' % (rerr, aerr))
        ok = ok and self.compare_vectors(solutions['TL'], solutions['UL'],
                                         label1='TLF',
                                         label2='ULF',
                                         allowed_error=rerr)

        ok = ok and self.compare_vectors(solutions['UL'], solutions['ULM'],
                                         label1='ULF',
                                         label2='ULF_mixed',
                                         allowed_error=rerr)

        return ok

def main():
    from sfepy.base.base import output
    from sfepy.base.conf import ProblemConf, get_standard_keywords
    from sfepy.discrete import Problem

    output.prefix = "therel:"

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

    problem = Problem.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.set_data(state_el())
    time_solver = problem.get_time_solver()
    time_solver()

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

    def from_conf(conf, options, cls=None):
        from sfepy.base.base import Struct
        from sfepy.base.conf import ProblemConf, get_standard_keywords
        from sfepy.applications.simple_app import assign_standard_hooks

        required, other = get_standard_keywords()
        input_name = op.join(op.dirname(__file__), conf.input_name)
        test_conf = ProblemConf.from_file(input_name, required, other)

        if cls is None:
            cls = TestInput
        test = cls(test_conf=test_conf, conf=conf, options=options)

        assign_standard_hooks(test, test_conf.options.get_default_attr,
                              test_conf.funmod)

        name = op.join(test.options.out_dir, test.get_output_name_trunk())
        test.solver_options = Struct(output_filename_trunk = name,
                                     output_format ='vtk',
                                     save_ebc = False, save_regions = False,
                                     save_regions_as_groups = False,
                                     save_field_meshes = False,
                                     save_region_field_meshes = False,
                                     solve_not = False)

        return test

    def create_app(filename, is_homog=False, **kwargs):
        from sfepy.base.conf import ProblemConf, get_standard_keywords
        from sfepy.homogenization.homogen_app import HomogenizationApp
        from sfepy.applications import PDESolverApp

        required, other = get_standard_keywords()
        if is_homog:
            required.remove('equations')

        conf = ProblemConf.from_file(filename, required, other,
                                     define_args=kwargs)
        options = Struct(output_filename_trunk=None,
                         save_ebc=False,
                         save_ebc_nodes=False,
                         save_regions=False,
                         save_regions_as_groups=False,
                         save_field_meshes=False,
                         solve_not=False)

        if is_homog:
            app = HomogenizationApp(conf, options, 'material_opt_micro:')

        else:
            app = PDESolverApp(conf, options, 'material_opt_macro:')

        app.conf.opt_data = {}
        opts = conf.options
        if hasattr(opts, 'parametric_hook'):  # Parametric study.
            parametric_hook = conf.get_function(opts.parametric_hook)
            app.parametrize(parametric_hook)

        return app

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 main():
    parser = OptionParser(usage=usage, version="%prog " + sfepy.__version__)
    parser.add_option("-o", "", metavar='filename', action="store",
                      dest="output_filename_trunk",
                      default=None, help=help['filename'])

    (options, args) = parser.parse_args()

    if (len(args) == 1):
        filename_in = args[0]
    else:
        parser.print_help(),
        return

    required, other = get_standard_keywords()
    required.remove('equations')

    conf = ProblemConf.from_file(filename_in, required, other)

    app = HomogenizationApp(conf, options, 'homogen:')
    opts = conf.options
    if hasattr(opts, 'parametric_hook'):  # Parametric study.
        parametric_hook = conf.get_function(opts.parametric_hook)
        app.parametrize(parametric_hook)
    app()

def main():
    from sfepy.base.base import output
    from sfepy.base.conf import ProblemConf, get_standard_keywords
    from sfepy.discrete import Problem
    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 = Problem.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:
        output('matplotlib cannot be imported, printing raw data!')
        output(displacements)
        output(load)
    else:
        legend = []
        for key, val in six.iteritems(displacements):
            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 assemble_matrices(define, mod, pars, set_wave_dir, options):
    """
    Assemble the blocks of dispersion eigenvalue problem matrices.
    """
    define_problem = functools.partial(define,
                                       filename_mesh=options.mesh_filename,
                                       pars=pars,
                                       approx_order=options.order,
                                       refinement_level=options.refine,
                                       solver_conf=options.solver_conf,
                                       plane=options.plane,
                                       post_process=options.post_process)

    conf = ProblemConf.from_dict(define_problem(), mod)

    pb = Problem.from_conf(conf)
    pb.dispersion_options = options
    pb.set_output_dir(options.output_dir)
    dim = pb.domain.shape.dim

    # Set the normalized wave vector direction to the material(s).
    wdir = nm.asarray(options.wave_dir[:dim], dtype=nm.float64)
    wdir = wdir / nm.linalg.norm(wdir)
    set_wave_dir(pb, wdir)

    bbox = pb.domain.mesh.get_bounding_box()
    size = (bbox[1] - bbox[0]).max()
    scaling0 = apply_unit_multipliers([1.0], ['length'],
                                      options.unit_multipliers)[0]
    scaling = scaling0
    if options.mesh_size is not None:
        scaling *= options.mesh_size / size
    output('scaling factor of periodic cell mesh coordinates:', scaling)
    output('new mesh size with applied unit multipliers:', scaling * size)
    pb.domain.mesh.coors[:] *= scaling
    pb.set_mesh_coors(pb.domain.mesh.coors, update_fields=True)

    bzone = 2.0 * nm.pi / (scaling * size)
    output('1. Brillouin zone size:', bzone * scaling0)
    output('1. Brillouin zone size with applied unit multipliers:', bzone)

    pb.time_update()
    pb.update_materials()

    # Assemble the matrices.
    mtxs = {}
    for key, eq in pb.equations.iteritems():
        mtxs[key] = mtx = pb.mtx_a.copy()
        mtx = eq.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx)
        mtx.eliminate_zeros()
        output_array_stats(mtx.data, 'nonzeros in %s' % key)

        output('symmetry checks:')
        output('%s - %s^T:' % (key, key), max_diff_csr(mtx, mtx.T))
        output('%s - %s^H:' % (key, key), max_diff_csr(mtx, mtx.H))

    return pb, wdir, bzone, mtxs

def get_homog_coefs_linear(ts, coor, mode,
                           micro_filename=None, regenerate=False,
                           coefs_filename=None, define_args=None):

    oprefix = output.prefix
    output.prefix = 'micro:'

    required, other = get_standard_keywords()
    required.remove( 'equations' )

    conf = ProblemConf.from_file(micro_filename, required, other,
                                 verbose=False, define_args=define_args)
    if coefs_filename is None:
        coefs_filename = conf.options.get('coefs_filename', 'coefs')
        coefs_filename = op.join(conf.options.get('output_dir', '.'),
                                 coefs_filename) + '.h5'

    if not regenerate:
        if op.exists( coefs_filename ):
            if not pt.is_hdf5_file( coefs_filename ):
                regenerate = True
        else:
            regenerate = True

    if regenerate:
        options = Struct( output_filename_trunk = None )

        app = HomogenizationApp( conf, options, 'micro:' )
        coefs = app()
        if type(coefs) is tuple:
            coefs = coefs[0]

        coefs.to_file_hdf5( coefs_filename )
    else:
        coefs = Coefficients.from_file_hdf5( coefs_filename )

    out = {}
    if mode == None:
        for key, val in six.iteritems(coefs.__dict__):
            out[key] = val

    elif mode == 'qp':
        for key, val in six.iteritems(coefs.__dict__):
            if type( val ) == nm.ndarray or type(val) == nm.float64:
                out[key] = nm.tile( val, (coor.shape[0], 1, 1) )
            elif type(val) == dict:
                for key2, val2 in six.iteritems(val):
                    if type(val2) == nm.ndarray or type(val2) == nm.float64:
                        out[key+'_'+key2] = \
                                          nm.tile(val2, (coor.shape[0], 1, 1))

    else:
        out = None

    output.prefix = oprefix

    return out

def main():
    version = open( op.join( init_sfepy.install_dir,
                             'VERSION' ) ).readlines()[0][:-1]

    parser = OptionParser( usage = usage, version = "%prog " + version )
    parser.add_option( "-o", "", metavar = 'filename',
                       action = "store", dest = "output_filename_trunk",
                       default = None, help = help['filename'] )
    parser.add_option( "", "--format", metavar = 'format',
                       action = "store", dest = "output_format",
                       default = "vtk", help = help['output_format'] )
    parser.add_option( "", "--save-ebc",
                       action = "store_true", dest = "save_ebc",
                       default = False, help = help['save_ebc'] )
    parser.add_option( "", "--save-regions",
                       action = "store_true", dest = "save_regions",
                       default = False, help = help['save_regions'] )
    parser.add_option( "", "--save-field-meshes",
                       action = "store_true", dest = "save_field_meshes",
                       default = False, help = help['save_field_meshes'] )
    parser.add_option( "", "--save-region-field-meshes",
                       action = "store_true", dest = "save_region_field_meshes",
                       default = False, help = help['save_region_field_meshes'] )
    parser.add_option( "", "--solve-not",
                       action = "store_true", dest = "solve_not",
                       default = False, help = help['solve_not'] )
    parser.add_option( "", "--list", metavar = 'what',
                       action = "store", dest = "_list",
                       default = None, help = help['list'] )

    options, args = parser.parse_args()
#    print options; pause()

    if (len( args ) == 1):
        filename_in = args[0];
    else:
        if options._list == 'terms':
            print_terms()
        else:
            parser.print_help(),
        return
    
    required, other = get_standard_keywords()
    if options.solve_not:
        required.remove( 'equations' )
        required.remove( 'solver_[0-9]+|solvers' )
        other.extend( ['equations'] )

    conf = ProblemConf.from_file( filename_in, required, other )
    opts = conf.options
    output_prefix = get_default_attr( opts, 'output_prefix', 'sfepy:' )

    app = SimpleApp( conf, options, output_prefix )
    if hasattr( opts, 'parametric_hook' ): # Parametric study.
        parametric_hook = getattr( conf, opts.parametric_hook )
        app.parametrize( parametric_hook )
    app()

    def test_stokes_slip_bc(self):
        import scipy.sparse as sp

        from sfepy.base.conf import ProblemConf
        from sfepy.discrete import Problem
        import examples.navier_stokes.stokes_slip_bc as ssb

        conf = ProblemConf.from_module(ssb)
        pb = Problem.from_conf(conf, init_solvers=False)
        pb.time_update()
        variables = pb.get_variables()

        adi = variables.adi
        lcdi = variables.lcdi
        mtx = variables.mtx_lcbc

        ok = adi.var_names == lcdi.var_names
        self.report('same adi-lcdi ordering:', ok)

        ublock = mtx[adi.indx['u']]
        ir, ic = ublock.nonzero()
        ir += adi.indx['u'].start

        i0, i1 = adi.indx['u'].start, adi.indx['u'].stop
        _ok0 = (i0 <= ir).all() and  (ir < i1).all()
        self.report('u block rows in [%d %d[: %s' % (i0, i1, _ok0))

        i0, i1 = lcdi.indx['u'].start, lcdi.indx['u'].stop
        _ok1 = (i0 <= ic).all() and  (ic < i1).all()
        self.report('u block cols in [%d %d[: %s' % (i0, i1, _ok1))

        ok = ok and _ok0 and _ok1

        pblock = mtx[adi.indx['p']]
        ir, ic, iv = sp.find(pblock)
        ir += adi.indx['p'].start

        i0, i1 = adi.indx['p'].start, adi.indx['p'].stop
        _ok0 = (i0 <= ir).all() and  (ir < i1).all()
        self.report('p block rows in [%d %d[: %s' % (i0, i1, _ok0))

        i0, i1 = lcdi.indx['p'].start, lcdi.indx['p'].stop
        _ok1 = (i0 <= ic).all() and  (ic < i1).all()
        self.report('p block cols in [%d %d[: %s' % (i0, i1, _ok1))

        ok = ok and _ok0 and _ok1

        _ok0 = (len(ir) == adi.n_dof['p'])
        self.report('p block size correct:', _ok0)
        _ok1 = ((ir - adi.indx['p'].start) == (ic - lcdi.indx['p'].start)).all()
        self.report('p block diagonal:', _ok1)
        _ok2 = (iv == 1.0).all()
        self.report('p block identity:', _ok2)

        ok = ok and _ok0 and _ok1 and _ok2

        return ok

    def from_conf( conf, options, cls = None ):
        from sfepy.base.conf import ProblemConf, get_standard_keywords

        required, other = get_standard_keywords()
        test_conf = ProblemConf.from_file( conf.input_name, required, other )

        if cls is None:
            cls = TestInput
        test = cls( test_conf = test_conf, conf = conf, options = options )

        return test

def solve_pde(conf, options=None, status=None, **app_options):
    """
    Solve a system of partial differential equations (PDEs).

    This function is a convenience wrapper that creates and runs an instance of
    :class:`PDESolverApp`.

    Parameters
    ----------
    conf : str or ProblemConf instance
        Either the name of the problem description file defining the PDEs,
        or directly the ProblemConf instance.
    options : options
        The command-line options.
    status : dict-like
        The object for storing the solver return status.
    app_options : kwargs
        The keyword arguments that can override application-specific options.
    """
    if not isinstance(conf, ProblemConf):
        required, other = get_standard_keywords()
        conf = ProblemConf.from_file(conf, required, other)

    opts = conf.options = (dict_to_struct(app_options, flag=(1,),
                                          constructor=type(conf.options))
                           + conf.options)

    output_prefix = opts.get('output_prefix', None)
    if output_prefix is None:
        output_prefix = output.prefix

    if options is None:
        options = Struct(output_filename_trunk=None,
                         save_ebc=False,
                         save_ebc_nodes=False,
                         save_regions=False,
                         save_field_meshes=False,
                         save_regions_as_groups=False,
                         solve_not=False)

    if conf.options.get('evps') is None:
        app = PDESolverApp(conf, options, output_prefix)

    else:
        from .evp_solver_app import EVPSolverApp

        app = EVPSolverApp(conf, options, output_prefix)

    if hasattr(opts, 'parametric_hook'): # Parametric study.
        parametric_hook = conf.get_function(opts.parametric_hook)
        app.parametrize(parametric_hook)

    return app(status=status)

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 main():
    parser = ArgumentParser()
    parser.add_argument("--version", action="version",
                        version="%(prog)s " + sfepy.__version__)
    parser.add_argument('--debug',
                        action='store_true', dest='debug',
                        default=False, help=helps['debug'])
    parser.add_argument("-o", metavar='filename',
                        action="store", dest="output_filename_trunk",
                        default=None, help=helps['filename'])
    parser.add_argument("-b", "--band-gaps",
                        action="store_true", dest="detect_band_gaps",
                        default=False, help=helps['detect_band_gaps'])
    parser.add_argument("-d", "--dispersion",
                        action="store_true", dest="analyze_dispersion",
                        default=False, help=helps['analyze_dispersion'])
    parser.add_argument("-p", "--plot",
                        action="store_true", dest="plot",
                        default=False, help=helps['plot'])
    parser.add_argument("--phase-velocity",
                        action="store_true", dest="phase_velocity",
                        default=False, help=helps['phase_velocity'])
    parser.add_argument("filename_in")
    options = parser.parse_args()

    if options.debug:
        from sfepy.base.base import debug_on_error; debug_on_error()

    if options.plot:
        if plt is None:
            output('matplotlib.pyplot cannot be imported, ignoring option -p!')
            options.plot = False
        elif options.analyze_dispersion == False:
            options.detect_band_gaps = True

    required, other = get_standard_keywords()
    required.remove('equations')
    if not options.analyze_dispersion:
        required.remove('solver_[0-9]+|solvers')
    if options.phase_velocity:
        required.remove('ebc_[0-9]+|ebcs')
    conf = ProblemConf.from_file(options.filename_in, required, other)

    app = AcousticBandGapsApp(conf, options, 'phonon:')
    opts = conf.options
    if hasattr(opts, 'parametric_hook'): # Parametric study.
        parametric_hook = conf.get_function(opts.parametric_hook)
        app.parametrize(parametric_hook)
    app()

def main():
    parser = OptionParser(usage=usage, version="%prog " + sfepy.__version__)
    parser.add_option("-o", "", metavar='filename',
                      action="store", dest="output_filename_trunk",
                      default=None, help=help['filename'])
    parser.add_option("-b", "--band-gaps",
                      action="store_true", dest="detect_band_gaps",
                      default=False, help=help['detect_band_gaps'])
    parser.add_option("-d", "--dispersion",
                      action="store_true", dest="analyze_dispersion",
                      default=False, help=help['analyze_dispersion'])
    parser.add_option("-p", "--plot",
                      action="store_true", dest="plot",
                      default=False, help=help['plot'])
    parser.add_option("--phase-velocity",
                      action="store_true", dest="phase_velocity",
                      default=False, help=help['phase_velocity'])

    options, args = parser.parse_args()
    if options.plot:
        if plt is None:
            output('matplotlib.pyplot cannot be imported, ignoring option -p!')
            options.plot = False
        elif options.analyze_dispersion == False:
            options.detect_band_gaps = True

    if (len(args) == 1):
        filename_in = args[0];
    else:
        parser.print_help(),
        return

    required, other = get_standard_keywords()
    required.remove('equations')
    if not options.analyze_dispersion:
        required.remove('solver_[0-9]+|solvers')
    if options.phase_velocity:
        required.remove('ebc_[0-9]+|ebcs')
    conf = ProblemConf.from_file(filename_in, required, other)

    app = AcousticBandGapsApp(conf, options, 'phonon:')
    opts = conf.options
    if hasattr(opts, 'parametric_hook'): # Parametric study.
        parametric_hook = conf.get_function(opts.parametric_hook)
        app.parametrize(parametric_hook)
    app()