def _save_dict_latex(self, adict, fd, names):
        fd.write( r'\begin{itemize}' )
        fd.write( '\n' )
        for key, val in ordered_iteritems(adict):
            if key.startswith('_a_'): continue

            try:
                lname = names[key]
            except:
                lname = self._escape_latex(key)
            fd.write( '\item %s:' % lname )
            fd.write( '\n' )

            if isinstance(val, dict):
                self._save_dict_latex(val, fd, names)

            elif isinstance(val, basestr):
                fd.write(self._escape_latex(val) + '\n')

            elif val.ndim == 0:
                fd.write('$' + self._format(val) + '$\n')

            elif val.ndim == 1:
                self._write1d(fd, val)

            elif val.ndim == 2:
                self._write2d(fd, val)

        fd.write( r'\end{itemize}' )
        fd.write( '\n\n' )

    def test_base_functions_delta(self):
        """
        Test :math:`\delta` property of base functions evaluated in the
        reference element nodes.
        """
        from sfepy.base.base import ordered_iteritems
        from sfepy.discrete.fem.poly_spaces import PolySpace

        ok = True

        for key, gel in ordered_iteritems(self.gels):
            for order in range(11):
                ps = PolySpace.any_from_args('aux', gel, order,
                                             base='lagrange',
                                             force_bubble=False)
                bf = ps.eval_base(ps.node_coors)
                _ok = nm.allclose(nm.eye(ps.n_nod),
                                  bf.squeeze(),
                                  rtol=0.0, atol=(order + 1) * 1e-14)

                self.report('%s order %d (n_nod: %d): %s'
                            % (key, order, ps.n_nod, _ok))

                if not _ok:
                    import pdb; pdb.set_trace()

            ok = ok and _ok

        return ok

        def poststep_fun(ts, vec):
            _poststep_fun(ts, vec)

            # Probe the solution.
            dvel_qp = ev('ev_diffusion_velocity.%d.Omega(m.diffusivity, T)'
                         % order, copy_materials=False, mode='qp')
            project_by_component(dvel, dvel_qp, component, order,
                                 nls_options=nls_options)

            all_results = []
            for ii, probe in enumerate(probes):
                fig, results = probe_results(ii, T, dvel, probe, labels[ii])

                all_results.append(results)

            plt.tight_layout()
            fig.savefig('time_poisson_interactive_probe_%s.png'
                        % (suffix % ts.step), bbox_inches='tight')

            for ii, results in enumerate(all_results):
                output('probe %d (%s):' % (ii, probes[ii].name))
                output.level += 2
                for key, res in ordered_iteritems(results):
                    output(key + ':')
                    val = res[1]
                    output('  min: %+.2e, mean: %+.2e, max: %+.2e'
                           % (val.min(), val.mean(), val.max()))
                output.level -= 2

    def test_base_functions_values(self):
        """
        Compare base function values and their gradients with correct
        data. Also test that sum of values over all element nodes gives one.
        """
        from sfepy.base.base import ordered_iteritems
        from sfepy.discrete.fem.poly_spaces import PolySpace

        ok = True

        for key, val in ordered_iteritems(test_bases):
            gel = self.gels[key[:3]]
            diff = key[-4:] == "grad"
            order = int(key[5])
            force_bubble = key[6:7] == "B"

            ps = PolySpace.any_from_args("aux", gel, order, base="lagrange", force_bubble=force_bubble)
            dim = ps.geometry.dim
            coors = nm.r_[ps.geometry.coors, [[0.2] * dim]]

            bf = ps.eval_base(coors, diff=diff)
            _ok = nm.allclose(val, bf, rtol=0.0, atol=1e-14)
            ## if not _ok:
            ##     nm.set_printoptions(threshold=1000000, linewidth=65)
            ##     print bf.__repr__()

            if not diff:
                _ok = _ok and nm.allclose(bf.sum(axis=2), 1.0, rtol=0.0, atol=1e-14)

            self.report("%s: %s" % (key, _ok))

            ok = ok and _ok

        return ok

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 plot_data(self, igs):
        send = self.plot_pipe.send

        ii = 0
        for ig, names in ordered_iteritems(self.data_names):
            if ig in igs:
                send(['ig', ig])
                send(['clear'])
                for name in names:
                    key = name_to_key(name, ii)
                    try:
                        send(['plot',
                              nm.array(self.x_values[ig]),
                              nm.array(self.data[key])])
                    except:
                        msg = "send failed! (%s, %s, %s)!" \
                              % (ii, name, self.data[key])
                        raise IOError(msg)
                    ii += 1

            else:
                ii += len(names)

        send(['legends'])
        send(['continue'])

    def test_weight_consistency(self):
        """
        Test if integral of 1 (= sum of weights) gives the domain volume.
        """
        from sfepy.fem.quadratures import quadrature_tables

        ok = True
        for geometry, qps in ordered_iteritems(quadrature_tables):
            self.report('geometry:', geometry)
            for order, qp in ordered_iteritems(qps):
                _ok = nm.allclose(qp.weights.sum(), qp.volume,
                                  rtol=0.0, atol=1e-15)
                self.report('order %d: %s' % (order, _ok))

                ok = ok and _ok

        return ok

    def test_weight_consistency(self):
        """
        Test if integral of 1 (= sum of weights) gives the domain volume.
        """
        from sfepy.discrete.quadratures import quadrature_tables

        ok = True
        for geometry, qps in ordered_iteritems(quadrature_tables):
            self.report("geometry:", geometry)
            for order, qp in ordered_iteritems(qps):
                diff = nm.abs(qp.weights.sum() - qp.volume)
                _ok = diff < 1e-14
                self.report("order %d: %s (%.2e)" % (order, _ok, diff))

                ok = ok and _ok

        return ok

def plot_partitioning(axs, field, cell_tasks, gfd, output_dir, size):
    """
    Plot the partitioning of the domain and field DOFs.
    """
    mesh = field.domain.mesh

    ax = pc.plot_wireframe(axs[0], mesh.cmesh)
    coors = field.get_coor()
    econn = field.econn
    ax = pd.plot_global_dofs(ax, coors, econn)
    ax.set_title('global DOFs')
    ax.figure.savefig(os.path.join(output_dir, 'global_dofs.png'),
                      bbox_inches='tight')

    ax = pc.plot_wireframe(axs[1], mesh.cmesh)
    fig = ax.figure
    coors = field.get_coor()
    econn = field.econn

    id_map = gfd.id_map

    colors = nm.zeros((field.n_nod, 4))
    for ir, dof_map in ordered_iteritems(gfd.dof_maps):
        aux = id_map[dof_map[0]]
        colors[aux] = [0, 0, float(ir + 1) / size, 0.6]
        for aux in dof_map[1]:
            colors[id_map[aux]] = [0, 0, float(ir + 1) / size, 0.9]

    from sfepy.discrete.fem.utils import prepare_translate
    aux = prepare_translate(id_map[econn], econn)
    ax = label_dofs(ax, coors[aux], id_map, colors)

    mark_subdomains(ax, mesh.cmesh, cell_tasks, size, 0, 0.7)

    ax.set_title('subdomains of tasks and PETSc DOFs')

    fig.savefig(os.path.join(output_dir, 'petsc_dofs.png'),
                bbox_inches='tight')

    ax.set_title('')

    axis = ax.axis()
    for ir, ocells in enumerate(gfd.overlap_cells):
        aux = nm.zeros_like(cell_tasks)
        aux[ocells] = 10
        aux[gfd.cell_parts[ir]] = 1

        ax = fig.add_axes(ax.get_position(), frameon=False, label='aux')
        mark_subdomains(ax, mesh.cmesh, aux, 11, 1, 0.3, True)
        ax.axis(axis)
        ax.set_title('overlap cells on task %d' % ir)
        fig.savefig(os.path.join(output_dir,
                                 'petsc_overlaps_%02d.png' % ir),
                    bbox_inches='tight')
        fig.delaxes(ax)

def generate_gallery_html(output_filename, gallery_dir,
                          rst_dir, thumbnails_dir, dir_map, link_prefix=''):
    """
    Generate the gallery html file with thumbnail images and links to
    examples.

    Parameters
    ----------
    output_filename : str
        The output html file name.
    gallery_dir : str
        The top level directory of gallery files.
    rst_dir : str
        The full path to rst files of examples within `gallery_dir`.
    thumbnails_dir : str
        The full path to thumbnail images within `gallery_dir`.
    dir_map : dict
        The directory mapping returned by `generate_rst_files()`
    link_prefix : str, optional
        The prefix to prepend to links to individual pages of examples.
    """
    output('generating %s...' % output_filename)

    lines = []
    for dirname, filenames in ordered_iteritems(dir_map):
        full_dirname = os.path.join(rst_dir, dirname)

        for ex_filename, rst_filename in filenames:
            full_rst_filename = os.path.join(full_dirname, rst_filename)

            ebase = full_rst_filename.replace(rst_dir, '')[1:]
            ebase = edit_filename(ebase, new_ext='.py')

            link_base = full_rst_filename.replace(gallery_dir, '')[1:]
            link = os.path.join(link_prefix,
                                os.path.splitext(link_base)[0] + '.html')

            _get_fig_filenames(ebase, thumbnails_dir).next()
            for thumbnail_filename in _get_fig_filenames(ebase,
                                                         thumbnails_dir):
                if not os.path.isfile(thumbnail_filename):
                    # Skip examples with no image (= failed examples).
                    continue

                thumbnail_name = thumbnail_filename.replace(gallery_dir,
                                                            '')[1:]
                line = _link_template % (link, thumbnail_name,
                                         os.path.splitext(ebase)[0])
                lines.append(line)

    fd = open(output_filename, 'w')
    fd.write(_gallery_template % (link_prefix, '\n'.join(lines)))
    fd.close()

    output('...done')

    def iter_names(self, igs=None):
        if igs is None:
            igs = nm.arange(self.n_gr)

        ii = iseq = 0
        for ig, names in ordered_iteritems(self.data_names):
            for name in names:
                if ig in igs:
                    yield ig, ii, iseq, name
                    iseq += 1
                ii += 1

def save_raw_bg_logs(filename, logs):
    """
    Save raw band gaps `logs` into the `filename` file.
    """
    out = {}

    iranges = nm.cumsum([0] + [len(ii) for ii in logs.freqs])
    out['iranges'] = iranges
    for key, log in ordered_iteritems(logs.to_dict()):
        out[key] = nm.concatenate(log, axis=0)

    with open(filename, 'w') as fd:
        nm.savez(fd, **out)

    def test_quadratures(self):
        """
        Test if the quadratures have orders they claim to have, using
        symbolic integration by sympy.
        """
        from sfepy.discrete.quadratures import quadrature_tables, tp_geometries
        from sfepy.discrete.integrals import Integral

        if sm is None:
            self.report("cannot import sympy, skipping")
            return True

        quad = Integral("test_integral")

        ok = True
        failed = []
        for geometry, qps in ordered_iteritems(quadrature_tables):
            self.report("geometry:", geometry)

            if geometry in tp_geometries:
                iter_qp = range(1, 11)

            elif geometry == "2_3":
                iter_qp = range(1, 25)

            elif geometry == "3_4":
                iter_qp = range(1, 12)

            else:
                iter_qp = sorted(qps.keys())

            for order in iter_qp:
                self.report("order:", order)

                aux = self._test_quadratures(quad, geometry, order)
                _ok, integral, val = aux

                if not _ok:
                    failed.append((geometry, order, float(integral), val))

                ok = ok and _ok

        if not ok:
            self.report("failed:")
            for aux in failed:
                self.report(aux, "%.1e" % abs(aux[2] - aux[3]))

        return ok

    def test_quadratures(self):
        """
        Test if the quadratures have orders they claim to have, using
        symbolic integration by sympy.
        """
        from sfepy.fem.quadratures import quadrature_tables, tp_geometries
        from sfepy.fem.integrals import Integral

        if sm is None:
            self.report('cannot import sympy, skipping')
            return True

        quad = Integral('test_integral')

        ok = True
        failed = []
        for geometry, qps in ordered_iteritems(quadrature_tables):
            self.report('geometry:', geometry)

            if geometry in tp_geometries:
                iter_qp = xrange(1, 11)

            elif geometry == '2_3':
                iter_qp = xrange(1, 25)

            elif geometry == '3_4':
                iter_qp = xrange(1, 12)

            else:
                iter_qp = sorted(qps.keys())

            for order in iter_qp:
                self.report('order:', order)

                aux = self._test_quadratures(quad, geometry, order)
                _ok, integral, val = aux

                if not _ok:
                    failed.append((geometry, order, float(integral), val))

                ok = ok and _ok

        if not ok:
            self.report('failed:')
            for aux in failed:
                self.report(aux, '%.1e' % abs(aux[2] - aux[3]))

        return ok

def write_log(output, log, info):
    xlabels, ylabels, yscales, names, plot_kwargs = zip(*info.values())
    _write_header(output, xlabels, ylabels, yscales, names, plot_kwargs)

    for ii, (xlabel, ylabel, yscale, names, plot_kwargs) \
        in ordered_iteritems(info):
        for ip, name in enumerate(names):
            xs, ys, vlines = log[name]

            for ir, x in enumerate(xs):
                output('{}: {}: {:.16e}'.format(name, x, ys[ir]))

                if x in vlines:
                    output(name + ': -----')

    output('# ended: %s' % time.asctime())

 def _save_dict( self, adict, fd, names, format ):
     for key, val in ordered_iteritems(adict):
         try:
             lname = names[key]
         except:
             lname = key
         fd.write( '%s:\n' % lname )
         if isinstance( val, dict ):
             self._save_dict( val, fd, names, format )
             fd.write( '\n' )
         elif isinstance(val, basestr):
             fd.write( val + '\n' )
         elif isinstance( val, float ):
             fd.write( '%e\n' % val )
         elif val.ndim == 0:
             fd.write( format % val )
             fd.write( '\n' )
         elif val.ndim == 1:
             for ic in xrange( val.shape[0] ):
                 fd.write( format % val[ic] )
                 if ic < (val.shape[0] - 1):
                     fd.write( ', ' )
                 else:
                     fd.write( '\n' )
         elif val.ndim == 2:
             for ir in xrange( val.shape[0] ):
                 for ic in xrange( val.shape[1] ):
                     fd.write( format % val[ir,ic] )
                     if ic < (val.shape[1] - 1):
                         fd.write( ', ' )
                     elif ir < (val.shape[0] - 1):
                         fd.write( ';\n' )
             fd.write( '\n' )
         elif val.ndim == 3:
             for ii in xrange( val.shape[0] ):
                 fd.write( '  step %d:\n' % ii )
                 for ir in xrange( val.shape[1] ):
                     for ic in xrange( val.shape[2] ):
                         fd.write( '  ' + format % val[ii,ir,ic] )
                         if ic < (val.shape[2] - 1):
                             fd.write( ', ' )
                         elif ir < (val.shape[1] - 1):
                             fd.write( ';\n' )
                 fd.write( '\n' )
             fd.write( '\n' )
         fd.write( '\n' )

    def _save_dict_latex(self, adict, fd, names, idx=None):
        fd.write( r'\begin{itemize}' )
        fd.write( '\n' )
        for key, val in ordered_iteritems(adict):
            if key.startswith('_a_'): continue

            try:
                lname = names[key]
            except:
                lname = self._escape_latex(key)
            fd.write( '\item %s:' % lname )
            fd.write( '\n' )

            if isinstance(val, list):
                if idx is not None:
                    val = val[idx]
                else:
                    raise NotImplementedError("'idx' must be set in the case "
                                              "of multi-coefficients!")

            if isinstance(val, dict):
                self._save_dict_latex(val, fd, names)

            elif isinstance(val, basestr):
                fd.write(self._escape_latex(val) + '\n')

            elif isinstance(val, float):
                fd.write('$' + self._format(val) + '$\n')

            elif isinstance(val, nm.ndarray):
                if val.ndim == 0:
                    fd.write('$' + self._format(val) + '$\n')

                elif val.ndim == 1:
                    self._write1d(fd, val)

                elif val.ndim == 2:
                    self._write2d(fd, val)

            else:
                fd.write('%s' % val)

        fd.write( r'\end{itemize}' )
        fd.write( '\n\n' )

    def test_term_call_modes(self):
        from sfepy.terms import term_table
        ok = True

        failed = []
        for domain in self.domains:
            self.report('domain: %s' % domain.name)

            domain_geometry = list(domain.geom_els.values())[0].name
            if domain.shape.dim != domain.shape.tdim:
                domain_geometry = '%d_%s' % (domain.shape.dim, domain_geometry)

            for _, term_cls in ordered_iteritems(term_table):
                if not domain_geometry in term_cls.geometries:
                    continue

                vint = ('volume', 'point', 'custom')
                rname = 'Omega' if term_cls.integration in vint else 'Gamma'

                self.report('<-- %s ...' % term_cls.name)

                if rname == 'Gamma' and domain.mesh.dim == 1:
                    self.report('--> 1D Gamma region: not tested!')

                elif term_cls.arg_shapes:
                    try:
                        _ok = self._test_single_term(term_cls, domain, rname)

                    except:
                        _ok = False

                    if not _ok:
                        failed.append((domain.name, term_cls.name))

                    ok = ok and _ok
                    self.report('--> ok: %s' % _ok)

                else:
                    self.report('--> not tested!')

        self.report('failed:', failed)

        return ok

    def test_term_call_modes(self):
        from sfepy.terms import term_table
        ok = True

        failed = []
        for domain in self.domains:
            self.report('domain: %s' % domain.name)

            for _, term_cls in ordered_iteritems(term_table):
                if not domain.groups[0].gel.name in term_cls.geometries:
                    continue

                rname = 'Omega' \
                        if term_cls.integration in ('volume', 'point') \
                        else 'Gamma'

                self.report('<-- %s ...' % term_cls.name)

                if rname == 'Gamma' and domain.mesh.dim == 1:
                    self.report('--> 1D Gamma region: not tested!')

                elif term_cls.arg_shapes:
                    try:
                        _ok = self._test_single_term(term_cls, domain, rname)

                    except:
                        _ok = False

                    if not _ok:
                        failed.append((domain.name, term_cls.name))

                    ok = ok and _ok
                    self.report('--> ok: %s' % _ok)

                else:
                    self.report('--> not tested!')

        self.report('failed:', failed)

        return ok

#.........这里部分代码省略.........
    pb.time_update(ebcs=Conditions(bc))
    pb.save_regions_as_groups('regions')

#ls = ScipyIterative({'method':'bicgstab','i_max':100,'eps_r':1e-10})


#   A = pb.mtx_a
#   M = spilu(A,fill_factor = 1)
    
    #conf = Struct(solvers ='ScipyIterative',method='bcgsl', sub_precond=None,
# i_max=1000, eps_r=1e-8)
        
#pb.set_conf_solvers(conf)
    vec = pb.solve()
    print nls_status
    file.write('TIME TO SOLVE\n')
    file.write(str(nls.status.time_stats['solve'])+'\n')
    file.write('TIME TO CREATE MATRIX\n')
    file.write(str(nls.status.time_stats['matrix'])+'\n')
    #out = post_process(out, pb, state, extend=False)
    ev = pb.evaluate
    out = vec.create_output_dict()
    strain = ev('ev_cauchy_strain.3.Omega(u)', mode='el_avg')
    stress = ev('ev_cauchy_stress.3.Omega(Mat.D, u)', mode='el_avg',
                copy_materials=False)

    out['cauchy_strain'] = Struct(name='output_data', mode='cell',
                                  data=strain, dofs=None)
    out['cauchy_stress'] = Struct(name='output_data', mode='cell',
                                  data=stress, dofs=None)


    # Postprocess the solution.
    #out = vec.create_output_dict()
    #out = stress_strain(out, pb, vec,lam,mu, extend=True)
    #pb.save_state('its2D_interactive.vtk', out=out)
    #print 'aqui estoy'
    pb.save_state('strain.vtk', out=out)
    #pb.save_state('disp.vtk', out=vec)
    #print 'ahora estoy aqui'
    #out = stress_strain(out, pb, vec, extend=True)
    #pb.save_state('out.vtk', out=out)
    print nls_status
    
    order = 3
    strain_qp = ev('ev_cauchy_strain.%d.Omega(u)' % order, mode='qp')
    stress_qp = ev('ev_cauchy_stress.%d.Omega(Mat.D, u)' % order,
                       mode='qp', copy_materials=False)

    file.close()
    options_probe=False
    if options_probe:
        # Probe the solution.
        probes, labels = gen_lines(pb)
        nls_options = {'eps_a':1e-8,'i_max':1}
        ls = ScipyDirect({})
        ls2 = ScipyIterative({'method':'bicgstab','i_max':5000,'eps_r':1e-20})
        order = 5
        sfield = Field.from_args('sym_tensor', nm.float64, (3,), omega,
                                approx_order=order-1)
        stress = FieldVariable('stress', 'parameter', sfield,
                               primary_var_name='(set-to-None)')
        strain = FieldVariable('strain', 'parameter', sfield,
                               primary_var_name='(set-to-None)')

        cfield = Field.from_args('component', nm.float64, 1, omega,
                                 approx_order=order-1)
        component = FieldVariable('component', 'parameter', cfield,
                                  primary_var_name='(set-to-None)')

        ev = pb.evaluate
        order = 2*(order - 1) #2 * (2- 1)
        print "before strain_qp"
        strain_qp = ev('ev_cauchy_strain.%d.Omega(u)' % order, mode='qp')
        stress_qp = ev('ev_cauchy_stress.%d.Omega(Mat.D, u)' % order,
                       mode='qp', copy_materials=False)
        print "before projections"
        print stress
        project_by_component(strain, strain_qp, component, order,ls2,nls_options)
        #print 'strain done'
        project_by_component(stress, stress_qp, component, order,ls2,nls_options)

        print "after projections"
        
        all_results = []
        for ii, probe in enumerate(probes):
            fig, results = probe_results2(u, strain, stress, probe, labels[ii])

            fig.savefig('test_probe_%d.png' % ii)
            all_results.append(results)

        for ii, results in enumerate(all_results):
            output('probe %d:' % ii)
            output.level += 2
            for key, res in ordered_iteritems(results):
                output(key + ':')
                val = res[1]
                output('  min: %+.2e, mean: %+.2e, max: %+.2e'
                       % (val.min(), val.mean(), val.max()))
            output.level -= 2