def process_options(options):
        """
        Application options setup. Sets default values for missing
        non-compulsory options.
        """
        get = options.get_default_attr

        output_dir = get('output_dir', '.')
        if not os.path.exists(output_dir):
            os.makedirs(output_dir)

        return Struct(save_results=get('save_results', True),
                      # Save each variable into a separate file, using
                      # the region of its definition only.
                      linearization=dict_to_struct(get('linearization',
                                                       {'kind' : 'strip'})),
                      file_per_var=get('file_per_var', False),
                      output_format=get('output_format', 'vtk'),
                      output_dir=output_dir,
                      # Called after each time step, can do anything, no
                      # return value.
                      step_hook=get('step_hook', None),
                      # Called after each time step.
                      post_process_hook=get('post_process_hook', None),
                      # Called after all time steps, or in the
                      # stationary case.
                      post_process_hook_final=get('post_process_hook_final',
                                                  None),
                      # Called in init process.
                      pre_process_hook=get('pre_process_hook', None),
                      use_equations=get('use_equations', 'equations'))

    def process_options(options):
        """
        Application options setup. Sets default values for missing
        non-compulsory options.
        """
        get = options.get

        output_dir = get("output_dir", ".")
        if not os.path.exists(output_dir):
            os.makedirs(output_dir)

        return Struct(
            save_results=get("save_results", True),
            # Save each variable into a separate file, using
            # the region of its definition only.
            linearization=dict_to_struct(get("linearization", {"kind": "strip"})),
            file_per_var=get("file_per_var", False),
            output_format=get("output_format", "vtk"),
            output_dir=output_dir,
            # Called after each time step, can do anything, no
            # return value.
            step_hook=get("step_hook", None),
            # Called after each time step.
            post_process_hook=get("post_process_hook", None),
            # Called after all time steps, or in the
            # stationary case.
            post_process_hook_final=get("post_process_hook_final", None),
            # Called in init process.
            pre_process_hook=get("pre_process_hook", None),
            use_equations=get("use_equations", "equations"),
        )

    def process_options(self):
        options = dict_to_struct(self.options)
        get = options.get

        return Struct(var_name=get('var_name', None,
                                   'missing "var_name" in options!'),
                      threshold=get('threshold', 1e-4),
                      threshold_is_relative=get('threshold_is_relative', True),
                      transform=get('transform', None))

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 parse_linearization(linearization):
    out = {}
    for item in linearization.split(","):
        key, val = item.split(":")
        if key == "eps":
            val = float(val)
        elif key in ("min_level", "max_level"):
            val = int(val)
        elif key == "kind":
            pass
        else:
            raise ValueError("wrong linearization option key! (%s)" % key)
        out[key] = val

    return dict_to_struct(out)

def parse_linearization(linearization):
    out = {}
    for item in linearization.split(','):
        key, val = item.split(':')
        if key == 'eps':
            val = float(val)
        elif key in ('min_level', 'max_level'):
            val = int(val)
        elif key == 'kind':
            pass
        else:
            raise ValueError('wrong linearization option key! (%s)'
                             % key)
        out[key] = val

    return dict_to_struct(out)

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_regions_as_groups = False,
                         solve_not = False)
        
    app = SimpleApp(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()

def transform_to_struct_10(adict):
    return dict_to_struct(adict, flag=(1,0))

def transform_to_struct_01(adict):
    return dict_to_struct(adict, flag=(0,1))

def transform_to_i_struct_1(adict):
    return dict_to_struct(adict, flag=(1,), constructor=IndexedStruct)

def create_stabil_mat(problem):
    """Using the stabilization material stub make it the true material."""
    from sfepy.base.base import dict_to_struct, debug
    from sfepy.fem.functions import Function

    # Identity map...
    ns = {'p' : 'p', 'q' : 'q',
          'u' : 'u', 'b' : 'b', 'v' : 'v',
          'fluid' : 'fluid', 'omega' : 'omega', 'i1' : 'i1', 'i2' : 'i2'}

    variables = problem.get_variables()

    # Indices to the state vector.
    ii = {}
    ii['u'] = variables.get_indx(ns['u'])
    ii['us'] = variables.get_indx(ns['u'], stripped=True)
    ii['ps'] = variables.get_indx(ns['p'], stripped=True)

    materials = problem.get_materials()
    stabil_mat = materials['stabil']
    stabil = dict_to_struct(stabil_mat.datas['special'], flag=(1,))

    # The viscosity.
    fluid_mat = materials[ns['fluid']]
    viscosity = fluid_mat.function()['viscosity']

    # The Friedrich's constant.
    c_friedrichs = problem.domain.get_diameter()
    sigma = 1e-12 # 1 / dt.

    # Element diameter modes.
    diameter_modes = {'edge' : 0, 'volume' : 1, 'max' : 2}

    def mat_fun(ts, coor, mode=None, term=None, problem=None,
                b_norm=1.0, **kwargs):
        if mode != 'qp': return

        region = term.region

        print '|b|_max (mat_fun):', b_norm
        gamma = viscosity + b_norm * c_friedrichs

        data = {}
        if stabil.gamma is None:
            data['gamma'] = stabil.gamma_mul * gamma
        else:
            data['gamma'] = nm.asarray( stabil.gamma_mul * stabil.gamma,
                                        dtype = nm.float64 )
        data['gamma'] = nm.tile(data['gamma'], (coor.shape[0], 1, 1))

        if stabil.delta is None:
            term = problem.equations['balance'].terms['dw_lin_convect']
            for ig in term.iter_groups():
                # This sets term.ig - for 1 group only!!!
                break
            var = variables[ns['u']]
            ap, vg = term.get_approximation(var)
            delta = 1.0
            mode = diameter_modes[stabil.diameter_mode]
            cells = region.get_cells( ig )
            diameters2 = problem.domain.get_element_diameters( ig, cells, vg,
                                                             mode )
            val1 = min( 1.0, 1.0 / sigma )
            val2 = sigma * c_friedrichs**2
            val3 = (b_norm**2) * min( (c_friedrichs**2) / viscosity, 1.0 / sigma )
#            print val1, gamma, val2, val3
            delta = stabil.delta_mul * val1 * diameters2 / (gamma + val2 + val3)

            n_qp = coor.shape[0] / diameters2.shape[0]
            data['diameters2'] = nm.repeat(diameters2, n_qp)
            data['diameters2'].shape = data['diameters2'].shape + (1, 1)

            data['delta'] = nm.repeat(delta, n_qp)
            data['delta'].shape = data['delta'].shape + (1, 1)
        else:
            val = stabil.delta_mul * stabil.delta
            data['delta'] = nm.tile(data['delta'], (coor.shape[0], 1, 1))
        
        if stabil.tau is None:
            data['tau'] = stabil.tau_red * data['delta']
        else:
            data['tau'] = nm.asarray( stabil.tau_mul * stabil.tau,
                                      dtype = nm.float64 )
            data['tau'] = nm.tile(data['tau'], (coor.shape[0], 1, 1))

        return data

    stabil_mat.set_function(Function('stabil', mat_fun))

    return stabil_mat, ns, ii

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('-d', '--dump', action='store_true', dest='dump',
                       default=False, help=help['dump'])
    parser.add_option('', '--same-dir', action='store_true', dest='same_dir',
                      default=False, help=help['same_dir'])
    parser.add_option('-l', '--linearization', metavar='options',
                      action='store', dest='linearization',
                      default=None, help=help['linearization'])
    parser.add_option('', '--times', action='store_true', dest='times',
                      default=False, help=help['times'])
    parser.add_option('-f', '--from', type=int, metavar='ii',
                      action='store', dest='step_from',
                      default=0, help=help['from'])
    parser.add_option('-t', '--to', type=int, metavar='ii',
                      action='store', dest='step_to',
                      default=None, help=help['to'])
    parser.add_option('-s', '--step', type=int, metavar='ii',
                      action='store', dest='step_by',
                      default=1, help=help['step'])
    parser.add_option('-e', '--extract', metavar='list',
                      action='store', dest='extract',
                      default=None, help=help['extract'])
    parser.add_option('-a', '--average', action='store_true', dest='average',
                      default=False, help=help['average'])

    (options, args) = parser.parse_args()

    nargs = len(args)
    if nargs == 1:
        filename_results = args[0]
        linearize = False

    elif nargs == 2:
        filename_in, filename_results = args
        linearize = True
        options.dump = True

    else:
        parser.print_help()
        return

    if options.times:
        steps, times, nts, dts = th.extract_times(filename_results)
        for ii, time in enumerate(times):
            step = steps[ii]
            print('%d %e %e %e' % (step, time, nts[ii], dts[ii]))

    if options.dump:
        trunk = get_default(options.output_filename_trunk,
                            get_trunk(filename_results))
        if options.same_dir:
            trunk = os.path.join(os.path.dirname(filename_results),
                                 os.path.basename(trunk))

        args = {}
        if linearize:
            problem = create_problem(filename_in)

            linearization = Struct(kind='adaptive', min_level=0,
                                   max_level=2, eps=1e-2)
            aux = problem.conf.options.get('linearization', None)
            linearization.update(aux)

            if options.linearization is not None:
                aux = parse_linearization(options.linearization)
                linearization.update(aux)

            args.update({'fields' : problem.fields,
                         'linearization' : linearization})

        if options.step_to is None:
            args.update({'step0' : options.step_from})

        else:
            args.update({'steps' : nm.arange(options.step_from,
                                             options.step_to + 1,
                                             options.step_by, dtype=nm.int)})

        th.dump_to_vtk(filename_results, output_filename_trunk=trunk, **args)

    if options.extract:
        ths, ts = th.extract_time_history(filename_results, options.extract)

        if options.average:
            ths = th.average_vertex_var_in_cells(ths)

        if options.output_filename_trunk:
            th.save_time_history(ths, ts, options.output_filename_trunk + '.h5')

        else:
            print(dict_to_struct(ths, flag=(1, 1, 1)).str_all())

def main():
    parser = ArgumentParser(description=__doc__, formatter_class=RawDescriptionHelpFormatter)
    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("-d", "--dump", action="store_true", dest="dump", default=False, help=helps["dump"])
    parser.add_argument("--same-dir", action="store_true", dest="same_dir", default=False, help=helps["same_dir"])
    parser.add_argument(
        "-l",
        "--linearization",
        metavar="options",
        action="store",
        dest="linearization",
        default=None,
        help=helps["linearization"],
    )
    parser.add_argument("--times", action="store_true", dest="times", default=False, help=helps["times"])
    parser.add_argument(
        "-f", "--from", type=int, metavar="ii", action="store", dest="step_from", default=0, help=helps["from"]
    )
    parser.add_argument(
        "-t", "--to", type=int, metavar="ii", action="store", dest="step_to", default=None, help=helps["to"]
    )
    parser.add_argument(
        "-s", "--step", type=int, metavar="ii", action="store", dest="step_by", default=1, help=helps["step"]
    )
    parser.add_argument(
        "-e", "--extract", metavar="list", action="store", dest="extract", default=None, help=helps["extract"]
    )
    parser.add_argument("-a", "--average", action="store_true", dest="average", default=False, help=helps["average"])
    parser.add_argument("input_file", nargs="?", default=None)
    parser.add_argument("results_file")
    options = parser.parse_args()

    if options.debug:
        from sfepy.base.base import debug_on_error

        debug_on_error()

    filename_in = options.input_file
    filename_results = options.results_file

    if filename_in is None:
        linearize = False
    else:
        linearize = True
        options.dump = True

    if options.times:
        steps, times, nts, dts = th.extract_times(filename_results)
        for ii, time in enumerate(times):
            step = steps[ii]
            print("%d %e %e %e" % (step, time, nts[ii], dts[ii]))

    if options.dump:
        trunk = get_default(options.output_filename_trunk, get_trunk(filename_results))
        if options.same_dir:
            trunk = os.path.join(os.path.dirname(filename_results), os.path.basename(trunk))

        args = {}
        if linearize:
            problem = create_problem(filename_in)

            linearization = Struct(kind="adaptive", min_level=0, max_level=2, eps=1e-2)
            aux = problem.conf.options.get("linearization", None)
            linearization.update(aux)

            if options.linearization is not None:
                aux = parse_linearization(options.linearization)
                linearization.update(aux)

            args.update({"fields": problem.fields, "linearization": linearization})

        if options.step_to is None:
            args.update({"step0": options.step_from})

        else:
            args.update({"steps": nm.arange(options.step_from, options.step_to + 1, options.step_by, dtype=nm.int)})

        th.dump_to_vtk(filename_results, output_filename_trunk=trunk, **args)

    if options.extract:
        ths, ts = th.extract_time_history(filename_results, options.extract)

        if options.average:
            ths = th.average_vertex_var_in_cells(ths)

        if options.output_filename_trunk:
            th.save_time_history(ths, ts, options.output_filename_trunk + ".h5")

        else:
            print(dict_to_struct(ths, flag=(1, 1, 1)).str_all())

def create_stabil_mat(problem):
    """Using the stabilization material stub make it the true material."""
    from sfepy.base.base import dict_to_struct, debug
    from sfepy.fem.functions import Function

    # Identity map...
    ns = {"p": "p", "q": "q", "u": "u", "b": "b", "v": "v", "fluid": "fluid", "omega": "omega", "i1": "i1", "i2": "i2"}

    variables = problem.get_variables()

    # Indices to the state vector.
    ii = {}
    ii["u"] = variables.get_indx(ns["u"])
    ii["us"] = variables.get_indx(ns["u"], stripped=True)
    ii["ps"] = variables.get_indx(ns["p"], stripped=True)

    stabil_mat = problem.materials["stabil"]
    stabil = dict_to_struct(stabil_mat.datas["special"], flag=(1,))

    # The viscosity.
    fluid_mat = problem.materials[ns["fluid"]]
    viscosity = fluid_mat.function()["viscosity"]

    # The Friedrich's constant.
    c_friedrichs = problem.domain.get_diameter()
    sigma = 1e-12  # 1 / dt.

    # Element diameter modes.
    diameter_modes = {"edge": 0, "volume": 1, "max": 2}

    def mat_fun(ts, coor, mode=None, region=None, ig=None, b_norm=1.0):
        if mode != "qp":
            return

        print "|b|_max (mat_fun):", b_norm
        gamma = viscosity + b_norm * c_friedrichs

        data = {}
        if stabil.gamma is None:
            data["gamma"] = stabil.gamma_mul * gamma
        else:
            data["gamma"] = nm.asarray(stabil.gamma_mul * stabil.gamma, dtype=nm.float64)
        data["gamma"] = nm.tile(data["gamma"], (coor.shape[0], 1, 1))

        if stabil.delta is None:
            term = problem.equations["balance"].terms["dw_lin_convect"]
            for ig in term.iter_groups():
                # This sets term.ig - for 1 group only!!!
                break
            var = variables[ns["u"]]
            ap, vg = term.get_approximation(var)
            delta = 1.0
            mode = diameter_modes[stabil.diameter_mode]
            cells = region.get_cells(ig)
            diameters2 = problem.domain.get_element_diameters(ig, cells, vg, mode)
            val1 = min(1.0, 1.0 / sigma)
            val2 = sigma * c_friedrichs ** 2
            val3 = (b_norm ** 2) * min((c_friedrichs ** 2) / viscosity, 1.0 / sigma)
            #            print val1, gamma, val2, val3
            delta = stabil.delta_mul * val1 * diameters2 / (gamma + val2 + val3)

            n_qp = coor.shape[0] / diameters2.shape[0]
            data["diameters2"] = nm.repeat(diameters2, n_qp)
            data["diameters2"].shape = data["diameters2"].shape + (1, 1)

            data["delta"] = nm.repeat(delta, n_qp)
            data["delta"].shape = data["delta"].shape + (1, 1)
        else:
            val = stabil.delta_mul * stabil.delta
            data["delta"] = nm.tile(data["delta"], (coor.shape[0], 1, 1))

        if stabil.tau is None:
            data["tau"] = stabil.tau_red * data["delta"]
        else:
            data["tau"] = nm.asarray(stabil.tau_mul * stabil.tau, dtype=nm.float64)
            data["tau"] = nm.tile(data["tau"], (coor.shape[0], 1, 1))

        return data

    stabil_mat.set_function(Function("stabil", mat_fun))

    return stabil_mat, ns, ii

    def test_semismooth_newton(self):
        import numpy as nm
        from sfepy.solvers import Solver

        ns = [0, 6, 2, 2]

        offsets = nm.cumsum(ns)
        nx = offsets[-1]

        iw = slice(offsets[0], offsets[1])
        ig = slice(offsets[1], offsets[2])
        il = slice(offsets[2], offsets[3])

        def fun_smooth(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]
            xl = vec_x[il]

            m = self.ms
            rw = m['A'] * xw - m['Bb'].T * xg - self.m['fs'] 
            rg = m['Bb'] * xw + xl * xg

            rwg = nm.r_[rw, rg]

            return rwg

        def fun_smooth_grad(vec_x):
            xw = vec_x[iw]
            xl = vec_x[il]
            xg = vec_x[ig]

            m = self.m

            mzl = nm.zeros((6, 2), dtype=nm.float64)

            mw = nm.c_[m['A'], -m['Bb'].T, mzl]
            mg = nm.c_[m['Bb'], nm.diag(xl), nm.diag(xg)]

            mx = nm.r_[mw, mg]

            mx = sps.csr_matrix(mx)
            return mx

        def fun_a(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]

            subsd = {}
            for ii, key in enumerate(self.w_names):
                subsd[key] = xw[ii]

            sn = eval_matrix(self.m['sn'], **subsd).squeeze()

            ra = nm.abs(xg) - fc * nm.abs(sn)

            return -ra

        def fun_a_grad(vec_x):
            xw = vec_x[iw]
            xg = vec_x[ig]
            xl = vec_x[il]

            subsd = {}
            for ii, key in enumerate(self.w_names):
                subsd[key] = xw[ii]

            md = eval_matrix(self.m['D'], **subsd)
            sn = eval_matrix(self.m['sn'], **subsd).squeeze()

            ma = nm.zeros((xl.shape[0], nx), dtype=nm.float64)
            ma[:,iw] = - fc * nm.sign(sn)[:,None] * md
            ma[:,ig] = nm.sign(xg)[:,None] * self.m['C']

            return -sps.csr_matrix(ma)

        def fun_b(vec_x):
            xl = vec_x[il]

            return xl

        def fun_b_grad(vec_x):
            xl = vec_x[il]

            mb = nm.zeros((xl.shape[0], nx), dtype=nm.float64)
            mb[:,il] = self.m['C']

            return sps.csr_matrix(mb)

        vec_x0 = 0.1 * nm.ones((nx,), dtype=nm.float64)

        lin_solver = Solver.any_from_conf(dict_to_struct(ls_conf))
        status = {}
        solver = Solver.any_from_conf(dict_to_struct(conf),
                                      fun_smooth=fun_smooth,
                                      fun_smooth_grad=fun_smooth_grad,
                                      fun_a=fun_a,
                                      fun_a_grad=fun_a_grad,
                                      fun_b=fun_b,
                                      fun_b_grad=fun_b_grad,
                                      lin_solver=lin_solver,
#.........这里部分代码省略.........

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( "-d", "--dump",
                       action = "store_true", dest = "dump",
                       default = False, help = help['dump'] )
    parser.add_option( "", "--same-dir",
                       action = "store_true", dest = "same_dir",
                       default = False, help = help['same_dir'] )
    parser.add_option( "-f", "--from", type = int, metavar = 'ii',
                       action = "store", dest = "step_from",
                       default = 0, help = help['from'] )
    parser.add_option( "-t", "--to", type = int, metavar = 'ii',
                       action = "store", dest = "step_to",
                       default = None, help = help['to'] )
    parser.add_option( "-s", "--step", type = int, metavar = 'ii',
                       action = "store", dest = "step_by",
                       default = 1, help = help['step'] )
    parser.add_option( "-e", "--extract", metavar = 'list',
                       action = "store", dest = "extract",
                       default = None, help = help['extract'] )
    parser.add_option( "-a", "--average",
                       action = "store_true", dest = "average",
                       default = False, help = help['average'] )

    (options, args) = parser.parse_args()

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

    if options.dump:
        trunk = get_default(options.output_filename_trunk,
                            get_trunk(filename_in))
        if options.same_dir:
            trunk = os.path.join(os.path.dirname(filename_in),
                                 os.path.basename(trunk))
        
        if options.step_to is None:
            dump_to_vtk(filename_in,
                        output_filename_trunk=trunk,
                        step0=options.step_from)

        else:
            dump_to_vtk(filename_in,
                        output_filename_trunk=trunk,
                        steps=nm.arange(options.step_from,
                                        options.step_to + 1,
                                        options.step_by, dtype=nm.int))

    if options.extract:
        ths, ts = extract_time_history(filename_in, options.extract)
##         print ths

        if options.average:
            ths = average_vertex_var_in_cells( ths )
##             print ths

        if options.output_filename_trunk:
            save_time_history(ths, ts, options.output_filename_trunk + '.h5')

        else:
            print dict_to_struct(ths, flag=(1, 1, 1)).str_all()