def save_image(self, filename):
        """Save a snapshot of the current scene."""
        name = os.path.join(self.output_dir, filename)
        ensure_path(name)

        output('saving %s...' % name)
        self.scene.scene.save(name)
        output('...done')

def generate_thumbnails(thumbnails_dir, images_dir, scale=0.3):
    """
    Generate thumbnails into `thumbnails_dir` corresponding to images in
    `images_dir`.
    """
    ensure_path(thumbnails_dir + os.path.sep)

    output('generating thumbnails...')
    filenames = glob.glob(os.path.join(images_dir, '*.png'))
    for fig_filename in filenames:
        ebase = fig_filename.replace(sfepy.data_dir, '').lstrip(os.path.sep)
        output('"%s"' % ebase)

        base = os.path.basename(fig_filename)
        thumb_filename = os.path.join(thumbnails_dir, base)

        image.thumbnail(fig_filename, thumb_filename, scale=scale)

    output('...done')

#.........这里部分代码省略.........
                        default=False, help=helps['3d'])
    parser.add_argument('--order', metavar='int', type=int,
                        action='store', dest='order',
                        default=1, help=helps['order'])
    options = parser.parse_args()

    dim = 3 if options.is_3d else 2
    dims = nm.array(eval(options.dims), dtype=nm.float64)[:dim]
    shape = nm.array(eval(options.shape), dtype=nm.int32)[:dim]
    centre = nm.array(eval(options.centre), dtype=nm.float64)[:dim]

    output('dimensions:', dims)
    output('shape:     ', shape)
    output('centre:    ', centre)

    mesh0 = gen_block_mesh(dims, shape, centre, name='block-fem',
                           verbose=True)
    domain0 = FEDomain('d', mesh0)

    bbox = domain0.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    eps = 1e-8 * (max_x - min_x)

    cnt = (shape[0] - 1) // 2
    g0 = 0.5 * dims[0]
    grading = nm.array([g0 / 2**ii for ii in range(cnt)]) + eps + centre[0] - g0

    domain, subs = refine_towards_facet(domain0, grading, 'x <')

    omega = domain.create_region('Omega', 'all')

    gamma1 = domain.create_region('Gamma1',
                                  'vertices in (x < %.10f)' % (min_x + eps),
                                  'facet')
    gamma2 = domain.create_region('Gamma2',
                                  'vertices in (x > %.10f)' % (max_x - eps),
                                  'facet')

    field = Field.from_args('fu', nm.float64, 1, omega,
                            approx_order=options.order)

    if subs is not None:
        field.substitute_dofs(subs)

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

    integral = Integral('i', order=2*options.order)

    t1 = Term.new('dw_laplace(v, u)',
                  integral, omega, v=v, u=u)
    eq = Equation('eq', t1)
    eqs = Equations([eq])

    def u_fun(ts, coors, bc=None, problem=None):
        """
        Define a displacement depending on the y coordinate.
        """
        if coors.shape[1] == 2:
            min_y, max_y = bbox[:, 1]
            y = (coors[:, 1] - min_y) / (max_y - min_y)

            val = (max_y - min_y) * nm.cos(3 * nm.pi * y)

        else:
            min_y, max_y = bbox[:, 1]
            min_z, max_z = bbox[:, 2]
            y = (coors[:, 1] - min_y) / (max_y - min_y)
            z = (coors[:, 2] - min_z) / (max_z - min_z)

            val = ((max_y - min_y) * (max_z - min_z)
                   * nm.cos(3 * nm.pi * y) * (1.0 + 3.0 * (z - 0.5)**2))

        return val

    bc_fun = Function('u_fun', u_fun)
    fix1 = EssentialBC('shift_u', gamma1, {'u.0' : bc_fun})
    fix2 = EssentialBC('fix2', gamma2, {'u.all' : 0.0})

    ls = ScipyDirect({})

    nls = Newton({}, lin_solver=ls)

    pb = Problem('heat', equations=eqs, nls=nls, ls=ls)

    pb.time_update(ebcs=Conditions([fix1, fix2]))
    state = pb.solve()

    if subs is not None:
        field.restore_dofs()

    filename = os.path.join(options.output_dir, 'hanging.vtk')
    ensure_path(filename)

    pb.save_state(filename, state)
    if options.order > 1:
        pb.save_state(filename, state, linearization=Struct(kind='adaptive',
                                                            min_level=0,
                                                            max_level=8,
                                                            eps=1e-3))

def generate_rst_files(rst_dir, examples_dir, images_dir):
    """
    Generate Sphinx rst files for examples in `examples_dir` with images
    in `images_dir` and put them into `rst_dir`.

    Returns
    -------
    dir_map : dict
        The directory mapping of examples and corresponding rst files.
    """
    ensure_path(rst_dir + os.path.sep)

    output('generating rst files...')

    dir_map = {}
    for ex_filename in locate_files('*.py', examples_dir):
        if _omit(ex_filename): continue

        ebase = ex_filename.replace(examples_dir, '')[1:]
        base_dir = os.path.dirname(ebase)

        rst_filename = os.path.basename(ex_filename).replace('.py', '.rst')

        dir_map.setdefault(base_dir, []).append((ex_filename, rst_filename))

    for dirname, filenames in dir_map.iteritems():
        filenames = sorted(filenames, cmp=lambda a, b: cmp(a[1], b[1]))
        dir_map[dirname ] = filenames

    # Main index.
    mfd = open(os.path.join(rst_dir, 'index.rst'), 'w')
    mfd.write(_index % ('sfepy', 'Examples', '=' * 8))

    for dirname, filenames in ordered_iteritems(dir_map):
        full_dirname = os.path.join(rst_dir, dirname)
        ensure_path(full_dirname + os.path.sep)

        # Subdirectory index.
        ifd = open(os.path.join(full_dirname, 'index.rst'), 'w')
        ifd.write(_index % (dirname, dirname, '=' * len(dirname)))

        for ex_filename, rst_filename in filenames:
            full_rst_filename = os.path.join(full_dirname, rst_filename)
            output('"%s"' % full_rst_filename.replace(rst_dir, '')[1:])
            rst_filename_ns = rst_filename.replace('.rst', '')
            ebase = ex_filename.replace(examples_dir, '')[1:]

            rst_ex_filename = _make_sphinx_path(ex_filename)
            docstring = get_default(import_file(ex_filename).__doc__,
                                    'missing description!')

            ifd.write('    %s\n' % rst_filename_ns)
            fig_include = ''
            fig_base = _get_fig_filenames(ebase, images_dir).next()
            for fig_filename in _get_fig_filenames(ebase, images_dir):
                rst_fig_filename = _make_sphinx_path(fig_filename)

                if os.path.exists(fig_filename):
                    fig_include += _image % rst_fig_filename + '\n'

            # Example rst file.
            fd = open(full_rst_filename, 'w')
            fd.write(_include % (fig_base, ebase, '=' * len(ebase),
                                 docstring,
                                 fig_include,
                                 rst_ex_filename, rst_ex_filename))
            fd.close()

        ifd.close()

        mfd.write('    %s/index\n' % dirname)

    mfd.close()

    output('...done')

    return dir_map

def generate_images(images_dir, examples_dir):
    """
    Generate images from results of running examples found in
    `examples_dir` directory.

    The generated images are stored to `images_dir`,
    """
    from sfepy.applications import solve_pde
    from sfepy.postprocess.viewer import Viewer
    from sfepy.postprocess.utils import mlab

    prefix = output.prefix

    output_dir = tempfile.mkdtemp()
    trunk = os.path.join(output_dir, 'result')
    options = Struct(output_filename_trunk=trunk,
                     output_format='vtk',
                     save_ebc=False,
                     save_ebc_nodes=False,
                     save_regions=False,
                     save_field_meshes=False,
                     save_regions_as_groups=False,
                     solve_not=False)
    default_views = {'' : {}}

    ensure_path(images_dir + os.path.sep)

    view = Viewer('', offscreen=False)

    for ex_filename in locate_files('*.py', examples_dir):
        if _omit(ex_filename): continue

        output.level = 0
        output.prefix = prefix
        ebase = ex_filename.replace(examples_dir, '')[1:]
        output('trying "%s"...' % ebase)

        try:
            problem, state = solve_pde(ex_filename, options=options)

        except KeyboardInterrupt:
            raise

        except:
            problem = None
            output('***** failed! *****')

        if problem is not None:
            if ebase in custom:
                views = custom[ebase]

            else:
                views = default_views

            tsolver = problem.get_time_solver()
            if tsolver.ts is None:
                suffix = None

            else:
                suffix = tsolver.ts.suffix % (tsolver.ts.n_step - 1)

            filename = problem.get_output_name(suffix=suffix)

            for suffix, kwargs in views.iteritems():
                fig_filename = _get_fig_filename(ebase, images_dir, suffix)

                fname = edit_filename(filename, suffix=suffix)
                output('displaying results from "%s"' % fname)
                disp_name = fig_filename.replace(sfepy.data_dir, '')
                output('to "%s"...' % disp_name.lstrip(os.path.sep))

                view.filename = fname
                view(scene=view.scene, show=False, is_scalar_bar=True, **kwargs)
                view.save_image(fig_filename)
                mlab.clf()

                output('...done')

            remove_files(output_dir)

        output('...done')

def main():
    parser = ArgumentParser(description=__doc__.rstrip(),
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('output_dir', help=helps['output_dir'])
    parser.add_argument('--dims', metavar='dims',
                        action='store', dest='dims',
                        default='1.0,1.0,1.0', help=helps['dims'])
    parser.add_argument('--shape', metavar='shape',
                        action='store', dest='shape',
                        default='11,11,11', help=helps['shape'])
    parser.add_argument('--centre', metavar='centre',
                        action='store', dest='centre',
                        default='0.0,0.0,0.0', help=helps['centre'])
    parser.add_argument('-2', '--2d',
                        action='store_true', dest='is_2d',
                        default=False, help=helps['2d'])
    parser.add_argument('--u-order', metavar='int', type=int,
                        action='store', dest='order_u',
                        default=1, help=helps['u-order'])
    parser.add_argument('--p-order', metavar='int', type=int,
                        action='store', dest='order_p',
                        default=1, help=helps['p-order'])
    parser.add_argument('--linearization', choices=['strip', 'adaptive'],
                        action='store', dest='linearization',
                        default='strip', help=helps['linearization'])
    parser.add_argument('--metis',
                        action='store_true', dest='metis',
                        default=False, help=helps['metis'])
    parser.add_argument('--silent',
                        action='store_true', dest='silent',
                        default=False, help=helps['silent'])
    parser.add_argument('--clear',
                        action='store_true', dest='clear',
                        default=False, help=helps['clear'])
    options, petsc_opts = parser.parse_known_args()

    comm = pl.PETSc.COMM_WORLD

    output_dir = options.output_dir

    filename = os.path.join(output_dir, 'output_log_%02d.txt' % comm.rank)
    if comm.rank == 0:
        ensure_path(filename)
    comm.barrier()

    output.prefix = 'sfepy_%02d:' % comm.rank
    output.set_output(filename=filename, combined=options.silent == False)

    output('petsc options:', petsc_opts)

    mesh_filename = os.path.join(options.output_dir, 'para.h5')

    if comm.rank == 0:
        from sfepy.mesh.mesh_generators import gen_block_mesh

        if options.clear:
            for _f in chain(*[glob.glob(os.path.join(output_dir, clean_pattern))
                              for clean_pattern
                              in ['*.h5', '*.txt', '*.png']]):
                output('removing "%s"' % _f)
                os.remove(_f)

        dim = 2 if options.is_2d else 3
        dims = nm.array(eval(options.dims), dtype=nm.float64)[:dim]
        shape = nm.array(eval(options.shape), dtype=nm.int32)[:dim]
        centre = nm.array(eval(options.centre), dtype=nm.float64)[:dim]

        output('dimensions:', dims)
        output('shape:     ', shape)
        output('centre:    ', centre)

        mesh = gen_block_mesh(dims, shape, centre, name='block-fem',
                              verbose=True)
        mesh.write(mesh_filename, io='auto')

    comm.barrier()

    output('field u order:', options.order_u)
    output('field p order:', options.order_p)

    solve_problem(mesh_filename, options, comm)

def run_test(conf_name, options, ifile):
    from sfepy.base.ioutils import ensure_path

    ensure_path(op.join(options.out_dir, "any"))

    if options.filter_none or options.raise_on_error:
        of = None
    elif options.filter_less:
        of = OutputFilter(["<<<", ">>>", "...", "!!!", "+++", "---"])
    elif options.filter_more:
        of = OutputFilter(["+++", "---"])
    else:
        of = OutputFilter(["<<<", "+++", "---"])

    print("<<< [%d] %s" % (ifile, conf_name))
    orig_prefix = output.get_output_prefix()
    output.set_output_prefix("[%d] %s" % (ifile, orig_prefix))

    _required, other = get_standard_keywords()
    required = ["Test"]

    num = 1
    test_time = 0.0
    try:
        conf = ProblemConf.from_file(conf_name, required, _required + other)
        test = conf.funmod.Test.from_conf(conf, options)
        num = test.get_number()
        ok = True
        print(">>> test instance prepared (%d test(s))" % num)
    except KeyboardInterrupt:
        print(">>> interrupted")
        sys.exit(0)
    except:
        print("--- test instance creation failed")
        if options.raise_on_error:
            raise
        ok, n_fail, n_total = False, num, num

    if ok:
        try:
            tt = time.clock()
            output.set_output_prefix(orig_prefix)
            ok, n_fail, n_total = test.run(debug=options.raise_on_error, ifile=ifile)
            output.set_output_prefix("[%d] %s" % (ifile, orig_prefix))
            test_time = time.clock() - tt
        except KeyboardInterrupt:
            print(">>> interrupted")
            sys.exit(0)
        except Exception as e:
            print(">>> %s" % e.__class__)
            if options.raise_on_error:
                raise
            ok, n_fail, n_total = False, num, num

    if ok:
        print(">>> all passed in %.2f s" % test_time)
    else:
        print("!!! %s test failed" % n_fail)

    if of is not None:
        of.stop()

    output.set_output_prefix(orig_prefix)

    return n_fail, n_total, test_time

def main():
    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-b', '--basis', metavar='name',
                      action='store', dest='basis',
                      default='lagrange', help=help['basis'])
    parser.add_option('-d', '--derivative', metavar='d', type=int,
                      action='store', dest='derivative',
                      default=0, help=help['derivative'])
    parser.add_option('-n', '--max-order', metavar='order', type=int,
                      action='store', dest='max_order',
                      default=2, help=help['max_order'])
    parser.add_option('-g', '--geometry', metavar='name',
                      action='store', dest='geometry',
                      default='2_4', help=help['geometry'])
    parser.add_option('-m', '--mesh', metavar='mesh',
                      action='store', dest='mesh',
                      default=None, help=help['mesh'])
    parser.add_option('', '--permutations', metavar='permutations',
                      action='store', dest='permutations',
                      default=None, help=help['permutations'])
    parser.add_option('', '--dofs', metavar='dofs',
                      action='store', dest='dofs',
                      default=None, help=help['dofs'])
    parser.add_option('-l', '--lin-options', metavar='options',
                      action='store', dest='lin_options',
                      default='min_level=2,max_level=5,eps=1e-3',
                      help=help['lin_options'])
    parser.add_option('', '--plot-dofs',
                      action='store_true', dest='plot_dofs',
                      default=False, help=help['plot_dofs'])
    options, args = parser.parse_args()

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

    output('polynomial space:', options.basis)
    output('max. order:', options.max_order)

    lin = Struct(kind='adaptive', min_level=2, max_level=5, eps=1e-3)
    for opt in options.lin_options.split(','):
        key, val = opt.split('=')
        setattr(lin, key, eval(val))

    if options.mesh is None:
        dim, n_ep = int(options.geometry[0]), int(options.geometry[2])
        output('reference element geometry:')
        output('  dimension: %d, vertices: %d' % (dim, n_ep))

        gel = GeometryElement(options.geometry)
        gps = PolySpace.any_from_args(None, gel, 1,
                                      base=options.basis)
        ps = PolySpace.any_from_args(None, gel, options.max_order,
                                     base=options.basis)

        n_digit, _format = get_print_info(ps.n_nod, fill='0')
        name_template = os.path.join(output_dir, 'bf_%s.vtk' % _format)
        for ip in get_dofs(options.dofs, ps.n_nod):
            output('shape function %d...' % ip)

            def eval_dofs(iels, rx):
                if options.derivative == 0:
                    bf = ps.eval_base(rx).squeeze()
                    rvals = bf[None, :, ip:ip+1]

                else:
                    bfg = ps.eval_base(rx, diff=True)
                    rvals = bfg[None, ..., ip]

                return rvals

            def eval_coors(iels, rx):
                bf = gps.eval_base(rx).squeeze()
                coors = nm.dot(bf, gel.coors)[None, ...]
                return coors

            (level, coors, conn,
             vdofs, mat_ids) = create_output(eval_dofs, eval_coors, 1,
                                             ps, min_level=lin.min_level,
                                             max_level=lin.max_level,
                                             eps=lin.eps)
            out = {
                'bf' : Struct(name='output_data',
                              mode='vertex', data=vdofs,
                              var_name='bf', dofs=None)
            }

            mesh = Mesh.from_data('bf_mesh', coors, None, [conn], [mat_ids],
                                  [options.geometry])

            name = name_template % ip
            ensure_path(name)
            mesh.write(name, out=out)

            output('...done (%s)' % name)

    else:
        mesh = Mesh.from_file(options.mesh)
#.........这里部分代码省略.........

def gen_misc_mesh(mesh_dir, force_create, kind, args, suffix='.mesh',
                  verbose=False):
    """
    Create sphere or cube mesh according to `kind` in the given
    directory if it does not exist and return path to it.
    """
    import os
    from sfepy import data_dir

    defdir = os.path.join(data_dir, 'meshes')
    if mesh_dir is None:
        mesh_dir = defdir

    def retype(args, types, defaults):
        args=list(args)
        args.extend(defaults[len(args):len(defaults)])
        return tuple([type(value) for type, value in zip(types, args) ])

    if kind == 'sphere':
        default = [5, 41, args[0]]
        args = retype(args, [float, int, float], default)
        mesh_pattern = os.path.join(mesh_dir, 'sphere-%.2f-%.2f-%i')

    else:
        assert_(kind == 'cube')

        args = retype(args,
                      (int, float, int, float, int, float),
                      (args[0], args[1], args[0], args[1], args[0], args[1]))
        mesh_pattern = os.path.join(mesh_dir, 'cube-%i_%.2f-%i_%.2f-%i_%.2f')

    if verbose:
        output(args)

    filename = mesh_pattern % args
    if not force_create:
        if os.path.exists(filename): return filename
        if os.path.exists(filename + '.mesh') : return filename + '.mesh'
        if os.path.exists(filename + '.vtk'): return filename + '.vtk'

    if kind == 'cube':
        filename = filename + suffix
        ensure_path(filename)

        output('creating new cube mesh')
        output('(%i nodes in %.2f) x (%i nodes in %.2f) x (%i nodes in %.2f)'
               % args)
        output('to file %s...' % filename)

        mesh = gen_block_mesh(args[1::2], args[0::2],
                              (0.0, 0.0, 0.0), name=filename)
        mesh.write(filename, io='auto')
        output('...done')

    else:
        import subprocess, shutil, tempfile
        filename = filename + '.mesh'
        ensure_path(filename)

        output('creating new sphere mesh (%i nodes, r=%.2f) and gradation %d'
               % args)
        output('to file %s...' % filename)

        f = open(os.path.join(defdir, 'quantum', 'sphere.geo'))
        tmp_dir = tempfile.mkdtemp()
        tmpfile = os.path.join(tmp_dir, 'sphere.geo.temp')
        ff = open(tmpfile, "w")
        ff.write("""
R = %i.0;
n = %i.0;
dens = %f;
""" % args)
        ff.write(f.read())
        f.close()
        ff.close()
        subprocess.call(['gmsh', '-3', tmpfile, '-format', 'mesh',
                         '-o', filename])
        shutil.rmtree(tmp_dir)
        output('...done')

    return filename

def main():
    parser = OptionParser(usage=usage, version='%prog ' + sfepy.__version__)
    parser.add_option('-c', '--conf', metavar='"key : value, ..."',
                      action='store', dest='conf', type='string',
                      default=None, help= help['conf'])
    parser.add_option('-O', '--options', metavar='"key : value, ..."',
                      action='store', dest='app_options', type='string',
                      default=None, help=help['options'])
    parser.add_option('-o', '', metavar='filename',
                      action='store', dest='output_filename_trunk',
                      default=None, help=help['filename'])
    parser.add_option('--create-mesh',
                      action='store_true', dest='create_mesh',
                      default=False, help=help['create_mesh'])
    parser.add_option('--2d',
                      action='store_true', dest='dim2',
                      default=False, help=help['dim'])
    parser.add_option('-m', '--mesh', metavar='filename',
                      action='store', dest='mesh',
                      default=None, help=help['mesh'])
    parser.add_option('--mesh-dir', metavar='dirname',
                      action='store', dest='mesh_dir',
                      default='tmp', help=help['mesh_dir'])
    parser.add_option('--oscillator',
                      action='store_true', dest='oscillator',
                      default=False, help=help['oscillator'])
    parser.add_option('--well',
                      action='store_true', dest='well',
                      default=False, help=help['well'])
    parser.add_option('--hydrogen',
                      action='store_true', dest='hydrogen',
                      default=False, help=help['hydrogen'])
    parser.add_option('--boron',
                      action='store_true', dest='boron',
                      default=False, help=help['boron'])

    options, args = parser.parse_args()

    if options.create_mesh and options.mesh:
        output('--create-mesh and --mesh options are mutually exclusive!')
        return

    if len(args) == 1:
        filename_in = args[0];
        auto_mesh_name = False

    elif len(args) == 0:
        auto_mesh_name = True

        mesh_filename = os.path.join(options.mesh_dir, 'mesh.vtk')
        ensure_path(mesh_filename)

        if options.oscillator:
            filename_in = fix_path("examples/quantum/oscillator.py")

        elif options.well:
            filename_in = fix_path("examples/quantum/well.py")

        elif options.hydrogen:
            filename_in = fix_path("examples/quantum/hydrogen.py")

        elif options.boron:
            filename_in = fix_path("examples/quantum/boron.py")

        elif options.create_mesh:
            output('generating mesh...')
            try:
                os.makedirs("tmp")
            except OSError, e:
                if e.errno != 17: # [Errno 17] File exists
                    raise
            if options.dim2:
                output("dimension: 2")
                gp = fix_path('meshes/quantum/square.geo')
                os.system("cp %s tmp/mesh.geo" % gp)
                os.system("gmsh -2 tmp/mesh.geo -format mesh")
                mtv = fix_path('script/mesh_to_vtk.py')
                os.system("%s tmp/mesh.mesh %s" % (mtv, mesh_filename))
            else:
                output("dimension: 3")
                import sfepy.geom as geom
                from sfepy.fem.mesh import Mesh
                try:
                    from site_cfg import tetgen_path
                except ImportError:
                    tetgen_path = '/usr/bin/tetgen'
                gp = fix_path('meshes/quantum/box.geo')
                os.system("gmsh -0 %s -o tmp/x.geo" % gp)
                g = geom.read_gmsh("tmp/x.geo")
                g.printinfo()
                geom.write_tetgen(g, "tmp/t.poly")
                geom.runtetgen("tmp/t.poly", a=0.03, Q=1.0,
                               quadratic=False, tetgenpath=tetgen_path)
                m = Mesh.from_file("tmp/t.1.node")
                m.write(mesh_filename, io="auto")
            output("...mesh written to %s" % mesh_filename)
            return

        else:
            parser.print_help()
#.........这里部分代码省略.........

    def __call__(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                 i_max=None, mtx=None, status=None, **kwargs):
        import os, sys, shutil, tempfile
        from sfepy import base_dir
        from sfepy.base.ioutils import ensure_path

        eps_a = get_default(eps_a, self.conf.eps_a)
        eps_r = get_default(eps_r, self.conf.eps_r)
        i_max = get_default(i_max, self.conf.i_max)
        eps_d = self.conf.eps_d

        petsc = self.petsc

        # There is no use in caching matrix in the solver - always set as new.
        pmtx, psol, prhs = self.set_matrix(mtx)

        ksp = self.ksp
        ksp.setOperators(pmtx)
        ksp.setFromOptions() # PETSc.Options() not used yet...
        ksp.setTolerances(atol=eps_a, rtol=eps_r, divtol=eps_d, max_it=i_max)

        output_dir = tempfile.mkdtemp()

        # Set PETSc rhs, solve, get solution from PETSc solution.
        if x0 is not None:
            psol[...] = x0
            sol0_filename = os.path.join(output_dir, 'sol0.dat')

        else:
            sol0_filename = ''

        prhs[...] = rhs

        script_filename = os.path.join(base_dir, 'solvers/petsc_worker.py')

        mtx_filename = os.path.join(output_dir, 'mtx.dat')
        rhs_filename = os.path.join(output_dir, 'rhs.dat')
        sol_filename = os.path.join(output_dir, 'sol.dat')
        status_filename = os.path.join(output_dir, 'status.txt')

        log_filename = os.path.join(self.conf.log_dir, 'sol.log')
        ensure_path(log_filename)

        output('storing system to %s...' % output_dir)
        tt = time.clock()
        view_mtx = petsc.Viewer().createBinary(mtx_filename, mode='w')
        view_rhs = petsc.Viewer().createBinary(rhs_filename, mode='w')
        pmtx.view(view_mtx)
        prhs.view(view_rhs)
        if sol0_filename:
            view_sol0 = petsc.Viewer().createBinary(sol0_filename, mode='w')
            psol.view(view_sol0)
        output('...done in %.2f s' % (time.clock() - tt))

        command = [
            'mpiexec -n %d' % self.conf.n_proc,
            sys.executable, script_filename,
            '-mtx %s' % mtx_filename, '-rhs %s' % rhs_filename,
            '-sol0 %s' % sol0_filename, '-sol %s' % sol_filename,
            '-status %s' % status_filename,
            '-ksp_type %s' % self.conf.method,
            '-pc_type %s' % self.conf.precond,
            '-sub_pc_type %s' % self.conf.sub_precond,
            '-ksp_atol %.3e' % self.conf.eps_a,
            '-ksp_rtol %.3e' % self.conf.eps_r,
            '-ksp_max_it %d' % self.conf.i_max,
            '-ksp_monitor %s' % log_filename,
            '-ksp_view %s' % log_filename,
        ]
        if self.conf.precond_side is not None:
            command.append('-ksp_pc_side %s' % self.conf.precond_side)

        out = os.system(" ".join(command))
        assert_(out == 0)

        output('reading solution...')
        tt = time.clock()
        view_sol = self.petsc.Viewer().createBinary(sol_filename, mode='r')
        psol = petsc.Vec().load(view_sol)

        fd = open(status_filename, 'r')
        line = fd.readline().split()
        reason = int(line[0])
        elapsed = float(line[1])
        fd.close()
        output('...done in %.2f s' % (time.clock() - tt))

        sol = psol[...].copy()
        output('%s(%s, %s/proc) convergence: %s (%s)'
               % (self.conf.method, self.conf.precond, self.conf.sub_precond,
                  reason, self.converged_reasons[reason]))
        output('elapsed: %.2f [s]' % elapsed)

        shutil.rmtree(output_dir)

        return sol

def generate_rst_files(rst_dir, examples_dir, images_dir):
    """
    Generate Sphinx rst files for examples in `examples_dir` with images
    in `images_dir` and put them into `rst_dir`.

    Returns
    -------
    dir_map : dict
        The directory mapping of examples and corresponding rst files.
    """
    ensure_path(rst_dir + os.path.sep)

    output("generating rst files...")

    dir_map = {}
    for ex_filename in locate_files("*.py", examples_dir):
        if _omit(ex_filename):
            continue

        ebase = ex_filename.replace(examples_dir, "")[1:]
        base_dir = os.path.dirname(ebase)

        rst_filename = os.path.basename(ex_filename).replace(".py", ".rst")

        dir_map.setdefault(base_dir, []).append((ex_filename, rst_filename))

    for dirname, filenames in dir_map.iteritems():
        filenames = sorted(filenames, cmp=lambda a, b: cmp(a[1], b[1]))
        dir_map[dirname] = filenames

    # Main index.
    mfd = open(os.path.join(rst_dir, "index.rst"), "w")
    mfd.write(_index % ("sfepy", "SfePy autogenerated gallery", "=" * 27))

    for dirname, filenames in ordered_iteritems(dir_map):
        full_dirname = os.path.join(rst_dir, dirname)
        ensure_path(full_dirname + os.path.sep)

        # Subdirectory index.
        ifd = open(os.path.join(full_dirname, "index.rst"), "w")
        ifd.write(_index % (dirname, dirname, "=" * len(dirname)))

        for ex_filename, rst_filename in filenames:
            full_rst_filename = os.path.join(full_dirname, rst_filename)
            output('"%s"' % full_rst_filename.replace(rst_dir, "")[1:])
            rst_filename_ns = rst_filename.replace(".rst", "")
            ebase = ex_filename.replace(examples_dir, "")[1:]

            rst_ex_filename = _make_sphinx_path(ex_filename)
            docstring = get_default(import_file(ex_filename).__doc__, "missing description!")

            ifd.write("    %s\n" % rst_filename_ns)
            fig_include = ""
            fig_base = _get_fig_filenames(ebase, images_dir).next()
            for fig_filename in _get_fig_filenames(ebase, images_dir):
                rst_fig_filename = _make_sphinx_path(fig_filename)

                if os.path.exists(fig_filename):
                    fig_include += _image % rst_fig_filename + "\n"

            # Example rst file.
            fd = open(full_rst_filename, "w")
            fd.write(
                _include % (fig_base, ebase, "=" * len(ebase), docstring, fig_include, rst_ex_filename, rst_ex_filename)
            )
            fd.close()

        ifd.close()

        mfd.write("    %s/index\n" % dirname)

    mfd.close()

    output("...done")

    return dir_map

#.........这里部分代码省略.........

    if options.conf is not None:
        mod = import_file(options.conf)

    else:
        mod = sys.modules[__name__]

    apply_units = mod.apply_units
    define = mod.define
    set_wave_dir = mod.set_wave_dir
    setup_n_eigs = mod.setup_n_eigs
    build_evp_matrices = mod.build_evp_matrices
    save_materials = mod.save_materials
    get_std_wave_fun = mod.get_std_wave_fun
    get_stepper = mod.get_stepper
    process_evp_results = mod.process_evp_results

    options.pars = [float(ii) for ii in options.pars.split(',')]
    options.unit_multipliers = [float(ii)
                                for ii in options.unit_multipliers.split(',')]
    options.wave_dir = [float(ii)
                        for ii in options.wave_dir.split(',')]
    aux = options.range.split(',')
    options.range = [float(aux[0]), float(aux[1]), int(aux[2])]
    options.solver_conf = dict_from_string(options.solver_conf)

    if options.clear:
        remove_files_patterns(output_dir,
                              ['*.h5', '*.vtk', '*.txt'],
                              ignores=['output_log.txt'],
                              verbose=True)

    filename = os.path.join(output_dir, 'options.txt')
    ensure_path(filename)
    save_options(filename, [('options', vars(options))],
                 quote_command_line=True)

    pars = apply_units(options.pars, options.unit_multipliers)
    output('material parameters with applied unit multipliers:')
    output(pars)

    if options.mode == 'omega':
        rng = copy(options.range)
        rng[:2] = apply_unit_multipliers(options.range[:2],
                                         ['wave_number', 'wave_number'],
                                         options.unit_multipliers)
        output('wave number range with applied unit multipliers:', rng)

    else:
        rng = copy(options.range)
        rng[:2] = apply_unit_multipliers(options.range[:2],
                                         ['frequency', 'frequency'],
                                         options.unit_multipliers)
        output('frequency range with applied unit multipliers:', rng)

    pb, wdir, bzone, mtxs = assemble_matrices(define, mod, pars, set_wave_dir,
                                              options)
    dim = pb.domain.shape.dim

    if dim != 2:
        options.plane = 'strain'

    if options.save_regions:
        pb.save_regions_as_groups(os.path.join(output_dir, 'regions'))

    if options.save_materials:

def save_basis_on_mesh(mesh, options, output_dir, lin,
                       permutations=None, suffix=''):
    domain = Domain('domain', mesh)

    if permutations is not None:
        for group in domain.iter_groups():
            perms = group.gel.get_conn_permutations()[permutations]
            offsets = nm.arange(group.shape.n_el) * group.shape.n_ep

            group.conn[:] = group.conn.take(perms + offsets[:, None])

        domain.setup_facets()

    omega = domain.create_region('Omega', 'all')
    field = Field.from_args('f', nm.float64, shape=1, region=omega,
                            approx_order=options.max_order,
                            poly_space_base=options.basis)
    var = FieldVariable('u', 'unknown', field, 1)

    if options.plot_dofs:
        import sfepy.postprocess.plot_dofs as pd
        group = domain.groups[0]
        ax = pd.plot_mesh(None, mesh.coors, mesh.conns[0], group.gel.edges)
        ax = pd.plot_global_dofs(ax, field.get_coor(), field.aps[0].econn)
        ax = pd.plot_local_dofs(ax, field.get_coor(), field.aps[0].econn)
        if options.dofs is not None:
            ax = pd.plot_nodes(ax, field.get_coor(), field.aps[0].econn,
                               field.aps[0].interp.poly_spaces['v'].nodes,
                               get_dofs(options.dofs, var.n_dof))
        pd.plt.show()

    output('dofs: %d' % var.n_dof)

    vec = nm.empty(var.n_dof, dtype=var.dtype)
    n_digit, _format = get_print_info(var.n_dof, fill='0')
    name_template = os.path.join(output_dir,
                                 'dof_%s_%s.vtk' % (_format, suffix))
    for ip in get_dofs(options.dofs, var.n_dof):
        output('dof %d...' % ip)

        vec.fill(0.0)
        vec[ip] = 1.0

        var.set_data(vec)

        if options.derivative == 0:
            out = var.create_output(vec, linearization=lin)

        else:
            out = create_expression_output('ev_grad.ie.Elements(u)',
                                           'u', 'f', {'f' : field}, None,
                                           Variables([var]),
                                           mode='qp', verbose=False,
                                           min_level=lin.min_level,
                                           max_level=lin.max_level,
                                           eps=lin.eps)

        name = name_template % ip
        ensure_path(name)
        out['u'].mesh.write(name, out=out)

        output('...done (%s)' % name)

def generate_images(images_dir, examples_dir):
    """
    Generate images from results of running examples found in
    `examples_dir` directory.

    The generated images are stored to `images_dir`,
    """
    from sfepy.applications import solve_pde
    from sfepy.postprocess import Viewer
    from sfepy.postprocess.utils import mlab

    prefix = output.prefix

    output_dir = tempfile.mkdtemp()
    trunk = os.path.join(output_dir, "result")
    options = Struct(
        output_filename_trunk=trunk,
        output_format="vtk",
        save_ebc=False,
        save_regions=False,
        save_field_meshes=False,
        save_regions_as_groups=False,
        solve_not=False,
    )

    ensure_path(images_dir + os.path.sep)

    view = Viewer("", output_dir=output_dir, offscreen=False)

    for ex_filename in locate_files("*.py", examples_dir):
        if _omit(ex_filename):
            continue

        output.level = 0
        output.prefix = prefix
        ebase = ex_filename.replace(examples_dir, "")[1:]
        output('trying "%s"...' % ebase)

        try:
            problem, state = solve_pde(ex_filename, options=options)

        except KeyboardInterrupt:
            raise

        except:
            problem = None
            output("***** failed! *****")

        if problem is not None:
            fig_filename = _get_fig_filename(ebase, images_dir)[1]

            if problem.ts_conf is None:
                filename = trunk + ".vtk"

            else:
                suffix = problem.ts.suffix % problem.ts.step
                filename = problem.get_output_name(suffix=suffix)

            output('displaying results from "%s"' % filename)
            output('to "%s"...' % fig_filename.replace(sfepy.data_dir, "")[1:])

            view.filename = filename
            view(scene=view.scene, show=False, is_scalar_bar=True, fig_filename=fig_filename)
            mlab.clf()

            output("...done")

            remove_files(output_dir)

        output("...done")

def main():
    parser = ArgumentParser(description=__doc__.rstrip(),
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('output_dir', help=helps['output_dir'])
    parser.add_argument('-d', '--dims', metavar='l,w,t',
                        action='store', dest='dims',
                        default='0.2,0.01,0.001', help=helps['dims'])
    parser.add_argument('-n', '--nx', metavar='start,stop,step',
                        action='store', dest='nx',
                        default='2,103,10', help=helps['nx'])
    parser.add_argument('-t', '--transform', choices=['none', 'bend', 'twist'],
                        action='store', dest='transform',
                        default='none', help=helps['transform'])
    parser.add_argument('--young', metavar='float', type=float,
                        action='store', dest='young',
                        default=210e9, help=helps['young'])
    parser.add_argument('--poisson', metavar='float', type=float,
                        action='store', dest='poisson',
                        default=0.3, help=helps['poisson'])
    parser.add_argument('--force', metavar='float', type=float,
                        action='store', dest='force',
                        default=-1.0, help=helps['force'])
    parser.add_argument('-p', '--plot',
                        ac