def __init__(self, conf, problem, **kwargs):
        NonlinearSolver.__init__(self, conf, **kwargs)

        conf = self.conf

        log = get_logging_conf(conf)
        conf.log = log = Struct(name='log_conf', **log)
        conf.is_any_log = (log.text is not None) or (log.plot is not None)

        conf.problem = problem

        conf = self.conf
        if conf.is_any_log:
            self.log = Log([[r'$||r||$'], ['iteration'],
                            [r'$\gamma$', r'$\max(\delta)$', r'$\max(\tau)$']],
                           xlabels=['', '', 'all iterations'],
                           ylabels=[r'$||r||$', 'iteration', 'stabilization'],
                           yscales=['log', 'linear', 'log'],
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%d'],
                                    ['%.8e', '%.8e', '%.8e']])

        else:
            self.log = None

def vary_incident_wave_dir( problem ):
    default_printer.prefix = 'vary_incident_wave_dir:'

    log_conf = {
        'is_plot' : True,
        'yscales' : ['linear'],
        'xaxes' : ['incident wave angle [degrees]'],
        'yaxes' : ['phase velocity [m/s]'],
    }
    dim = problem.domain.mesh.dim
    log = Log.from_conf( log_conf,
                         ['phase velocity %d' % ii for ii in range( dim )] )
    
    alphas = nm.linspace( 0.0, 360.0, 37 )
    output( 'running for angles:', alphas )
    pause()
    for ii, alpha in enumerate( alphas ):
        output( 'iteration %d: alpha %2f' % (ii, alpha) )
        opts = problem.conf.options

        ra = nm.pi * alpha / 180.0
        iwd = nm.zeros( (dim,), dtype = nm.float64 )
        iwd[0:2] = [nm.cos( ra ), nm.sin( ra )]
        opts.incident_wave_dir = iwd

        out = []
        yield problem, out
        
        #You have: sqrt( 10^10Pa / 10^4kg * m^3 )
        #You want:
        #Definition: 1000 m / s
        convert_to_ms = 1000.0
        phase_velocity = out[0] * convert_to_ms # to m/s.
        log( x = [alpha], *phase_velocity )

        output( 'phase velocity [m/s]:', phase_velocity )
        if opts.homogeneous:
            mat = problem.materials['matrix']
            # dilatation wave: vp = \sqrt{ (\lambda + 2\mu) / \rho }
            vd = nm.sqrt( (mat.D[0,0]) / mat.density) * convert_to_ms
            # shear wave: vs = \sqrt{ \mu / \rho }
            vs = nm.sqrt( mat.D[-1,-1] / mat.density ) * convert_to_ms
            output( 'analytical dilatation: %f, shear: %f' % (vd, vs) )
#            pause()

        yield None

    print log
    pause()

    if opts.homogeneous:
        name = 'homogeneous'
    else:
        name = 'perforated'
    fig_name = os.path.join( problem.output_dir,
                             'phase_velocity_%s%s' % (name, opts.fig_suffix) )
    log( save_figure = fig_name, finished = True )

    def __init__( self, conf, **kwargs ):
        OptimizationSolver.__init__( self, conf, **kwargs )

        conf = self.conf
        if conf.is_any_log:
            self.log = Log([[r'$||\Psi||$'], [r'$||\nabla \Psi||$'],
                            [r'$\alpha$'], ['iteration']],
                           xlabels=['', '', 'all iterations', 'all iterations'],
                           yscales=conf.yscales,
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%.3e'], ['%.3e'], ['%d']])
        else:
            self.log = None

    def __init__(self, conf, **kwargs):
        NonlinearSolver.__init__(self, conf, **kwargs)

        conf = self.conf
        if conf.is_any_log:
            self.log = Log([[r'$||r||$'], ['iteration']],
                           xlabels=['', 'all iterations'],
                           ylabels=[r'$||r||$', 'iteration'],
                           yscales=['log', 'linear'],
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%d']])

        else:
            self.log = None

    def __init__(self, conf, **kwargs):
        NonlinearSolver.__init__(self, conf, **kwargs)

        conf = self.conf
        if conf.is_any_log:
            self.log = Log(
                [[r"$||r||$"], ["iteration"]],
                xlabels=["", "all iterations"],
                ylabels=[r"$||r||$", "iteration"],
                yscales=["log", "linear"],
                is_plot=conf.log.plot is not None,
                log_filename=conf.log.text,
                formats=[["%.8e"], ["%d"]],
            )

        else:
            self.log = None

    def __init__(self, conf, **kwargs):
        NonlinearSolver.__init__(self, conf, **kwargs)

        conf = self.conf

        log = get_logging_conf(conf)
        conf.log = log = Struct(name='log_conf', **log)
        conf.is_any_log = (log.text is not None) or (log.plot is not None)

        if conf.is_any_log:
            self.log = Log([[r'$||r||$'], ['iteration']],
                           xlabels=['', 'all iterations'],
                           ylabels=[r'$||r||$', 'iteration'],
                           yscales=['log', 'linear'],
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%d']])

        else:
            self.log = None

def main():
    cwd = os.path.split(os.path.join(os.getcwd(), __file__))[0]

    log = Log(
        (["sin(x)", "cos(x)"], ["exp(x)"]),
        yscales=["linear", "log"],
        xlabels=["angle", None],
        ylabels=[None, "a function"],
        log_filename=os.path.join(cwd, "live_plot.log"),
    )
    log2 = Log(
        [["x^3"]],
        yscales=["linear"],
        xlabels=["x"],
        ylabels=["a cubic function"],
        log_filename=os.path.join(cwd, "live_plot2.log"),
    )

    added = 0
    for x in nm.linspace(0, 4.0 * nm.pi, 200):
        output("x: ", x)

        if x < (2.0 * nm.pi):
            log(nm.sin(x), nm.cos(x), nm.exp(x), x=[x, None])

        else:
            if added:
                log(nm.sin(x), nm.cos(x), nm.exp(x), x ** 2, x=[x, None, x])
            else:
                log.plot_vlines(color="r", linewidth=2)
                log.add_group(["x^2"], "linear", "new x", "square", formats=["%+g"])
            added += 1

        if (added == 20) or (added == 50):
            log.plot_vlines([2], color="g", linewidth=2)

        log2(x * x * x, x=[x])

    print log
    print log2
    pause()

    log(finished=True)
    log2(finished=True)

def main():

    log_conf = {
        'is_plot' : True,
        'aggregate' : 200,
        'yscales' : ['linear', 'log'],
        'xaxes' : ['angle', None],
        'yaxes' : [None, 'a function'],
    }

    log = Log.from_conf( log_conf, (['sin( x )', 'cos( x )'],['exp( x )']) )

    for x in nm.linspace( 0, 4.0 * nm.pi, 200 ):
        output( 'x: ', x )
        log( nm.sin( x ), nm.cos( x ), nm.exp( x ), x = [x, None] )

    print log
    pause()

    log( finished = True )

def main():
    cwd = os.path.split(os.path.join(os.getcwd(), __file__))[0]

    log = Log((['sin(x) + i sin(x**2)', 'cos(x)'], ['exp(x)']),
              yscales=['linear', 'log'],
              xlabels=['angle', None], ylabels=[None, 'a function'],
              log_filename=os.path.join(cwd, 'live_plot.log'))
    log2 = Log([['x^3']],
               yscales=['linear'],
               xlabels=['x'], ylabels=['a cubic function'],
               aggregate=50, sleep=0.05,
               log_filename=os.path.join(cwd, 'live_plot2.log'),
               formats=[['{:.5e}']])

    added = 0
    for x in nm.linspace(0, 4.0 * nm.pi, 200):
        output('x: ', x)

        if x < (2.0 * nm.pi):
            log(nm.sin(x)+1j*nm.sin(x**2), nm.cos(x), nm.exp(x), x=[x, None])

        else:
            if added:
                log(nm.sin(x)+1j*nm.sin(x**2), nm.cos(x), nm.exp(x), x**2,
                    x=[x, None, x])
            else:
                log.plot_vlines(color='r', linewidth=2)
                log.add_group(['x^2'], yscale='linear', xlabel='new x',
                              ylabel='square', formats=['%+g'])
            added += 1

        if (added == 20) or (added == 50):
            log.plot_vlines([2], color='g', linewidth=2)

        log2(x*x*x, x=[x])

    print(log)
    print(log2)
    pause()

    log(finished=True)
    log2(finished=True)

class FMinSteepestDescent( OptimizationSolver ):
    name = 'opt.fmin_sd'

    def process_conf( conf ):
        """
        Missing items are set to default values.
        
        Example configuration, all items::
        
            solver_0 = {
                'name'      : 'fmin_sd',
                'kind'      : 'opt.fmin_sd',

                'i_max'      : 10,
                'eps_rd'     : 1e-5, # Relative delta of objective function
                'eps_of'     : 1e-4,
                'eps_ofg'    : 1e-8,
                'norm'      : nm.Inf,
                'ls'        : True, # Linesearch.
                'ls_method'  : 'backtracking', # 'backtracking' or 'full'
                'ls0'       : 0.25,
                'ls_red'     : 0.5,
                'ls_red_warp' : 0.1,
                'ls_on'      : 0.99999,
                'ls_min'     : 1e-5,
                'check'     : 0,
                'delta'     : 1e-6,
                'output'    : None, # 'itc'
                'log'       : {'text' : 'output/log.txt',
                               'plot' : 'output/log.png'},
                'yscales'   : ['linear', 'log', 'log', 'linear'],
            }
        """
        get = conf.get_default_attr

        i_max = get( 'i_max', 10 )
        eps_rd = get( 'eps_rd', 1e-5 )
        eps_of = get( 'eps_of', 1e-4 )
        eps_ofg = get( 'eps_ofg', 1e-8 )
        norm = get( 'norm', nm.Inf )
        ls = get( 'ls', True )
        ls_method = get( 'ls_method', 'backtracking' )
        ls0 = get( 'ls0', 0.25 )
        ls_red = get( 'ls_red', 0.5 )
        ls_red_warp = get( 'ls_red_warp', 0.1 )
        ls_on = get( 'ls_on', 0.99999 )
        ls_min = get( 'ls_min', 1e-5 )
        check = get( 'check', 0 )
        delta = get( 'delta', 1e-6)
        output = get( 'output', None )
        yscales = get( 'yscales', ['linear', 'log', 'log', 'linear'] )

        log = get_logging_conf(conf)
        log = Struct(name='log_conf', **log)
        is_any_log = (log.text is not None) or (log.plot is not None)

        common = OptimizationSolver.process_conf( conf )
        return Struct( **locals() ) + common
    process_conf = staticmethod( process_conf )

    ##
    # 17.10.2007, c
    def __init__( self, conf, **kwargs ):
        OptimizationSolver.__init__( self, conf, **kwargs )

        conf = self.conf
        if conf.is_any_log:
            self.log = Log([[r'$||\Psi||$'], [r'$||\nabla \Psi||$'],
                            [r'$\alpha$'], ['iteration']],
                           xlabels=['', '', 'all iterations', 'all iterations'],
                           yscales=conf.yscales,
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%.3e'], ['%.3e'], ['%d']])
        else:
            self.log = None

    ##
    # 19.04.2006, c
    # 20.04.2006
    # 21.04.2006
    # 26.04.2006
    # 06.06.2006
    # 07.06.2006
    # 04.09.2006
    # 21.03.2007
    # 17.10.2007, from fmin_sd()
    def __call__( self, x0, conf = None, obj_fun = None, obj_fun_grad = None,
                  status = None, obj_args = None ):
#    def fmin_sd( conf, x0, fn_of, fn_ofg, args = () ):

        conf = get_default( conf, self.conf )
        obj_fun = get_default( obj_fun, self.obj_fun )
        obj_fun_grad = get_default( obj_fun_grad, self.obj_fun_grad )
        status = get_default( status, self.status )
        obj_args = get_default( obj_args, self.obj_args )

        if conf.output:
            globals()['output'] = conf.output

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

    def solve_eigen_problem_n( self ):
        opts = self.app_options
        pb = self.problem

        dim = pb.domain.mesh.dim

        pb.set_equations( pb.conf.equations )
        pb.select_bcs( ebc_names = ['ZeroSurface'] )

        output( 'assembling rhs...' )
        tt = time.clock()
        mtx_b = pb.evaluate(pb.conf.equations['rhs'], mode='weak',
                            auto_init=True, dw_mode='matrix')
        output( '...done in %.2f s' % (time.clock() - tt) )
        assert_( nm.alltrue( nm.isfinite( mtx_b.data ) ) )

        ## mtx_b.save( 'b.txt', format='%d %d %.12f\n' )

        aux = pb.create_evaluable(pb.conf.equations['lhs'], mode='weak',
                                  dw_mode='matrix')
        mtx_a_equations, mtx_a_variables = aux

        if self.options.plot:
            log_conf = {
                'is_plot' : True,
                'aggregate' : 1,
                'yscales' : ['linear', 'log'],
            }
        else:
            log_conf = {
                'is_plot' : False,
            }
        log =  Log.from_conf( log_conf, ([r'$|F(x)|$'], [r'$|F(x)-x|$']) )

        file_output = Output('', opts.log_filename, combined = True)

        eig_conf = pb.get_solver_conf( opts.eigen_solver )
        eig_solver = Solver.any_from_conf( eig_conf )

        # Just to get the shape. Assumes one element group only!!!
        v_hxc_qp = pb.evaluate('dq_state_in_volume_qp.i1.Omega(Psi)')
        v_hxc_qp.fill(0.0)
        self.qp_shape = v_hxc_qp.shape
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)

        self.norm_v_hxc0 = nla.norm(vec_v_hxc)
        self.itercount = 0
        aux = wrap_function(self.iterate,
                            (eig_solver,
                             mtx_a_equations, mtx_a_variables,
                             mtx_b, log, file_output))
        ncalls, times, nonlin_v, results = aux

        # Create and call the DFT solver.
        dft_conf = pb.get_solver_conf(opts.dft_solver)
        dft_status = {}
        dft_solver = Solver.any_from_conf(dft_conf,
                                          fun = nonlin_v,
                                          status = dft_status)
        v_hxc_qp = dft_solver(v_hxc_qp.ravel())

        v_hxc_qp = nm.array(v_hxc_qp, dtype=nm.float64)
        v_hxc_qp.shape = self.qp_shape
        eigs, mtx_s_phi, vec_n, vec_v_h, v_ion_qp, v_xc_qp, v_hxc_qp = results
        output( 'DFT iteration time [s]:', dft_status['time_stats'] )

        fun = pb.materials['mat_v'].function
        variable = self.problem.create_variables(['scalar'])['scalar']
        vec_v_ion = fun(None, variable.field.get_coor(),
                        mode='qp')['V_ion'].squeeze()

        vec_v_xc = self._interp_to_nodes(v_xc_qp)
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)
        vec_v_sum = self._interp_to_nodes(v_hxc_qp + v_ion_qp)

        coor = pb.domain.get_mesh_coors()
        r2 = norm_l2_along_axis(coor, squared=True)
        vec_nr2 = vec_n * r2

        pb.select_bcs( ebc_names = ['ZeroSurface'] )
        mtx_phi = self.make_full( mtx_s_phi )

        out = {}
        update_state_to_output(out, pb, vec_n, 'n')
        update_state_to_output(out, pb, vec_nr2, 'nr2')
        update_state_to_output(out, pb, vec_v_h, 'V_h')
        update_state_to_output(out, pb, vec_v_xc, 'V_xc')
        update_state_to_output(out, pb, vec_v_ion, 'V_ion')
        update_state_to_output(out, pb, vec_v_hxc, 'V_hxc')
        update_state_to_output(out, pb, vec_v_sum, 'V_sum')
        self.save_results(eigs, mtx_phi, out=out)

        if self.options.plot:
            log( save_figure = opts.iter_fig_name )
            pause()
            log(finished=True)

        return Struct( pb = pb, eigs = eigs, mtx_phi = mtx_phi,
                       vec_n = vec_n, vec_nr2 = vec_nr2,
                       vec_v_h = vec_v_h, vec_v_xc = vec_v_xc )

def solve_eigen_problem_n( conf, options ):

    pb = ProblemDefinition.from_conf( conf )
    dim = pb.domain.mesh.dim

    pb.time_update()

    dummy = pb.create_state_vector()

    output( 'assembling rhs...' )
    tt = time.clock()
    mtx_b = eval_term_op( dummy, conf.equations['rhs'], pb,
                       dw_mode = 'matrix', tangent_matrix = pb.mtx_a.copy() )
    output( '...done in %.2f s' % (time.clock() - tt) )

    #mtxA.save( 'tmp/a.txt', format='%d %d %.12f\n' )
    #mtxB.save( 'tmp/b.txt', format='%d %d %.12f\n' )

    try:
        n_eigs = conf.options.n_eigs
    except AttributeError:
        n_eigs = mtx_a.shape[0]

    if n_eigs is None:
        n_eigs = mtx_a.shape[0]

##     mtx_a.save( 'a.txt', format='%d %d %.12f\n' )
##     mtx_b.save( 'b.txt', format='%d %d %.12f\n' )

    if options.plot:
        log_conf = {
            'is_plot' : True,
            'aggregate' : 1,
            'yscales' : ['linear', 'log'],
        }
    else:
        log_conf = {
            'is_plot' : False,
        }
    log =  Log.from_conf( log_conf, ([r'$|F(x)|$'], [r'$|F(x)|$']) )

    eig_conf = pb.get_solver_conf( conf.options.eigen_solver )
    eig_solver = Solver.any_from_conf( eig_conf )
    vec_vhxc = nm.zeros( (pb.variables.di.ptr[-1],), dtype = nm.float64 )
    aux = wrap_function( iterate,
                         (pb, conf, eig_solver, n_eigs, mtx_b, log) )
    ncalls, times, nonlin_v = aux
    vec_vhxc = broyden3( nonlin_v, vec_vhxc, verbose = True )
    out = iterate( vec_vhxc, pb, conf, eig_solver, n_eigs, mtx_b )
    eigs, mtx_s_phi, vec_n, vec_vh, vec_vxc = out

    if options.plot:
        log( finished = True )
        pause()
        
    coor = pb.domain.get_mesh_coors()
    r = coor[:,0]**2 + coor[:,1]**2 + coor[:,2]**2
    vec_nr2 = vec_n * r

    n_eigs = eigs.shape[0]

    mtx_phi = nm.empty( (pb.variables.di.ptr[-1], mtx_s_phi.shape[1]),
                       dtype = nm.float64 )
    for ii in xrange( n_eigs ):
        mtx_phi[:,ii] = pb.variables.make_full_vec( mtx_s_phi[:,ii] )

    save = get_default_attr( conf.options, 'save_eig_vectors', None )
    out = {}
    for ii in xrange( n_eigs ):
        if save is not None:
            if (ii >= save[0]) and (ii < (n_eigs - save[1])): continue
        aux = pb.state_to_output( mtx_phi[:,ii] )
        key = aux.keys()[0]
        out[key+'%03d' % ii] = aux[key]

    update_state_to_output( out, pb, vec_n, 'n' )
    update_state_to_output( out, pb, vec_nr2, 'nr2' )
    update_state_to_output( out, pb, vec_vh, 'vh' )
    update_state_to_output( out, pb, vec_vxc, 'vxc' )

    ofn_trunk = options.output_filename_trunk
    pb.domain.mesh.write( ofn_trunk + '.vtk', io = 'auto', out = out )

    fd = open( ofn_trunk + '_eigs.txt', 'w' )
    eigs.tofile( fd, ' ' )
    fd.close()

    return Struct( pb = pb, eigs = eigs, mtx_phi = mtx_phi )

class Oseen(NonlinearSolver):
    """
    The Oseen solver for Navier-Stokes equations.
    """
    name = 'nls.oseen'

    __metaclass__ = SolverMeta

    _parameters = [
        ('stabil_mat', 'str', None, True,
         'The name of stabilization material.'),
        ('adimensionalize', 'bool', False, False,
         'If True, adimensionalize the problem (not implemented!).'),
        ('check_navier_stokes_rezidual', 'bool', False, False,
         'If True, check the Navier-Stokes rezidual after the nonlinear loop.'),
        ('i_max', 'int', 1, False,
         'The maximum number of iterations.'),
        ('eps_a', 'float', 1e-10, False,
         'The absolute tolerance for the residual, i.e. :math:`||f(x^i)||`.'),
        ('eps_r', 'float', 1.0, False,
         """The relative tolerance for the residual, i.e. :math:`||f(x^i)|| /
            ||f(x^0)||`."""),
        ('macheps', 'float', nm.finfo(nm.float64).eps, False,
         'The float considered to be machine "zero".'),
        ('lin_red', 'float', 1.0, False,
         """The linear system solution error should be smaller than (`eps_a` *
            `lin_red`), otherwise a warning is printed."""),
        ('lin_precision', 'float or None', None, False,
         """If not None, the linear system solution tolerances are set in each
            nonlinear iteration relative to the current residual norm by the
            `lin_precision` factor. Ignored for direct linear solvers."""),
    ]

    def __init__(self, conf, problem, **kwargs):
        NonlinearSolver.__init__(self, conf, **kwargs)

        conf = self.conf

        log = get_logging_conf(conf)
        conf.log = log = Struct(name='log_conf', **log)
        conf.is_any_log = (log.text is not None) or (log.plot is not None)

        conf.problem = problem

        conf = self.conf
        if conf.is_any_log:
            self.log = Log([[r'$||r||$'], ['iteration'],
                            [r'$\gamma$', r'$\max(\delta)$', r'$\max(\tau)$']],
                           xlabels=['', '', 'all iterations'],
                           ylabels=[r'$||r||$', 'iteration', 'stabilization'],
                           yscales=['log', 'linear', 'log'],
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%d'],
                                    ['%.8e', '%.8e', '%.8e']])

        else:
            self.log = None

    def __call__(self, vec_x0, conf=None, fun=None, fun_grad=None,
                 lin_solver=None, status=None, problem=None):
        """
        Oseen solver is problem-specific - it requires a Problem instance.
        """
        conf = get_default(conf, self.conf)
        fun = get_default(fun, self.fun)
        fun_grad = get_default(fun_grad, self.fun_grad)
        lin_solver = get_default(lin_solver, self.lin_solver)
        status = get_default(status, self.status)
        problem = get_default(problem, conf.problem,
                              '`problem` parameter needs to be set!')

        time_stats = {}

        stabil = problem.get_materials()[conf.stabil_mat]
        ns, ii = stabil.function.function.get_maps()

        variables = problem.get_variables()
        update_var = variables.set_data_from_state
        make_full_vec = variables.make_full_vec

        output('problem size:')
        output('    velocity: %s' % ii['us'])
        output('    pressure: %s' % ii['ps'])

        vec_x = vec_x0.copy()
        vec_x_prev = vec_x0.copy()
        vec_dx = None

        if self.log is not None:
            self.log.plot_vlines(color='r', linewidth=1.0)

        err0 = -1.0
        it = 0
        while 1:
            vec_x_prev_f = make_full_vec(vec_x_prev)
            update_var(ns['b'], vec_x_prev_f, ns['u'])

            vec_b = vec_x_prev_f[ii['u']]
            b_norm = nla.norm(vec_b, nm.inf)
#.........这里部分代码省略.........

class FMinSteepestDescent(OptimizationSolver):
    """
    Steepest descent optimization solver.
    """
    name = 'opt.fmin_sd'

    __metaclass__ = SolverMeta

    _parameters = [
        ('i_max', 'int', 10, False,
         'The maximum number of iterations.'),
        ('eps_rd', 'float', 1e-5, False,
         'The relative delta of the objective function.'),
        ('eps_of', 'float', 1e-4, False,
         'The tolerance for the objective function.'),
        ('eps_ofg', 'float', 1e-8, False,
         'The tolerance for the objective function gradient.'),
        ('norm', 'numpy norm', nm.Inf, False,
         'The norm to be used.'),
        ('ls', 'bool', True, False,
         'If True, use a line-search.'),
        ('ls_method', "{'backtracking', 'full'}", 'backtracking', False,
         'The line-search method.'),
        ('ls_on', 'float', 0.99999, False,
         """Start the backtracking line-search by reducing the step, if
            :math:`||f(x^i)|| / ||f(x^{i-1})||` is larger than `ls_on`."""),
        ('ls0', '0.0 < float < 1.0', 1.0, False,
         'The initial step.'),
        ('ls_red', '0.0 < float < 1.0', 0.5, False,
         'The step reduction factor in case of correct residual assembling.'),
        ('ls_red_warp', '0.0 < float < 1.0', 0.1, False,
         """The step reduction factor in case of failed residual assembling
            (e.g. the "warp violation" error caused by a negative volume
            element resulting from too large deformations)."""),
        ('ls_min', '0.0 < float < 1.0', 1e-5, False,
         'The minimum step reduction factor.'),
        ('check', '0, 1 or 2', 0, False,
         """If >= 1, check the tangent matrix using finite differences.  If 2,
            plot the resulting sparsity patterns."""),
        ('delta', 'float', 1e-6, False,
         r"""If `check >= 1`, the finite difference matrix is taken as
            :math:`A_{ij} = \frac{f_i(x_j + \delta) - f_i(x_j - \delta)}{2
            \delta}`."""),
        ('output', 'function', None, False,
         """If given, use it instead of :func:`output()
            <sfepy.base.base.output()>` function."""),
        ('yscales', 'list of str', ['linear', 'log', 'log', 'linear'], False,
         'The list of four convergence log subplot scales.'),
        ('log', 'dict or None', None, False,
         """If not None, log the convergence according to the configuration in
            the following form: ``{'text' : 'log.txt', 'plot' : 'log.pdf'}``.
            Each of the dict items can be None."""),
    ]

    def __init__(self, conf, **kwargs):
        OptimizationSolver.__init__(self, conf, **kwargs)

        conf = self.conf

        log = get_logging_conf(conf)
        conf.log = log = Struct(name='log_conf', **log)
        conf.is_any_log = (log.text is not None) or (log.plot is not None)

        if conf.is_any_log:
            self.log = Log([[r'$||\Psi||$'], [r'$||\nabla \Psi||$'],
                            [r'$\alpha$'], ['iteration']],
                           xlabels=['', '', 'all iterations', 'all iterations'],
                           yscales=conf.yscales,
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%.3e'], ['%.3e'], ['%d']])
        else:
            self.log = None

    def __call__(self, x0, conf=None, obj_fun=None, obj_fun_grad=None,
                 status=None, obj_args=None):
        conf = get_default(conf, self.conf)
        obj_fun = get_default(obj_fun, self.obj_fun)
        obj_fun_grad = get_default(obj_fun_grad, self.obj_fun_grad)
        status = get_default(status, self.status)
        obj_args = get_default(obj_args, self.obj_args)

        if conf.output:
            globals()['output'] = conf.output

        output('entering optimization loop...')

        nc_of, tt_of, fn_of = wrap_function(obj_fun, obj_args)
        nc_ofg, tt_ofg, fn_ofg = wrap_function(obj_fun_grad, obj_args)

        time_stats = {'of' : tt_of, 'ofg': tt_ofg, 'check' : []}

        ofg = None

        it = 0
        xit = x0.copy()
        while 1:

            of = fn_of(xit)

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

class Newton( NonlinearSolver ):
    name = 'nls.newton'

    def process_conf( conf ):
        """
        Missing items are set to default values for a linear problem.
        
        Example configuration, all items::
        
            solver_1 = {
                'name' : 'newton',
                'kind' : 'nls.newton',

                'i_max'      : 2,
                'eps_a'      : 1e-8,
                'eps_r'      : 1e-2,
                'macheps'   : 1e-16,
                'lin_red'    : 1e-2, # Linear system error < (eps_a * lin_red).
                'ls_red'     : 0.1,
                'ls_red_warp' : 0.001,
                'ls_on'      : 0.99999,
                'ls_min'     : 1e-5,
                'check'     : 0,
                'delta'     : 1e-6,
                'is_plot'    : False,
                'log'        : None,
                 # 'nonlinear' or 'linear' (ignore i_max)
                'problem'   : 'nonlinear',
            }
        """
        get = conf.get_default_attr

        i_max = get( 'i_max', 1 )
        eps_a = get( 'eps_a', 1e-10 )
        eps_r = get( 'eps_r', 1.0 )
        macheps = get( 'macheps', nm.finfo( nm.float64 ).eps )
        lin_red = get( 'lin_red', 1.0 )
        ls_red = get( 'ls_red', 0.1 )
        ls_red_warp = get( 'ls_red_warp', 0.001 )
        ls_on = get( 'ls_on', 0.99999 )
        ls_min = get( 'ls_min', 1e-5 )
        check = get( 'check', 0 )
        delta = get( 'delta', 1e-6)
        is_plot = get( 'is_plot', False )
        problem = get( 'problem', 'nonlinear' )

        log = get_logging_conf(conf)
        log = Struct(name='log_conf', **log)
        is_any_log = (log.text is not None) or (log.plot is not None)

        common = NonlinearSolver.process_conf( conf )
        return Struct( **locals() ) + common
    process_conf = staticmethod( process_conf )

    def __init__(self, conf, **kwargs):
        NonlinearSolver.__init__( self, conf, **kwargs )

        conf = self.conf
        if conf.is_any_log:
            self.log = Log([[r'$||r||$'], ['iteration']],
                           xlabels=['', 'all iterations'],
                           ylabels=[r'$||r||$', 'iteration'],
                           yscales=['log', 'linear'],
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%d']])

        else:
            self.log = None

    ##
    # c: 02.12.2005, r: 04.04.2008
    # 10.10.2007, from newton()
    def __call__( self, vec_x0, conf = None, fun = None, fun_grad = None,
                  lin_solver = None, status = None ):
        """setting conf.problem == 'linear' means 1 iteration and no rezidual
        check!
        """
        import sfepy.base.plotutils as plu
        conf = get_default( conf, self.conf )
        fun = get_default( fun, self.fun )
        fun_grad = get_default( fun_grad, self.fun_grad )
        lin_solver = get_default( lin_solver, self.lin_solver )
        status = get_default( status, self.status )

        time_stats = {}

        vec_x = vec_x0.copy()
        vec_x_last = vec_x0.copy()
        vec_dx = None

        if self.log is not None:
            self.log.plot_vlines(color='r', linewidth=1.0)

        err0 = -1.0
        err_last = -1.0
        it = 0
        while 1:

            ls = 1.0
#.........这里部分代码省略.........

def solve_navier_stokes(conf, options):
    opts = conf.options

    dpb = ProblemDefinition.from_conf(conf, init_equations=False)
    equations = getattr(conf, "_".join(("equations_direct", opts.problem)))
    dpb.set_equations(equations)

    ls_conf = dpb.get_solver_conf(opts.ls)
    nls_conf = dpb.get_solver_conf(opts.nls_direct)

    method = opts.direct_method
    if method == "stationary":
        data = {}
        dpb.time_update(None)
        state_dp = dpb.solve(nls_conf=nls_conf)

    elif method == "transient":
        ls = Solver.any_from_conf(ls_conf)
        ts_conf = dpb.get_solver_conf(opts.ts_direct)

        data = {"ts": Struct(dt=ts_conf.dt)}

        # Plug in mass term.
        mequations = {}
        for key, eq in equations.iteritems():
            if "dw_div_grad" in eq:
                eq = "+".join((ts_conf.mass_term, eq)).replace("++", "+")
            mequations[key] = eq

        if ts_conf.stokes_init:
            state_dp0 = solve_stokes(dpb, conf.equations_direct_stokes, nls_conf)
            dpb.set_equations(mequations)
        else:
            dpb.set_equations(mequations)
            state_dp0 = dpb.create_state()
            dpb.time_update(None)
            state_dp0.apply_ebc()

        from sfepy.base.log import Log

        log = Log.from_conf(Struct(is_plot=True), ([r"$||u||$"], [r"$||p||$"]))

        output("Navier-Stokes...")
        ev = BasicEvaluator(dpb, ts=Struct(dt=ts_conf.dt))
        nls = Solver.any_from_conf(nls_conf, evaluator=ev, lin_solver=ls)

        n_step = ts_conf.n_step
        step = 0
        while 1:
            for ii in xrange(n_step):
                output(step)

                vec_u = state_dp0("w")
                vec_p = state_dp0("r")
                log(nm.linalg.norm(vec_u), nm.linalg.norm(vec_p))

                dpb.variables.set_data_from_state("w_0", state_dp0(), "w")
                vec_dp = nls(state_dp0())

                step += 1
                state_dp = state_dp0.copy()
                state_dp.set_reduced(vec_dp)

                state_dp0 = state_dp

            if ts_conf.interactive:
                try:
                    n_step = int(raw_input("continue: "))
                    if n_step <= 0:
                        break
                except:
                    break

        vec_u = state_dp("w")
        vec_p = state_dp("r")
        log(nm.linalg.norm(vec_u), nm.linalg.norm(vec_p), finished=True)

    else:
        raise "unknown Navier-Stokes solution method (%s)!" % method

    return dpb, state_dp, data

class Newton(NonlinearSolver):
    r"""
    Solves a nonlinear system :math:`f(x) = 0` using the Newton method with
    backtracking line-search, starting with an initial guess :math:`x^0`.
    """
    name = 'nls.newton'

    __metaclass__ = SolverMeta

    _parameters = [
        ('i_max', 'int', 1, False,
         'The maximum number of iterations.'),
        ('eps_a', 'float', 1e-10, False,
         'The absolute tolerance for the residual, i.e. :math:`||f(x^i)||`.'),
        ('eps_r', 'float', 1.0, False,
         """The relative tolerance for the residual, i.e. :math:`||f(x^i)|| /
            ||f(x^0)||`."""),
        ('eps_mode', "'and' or 'or'", 'and', False,
         """The logical operator to use for combining the absolute and relative
            tolerances."""),
        ('macheps', 'float', nm.finfo(nm.float64).eps, False,
         'The float considered to be machine "zero".'),
        ('lin_red', 'float', 1.0, False,
         """The linear system solution error should be smaller than (`eps_a` *
            `lin_red`), otherwise a warning is printed."""),
        ('lin_precision', 'float or None', None, False,
         """If not None, the linear system solution tolerances are set in each
            nonlinear iteration relative to the current residual norm by the
            `lin_precision` factor. Ignored for direct linear solvers."""),
        ('ls_on', 'float', 0.99999, False,
         """Start the backtracking line-search by reducing the step, if
            :math:`||f(x^i)|| / ||f(x^{i-1})||` is larger than `ls_on`."""),
        ('ls_red', '0.0 < float < 1.0', 0.1, False,
         'The step reduction factor in case of correct residual assembling.'),
        ('ls_red_warp', '0.0 < float < 1.0', 0.001, False,
         """The step reduction factor in case of failed residual assembling
            (e.g. the "warp violation" error caused by a negative volume
            element resulting from too large deformations)."""),
        ('ls_min', '0.0 < float < 1.0', 1e-5, False,
         'The minimum step reduction factor.'),
        ('give_up_warp', 'bool', False, False,
         'If True, abort on the "warp violation" error.'),
        ('check', '0, 1 or 2', 0, False,
         """If >= 1, check the tangent matrix using finite differences.  If 2,
            plot the resulting sparsity patterns."""),
        ('delta', 'float', 1e-6, False,
         r"""If `check >= 1`, the finite difference matrix is taken as
            :math:`A_{ij} = \frac{f_i(x_j + \delta) - f_i(x_j - \delta)}{2
            \delta}`."""),
        ('log', 'dict or None', None, False,
         """If not None, log the convergence according to the configuration in
            the following form: ``{'text' : 'log.txt', 'plot' : 'log.pdf'}``.
            Each of the dict items can be None."""),
        ('is_linear', 'bool', False, False,
         'If True, the problem is considered to be linear.'),
    ]

    def __init__(self, conf, **kwargs):
        NonlinearSolver.__init__(self, conf, **kwargs)

        conf = self.conf

        log = get_logging_conf(conf)
        conf.log = log = Struct(name='log_conf', **log)
        conf.is_any_log = (log.text is not None) or (log.plot is not None)

        if conf.is_any_log:
            self.log = Log([[r'$||r||$'], ['iteration']],
                           xlabels=['', 'all iterations'],
                           ylabels=[r'$||r||$', 'iteration'],
                           yscales=['log', 'linear'],
                           is_plot=conf.log.plot is not None,
                           log_filename=conf.log.text,
                           formats=[['%.8e'], ['%d']])

        else:
            self.log = None

    def __call__(self, vec_x0, conf=None, fun=None, fun_grad=None,
                 lin_solver=None, iter_hook=None, status=None):
        """
        Nonlinear system solver call.

        Solves a nonlinear system :math:`f(x) = 0` using the Newton method with
        backtracking line-search, starting with an initial guess :math:`x^0`.

        Parameters
        ----------
        vec_x0 : array
            The initial guess vector :math:`x_0`.
        conf : Struct instance, optional
            The solver configuration parameters,
        fun : function, optional
            The function :math:`f(x)` whose zero is sought - the residual.
        fun_grad : function, optional
            The gradient of :math:`f(x)` - the tangent matrix.
        lin_solver : LinearSolver instance, optional
            The linear solver for each nonlinear iteration.
        iter_hook : function, optional
            User-supplied function to call before each iteration.
#.........这里部分代码省略.........

class Newton(NonlinearSolver):
    name = "nls.newton"

    @staticmethod
    def process_conf(conf, kwargs):
        """
        Missing items are set to default values for a linear problem.

        Example configuration, all items::

            solver_1 = {
                'name' : 'newton',
                'kind' : 'nls.newton',

                'i_max' : 2,
                'eps_a' : 1e-8,
                'eps_r' : 1e-2,
                'macheps' : 1e-16,