def get_conductivity(ts, coors, problem, equations = None, mode = None, **kwargs):
    """
    Calculates the conductivity with a constitutive k(psi) relation,
    where psi = -h.
    """
    if mode == 'qp':

        ## Get pressure values
        h_values = problem.evaluate('ev_volume_integrate.i.Omega(h)',
                                    mode = 'qp', verbose = False) * scaling
        psi_values = -h_values

        # van Genuchten
        val = vanGenuchten(ksat = silt_ksat, aVG = silt_aVG, nVG = silt_nVG,
                           mVG = silt_mVG, lVG = silt_lVG, psi = psi_values)

        # Brooks and Corey
        #val = brooksCorey(ksat = silt_ksat, aev = silt_aev, lBC = silt_lBC,
        #                  psi = psi_values)

        # Reshape the val vector to match SfePy expectations
        val.shape = (val.shape[0] * val.shape[1], 1, 1)

        # Check output
        output('h_values: min:', h_values.min(), 'max:', h_values.max())
        output('conductivity: min:', val.min(), 'max:', val.max())

        return {'val' : val}

    def create_conn_graph(self, verbose=True):
        """
        Create a graph of mesh connectivity.

        Returns
        -------
        graph : csr_matrix
            The mesh connectivity graph as a SciPy CSR matrix.
        """
        from extmods.cmesh import create_mesh_graph

        shape = (self.n_nod, self.n_nod)
        output('graph shape:', shape, verbose=verbose)
        if nm.prod(shape) == 0:
            output('no graph (zero size)!', verbose=verbose)
            return None

        output('assembling mesh graph...', verbose=verbose)
        tt = time.clock()

        nnz, prow, icol = create_mesh_graph(shape[0], shape[1],
                                            len(self.conns),
                                            self.conns, self.conns)
        output('...done in %.2f s' % (time.clock() - tt), verbose=verbose)
        output('graph nonzeros: %d (%.2e%% fill)' \
               % (nnz, float(nnz) / nm.prod(shape)))

        data = nm.ones((nnz,), dtype=nm.bool)
        graph = sp.csr_matrix((data, icol, prow), shape)

        return graph

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

    output.prefix = "therel:"

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

    problem = Problem.from_conf(conf, init_equations=False)

    # Setup output directory according to options above.
    problem.setup_default_output()

    # First solve the stationary electric conduction problem.
    problem.set_equations({"eq": conf.equations["1"]})
    problem.time_update()
    state_el = problem.solve()
    problem.save_state(problem.get_output_name(suffix="el"), state_el)

    # Then solve the evolutionary heat conduction problem, using state_el.
    problem.set_equations({"eq": conf.equations["2"]})
    phi_var = problem.get_variables()["phi_known"]
    phi_var.set_data(state_el())
    time_solver = problem.get_time_solver()
    time_solver()

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

def refine_mesh(filename, level):
    """
    Uniformly refine `level`-times a mesh given by `filename`.

    The refined mesh is saved to a file with name constructed from base
    name of `filename` and `level`-times appended `'_r'` suffix.

    Parameters
    ----------
    filename : str
        The mesh file name.
    level : int
        The refinement level.
    """
    import os
    from sfepy.base.base import output
    from sfepy.fem import Mesh, Domain

    if level > 0:
        mesh = Mesh.from_file(filename)
        domain = Domain(mesh.name, mesh)
        for ii in range(level):
            output('refine %d...' % ii)
            domain = domain.refine()
            output('... %d nodes %d elements'
                   % (domain.shape.n_nod, domain.shape.n_el))

        suffix = os.path.splitext(filename)[1]
        filename = domain.name + suffix

        domain.mesh.write(filename, io='auto')

    return filename

def main():
    from sfepy.base.base import output
    from sfepy.base.conf import ProblemConf, get_standard_keywords
    from sfepy.fem import ProblemDefinition
    from sfepy.applications import solve_evolutionary

    output.prefix = 'therel:'

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

    problem = ProblemDefinition.from_conf(conf, init_equations=False)

    # Setup output directory according to options above.
    problem.setup_default_output()

    # First solve the stationary electric conduction problem.
    problem.set_equations({'eq' : conf.equations['1']})
    problem.time_update()
    state_el = problem.solve()
    problem.save_state(problem.get_output_name(suffix = 'el'), state_el)

    # Then solve the evolutionary heat conduction problem, using state_el.
    problem.set_equations({'eq' : conf.equations['2']})
    phi_var = problem.get_variables()['phi_known']
    phi_var.data_from_any(state_el())
    solve_evolutionary(problem)

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

    def __call__(self, rhs, x0=None, conf=None, eps_a=None, eps_r=None,
                 i_max=None, mtx=None, status=None, **kwargs):

        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

        # 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)

        # Set PETSc rhs, solve, get solution from PETSc solution.
        if x0 is not None:
            psol[...] = x0
            ksp.setInitialGuessNonzero(True)
        prhs[...] = rhs
        ksp.solve(prhs, psol)
        sol = psol[...].copy()
        output('%s(%s) convergence: %s (%s)'
               % (self.conf.method, self.conf.precond,
                  ksp.reason, self.converged_reasons[ksp.reason]))

        return sol

    def __call__(self, mtx_a, mtx_b=None, n_eigs=None,
                 eigenvectors=None, status=None, conf=None):
        from pysparse import jdsym, itsolvers, precon

        output("loading...")
        A = self._convert_mat(mtx_a)
        output("...done")
        if mtx_b is not None:
            M = self._convert_mat(mtx_b)

        output("solving...")
        Atau=A.copy()
        Atau.shift(-conf.tau,M)
        K=precon.jacobi(Atau)
        A=A.to_sss();
        if mtx_b is not None:
            M=M.to_sss();

        method = getattr(itsolvers, conf.method)
        kconv, lmbd, Q, it, it_in = jdsym.jdsym(A, M, K, n_eigs, conf.tau,
                                                conf.eps_a, conf.i_max,
                                                method,
                                                clvl=conf.verbosity,
                                                strategy=conf.strategy)

        output("number of converged eigenvalues:", kconv)
        output("...done")

        if status is not None:
            status['q'] = Q
            status['it'] = it
            status['it_in'] = it_in

        return lmbd, Q

def main():
    parser = ArgumentParser(description=__doc__)
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('--eps', action='store', dest='eps',
                        default=1e-12, help=helps['eps'])
    parser.add_argument('-o', '--filename-out',
                        action='store', dest='filename_out',
                        default=None, help=helps['filename-out'])
    parser.add_argument('filename')
    options = parser.parse_args()

    filename = options.filename

    mesh = Mesh.from_file(filename)
    mesh_out = extract_edges(mesh, eps=float(options.eps))
    mesh_out = merge_lines(mesh_out)

    filename_out = options.filename_out
    if filename_out is None:
        filename_out = edit_filename(filename, prefix='edge_', new_ext='.vtk')

    output('Outline mesh - vertices: %d, edges: %d, output filename: %s'
           % (mesh_out[0].shape[0], mesh_out[2][0].shape[0], filename_out))

    # hack to write '3_2' elements - edges
    io = VTKMeshIO(None)
    aux_mesh = Struct()
    aux_mesh._get_io_data = lambda: mesh_out
    aux_mesh.n_el = mesh_out[2][0].shape[0]
    io.write(filename_out, aux_mesh)

    def __call__(self, state0=None, save_results=True, step_hook=None,
                 post_process_hook=None, nls_status=None):
        """
        Solve the time-dependent problem.
        """
        problem = self.problem
        ts = self.ts

        suffix, is_save = prepare_save_data(ts, problem.conf)

        if state0 is None:
            state0 = get_initial_state(problem)

        ii = 0
        for step, time in ts:
            output(self.format % (time, step + 1, ts.n_step))

            state = self.solve_step(ts, state0, nls_status=nls_status)
            state0 = state.copy(deep=True)

            if step_hook is not None:
                step_hook(problem, ts, state)

            if save_results and (is_save[ii] == ts.step):
                filename = problem.get_output_name(suffix=suffix % ts.step)
                problem.save_state(filename, state,
                                   post_process_hook=post_process_hook,
                                   file_per_var=None,
                                   ts=ts)
                ii += 1

            yield step, time, state

            problem.advance(ts)

def conv_test( conf, it, of, of0, ofg_norm = None ):
    """
    Returns
    -------
    flag : int
        * -1 ... continue
        *  0 ... small OF -> stop
        *  1 ... i_max reached -> stop
        *  2 ... small OFG -> stop
        *  3 ... small relative decrase of OF
     """

    status = -1
    output( 'opt: iter: %d, of: %e (||ofg||: %e)' % (it, of, ofg_norm) )
#    print (of0 - of), (conf.eps_rd * of0)

    if (abs( of ) < conf.eps_of):
        status = 0
    elif ofg_norm and (ofg_norm < conf.eps_ofg):
        status = 2
    elif (it > 0) and (abs(of0 - of) < (conf.eps_rd * abs( of0 ))):
        status = 3
        
    if (status == -1) and (it >= conf.i_max):
        status = 1

    return status

def setup_dof_conns(conn_info, dof_conns=None,
                    make_virtual=False, verbose=True):
    """
    Dof connectivity key:
        (field.name, var.n_components, region.name, type, ig)
    """
    if verbose:
        output('setting up dof connectivities...')
        tt = time.clock()

    dof_conns = get_default(dof_conns, {})

    for key, ii, info in iter_dict_of_lists(conn_info, return_keys=True):

        if info.primary is not None:
            var = info.primary
            field = var.get_field()
            field.setup_extra_data(info.ps_tg, info, info.is_trace)
            field.setup_dof_conns(dof_conns, var.n_components,
                                  info.dc_type, info.get_region(),
                                  info.is_trace)

        if info.has_virtual and not info.is_trace:
            # This is needed regardless make_virtual.
            var = info.virtual
            field = var.get_field()
            field.setup_extra_data(info.v_tg, info, False)
            field.setup_dof_conns(dof_conns, var.n_components,
                                  info.dc_type,
                                  info.get_region(can_trace=False))

    if verbose:
        output('...done in %.2f s' % (time.clock() - tt))

    return dof_conns

    def __call__(self, vec0=None, nls=None, init_fun=None, prestep_fun=None,
                 poststep_fun=None, status=None, **kwargs):
        """
        Solve elastodynamics problems by the Newmark method.
        """
        conf = self.conf
        nls = get_default(nls, self.nls)

        vec, unpack, pack = self.get_initial_vec(
            nls, vec0, init_fun, prestep_fun, poststep_fun)

        ts = self.ts
        for step, time in ts.iter_from(ts.step):
            output(self.format % (time, step + 1, ts.n_step),
                   verbose=self.verbose)
            dt = ts.dt

            prestep_fun(ts, vec)
            ut, vt, at = unpack(vec)

            nlst = self.create_nlst(nls, dt, conf.gamma, conf.beta, ut, vt, at)
            atp = nlst(at)
            vtp = nlst.v(atp)
            utp = nlst.u(atp)

            vect = pack(utp, vtp, atp)
            poststep_fun(ts, vect)

            vec = vect

        return vec

def make_implicit_step(ts, state0, problem, nls_status=None):
    """
    Make a step of an implicit time stepping solver.
    """
    if ts.step == 0:
        state0.apply_ebc()
        state = state0.copy(deep=True)

        if not ts.is_quasistatic:
            ev = problem.get_evaluator()
            try:
                vec_r = ev.eval_residual(state(), is_full=True)
            except ValueError:
                output('initial residual evaluation failed, giving up...')
                raise
            else:
                err = nm.linalg.norm(vec_r)
                output('initial residual: %e' % err)

        if problem.is_linear():
            mtx = prepare_matrix(problem, state)
            problem.try_presolve(mtx)

        # Initialize variables with history.
        state0.init_history()
        if ts.is_quasistatic:
            # Ordinary solve.
            state = problem.solve(state0=state0, nls_status=nls_status)

    else:
        problem.time_update(ts)
        state = problem.solve(state0=state0, nls_status=nls_status)

    return state

    def _check_region(self, region):
        """
        Check whether the `region` can be used for the
        field. Non-surface fields require the region to span whole
        element groups.

        Returns
        -------
        ok : bool
            True if the region is usable for the field.
        """
        ok = True
        domain = region.domain
        for ig in region.igs:
            if domain.groups[ig].gel.dim != domain.shape.tdim:
                output('cells with a bad topological dimension! (%d == %d)'
                       % (domain.groups[ig].gel.dim, domain.shape.tdim))
                ok = False
                break
            shape = domain.groups[ig].shape
            if region.shape[ig].n_vertex < shape.n_vertex:
                output('region does not span a whole element group!')
                ok = False
                break

        return ok

    def get_actual_order(self, geometry):
        """
        Return the actual integration order for given geometry.

        Parameters
        ----------
        geometry : str
            The geometry key describing the integration domain,
            see the keys of `sfepy.fem.quadratures.quadrature_tables`.

        Returns
        -------
        order : int
            If `self.order` is in quadrature tables it is this
            value. Otherwise it is the closest higher order. If no
            higher order is available, a warning is printed and the
            highest available order is used.
        """
        table = quadrature_tables[geometry]
        if self.order in table:
            order = self.order

        else:
            orders = table.keys()
            ii = nm.searchsorted(orders, self.order)
            if ii >= len(orders):
                order = max(orders)
                output(self._msg1 % (self.order, geometry))
                output(self._msg2 % order)

            else:
                order = orders[ii]

        return order

    def __init__(self, conf, **kwargs):
        LinearSolver.__init__(self, conf, solve=None, **kwargs)
        um = self.sls = None

        aux = try_imports(['import scipy.linsolve as sls',
                           'import scipy.splinalg.dsolve as sls',
                           'import scipy.sparse.linalg.dsolve as sls'],
                          'cannot import scipy sparse direct solvers!')
        self.sls = aux['sls']
        aux = try_imports(['import scipy.linsolve.umfpack as um',
                           'import scipy.splinalg.dsolve.umfpack as um',
                           'import scipy.sparse.linalg.dsolve.umfpack as um',
                           'import scikits.umfpack as um'])
        if 'um' in aux:
            um = aux['um']

        if um is not None:
            is_umfpack = hasattr(um, 'UMFPACK_OK')
        else:
            is_umfpack = False

        method = self.conf.method
        if method == 'superlu':
            self.sls.use_solver(useUmfpack=False)
        elif method == 'umfpack':
            if not is_umfpack and self.conf.warn:
                output('umfpack not available, using superlu!')
        elif method != 'auto':
            raise ValueError('uknown solution method! (%s)' % method)

        if method != 'superlu' and is_umfpack:
            self.sls.use_solver(useUmfpack=True,
                                assumeSortedIndices=True)

def post_process(out, pb, state, extend=False):
    """
    Calculate :math:`\nabla t` and compute boundary fluxes.
    """
    dv = pb.evaluate('ev_diffusion_velocity.i.Omega(m.K, t)', mode='el_avg',
                     verbose=False)
    out['dv'] = Struct(name='output_data', mode='cell',
                       data=dv, dofs=None)

    totals = nm.zeros(3)
    for gamma in ['Gamma_N', 'Gamma_N0', 'Gamma_D']:

        flux = pb.evaluate('d_surface_flux.i.%s(m.K, t)' % gamma,
                           verbose=False)
        area = pb.evaluate('d_surface.i.%s(t)' % gamma, verbose=False)

        flux_data = (gamma, flux, area, flux / area)
        totals += flux_data[1:]

        output('%8s flux: % 8.3f length: % 8.3f flux/length: % 8.3f'
               % flux_data)

    totals[2] = totals[0] / totals[1]
    output('   total flux: % 8.3f length: % 8.3f flux/length: % 8.3f'
           % tuple(totals))

    return out

def get_homog_mat(ts, coors, mode, term=None, problem=None, **kwargs):
    if problem.update_materials_flag == 2 and mode == 'qp':
        out = hyperelastic_data['homog_mat']
        return {k: nm.array(v) for k, v in six.iteritems(out)}
    elif problem.update_materials_flag == 0 or not mode == 'qp':
        return

    output('get_homog_mat')
    dim = problem.domain.mesh.dim

    update_var = problem.conf.options.mesh_update_variables[0]
    state_u = problem.equations.variables[update_var]
    state_u.field.clear_mappings()
    family_data = problem.family_data(state_u, term.region,
                                      term.integral, term.integration)

    mtx_f = family_data.mtx_f.reshape((coors.shape[0],)
                                      + family_data.mtx_f.shape[-2:])
    out = get_homog_coefs_nonlinear(ts, coors, mode, mtx_f,
                                    term=term, problem=problem,
                                    iteration=problem.iiter, **kwargs)

    out['E'] = 0.5 * (la.dot_sequences(mtx_f, mtx_f, 'ATB') - nm.eye(dim))

    hyperelastic_data['time'] = ts.step
    hyperelastic_data['homog_mat_shape'] = family_data.det_f.shape[:2]
    hyperelastic_data['homog_mat'] = \
        {k: nm.array(v) for k, v in six.iteritems(out)}

    return out

    def __call__( self, problem = None, data = None, save_hook = None ):
        """data: corrs_pressure, evp, optionally vec_g"""
        problem = get_default( problem, self.problem )
        ts = problem.get_time_solver().ts

        corrs, evp = [data[ii] for ii in self.requires[:2]]
        if len(self.requires) == 3:
            vec_g = data[self.requires[2]]
        else:
            vec_g = None

        assert_( evp.ebcs == self.ebcs )
        assert_( evp.epbcs == self.epbcs )

        filename = self.get_dump_name()
        savename = self.get_save_name()

        self.setup_equations(self.equations)

        solve = self.compute_correctors
        solve(evp, 1.0, corrs.state, ts, filename, savename, vec_g=vec_g)

        if self.check:
            self.setup_equations(self.verify_equations)
            self.init_solvers(problem)

            output( 'verifying correctors %s...' % self.name )
            verify = self.verify_correctors
            ok = verify(1.0, corrs.state, filename)
            output( '...done, ok: %s' % ok )

        return Struct( name = self.name,
                       filename = filename )

def get_std_wave_fun(pb, options):
    stiffness = pb.evaluate('ev_volume_integrate_mat.2.Omega(m.D, u)',
                            mode='el_avg', copy_materials=False, verbose=False)
    young, poisson = mc.youngpoisson_from_stiffness(stiffness,
                                                    plane=options.plane)
    density = pb.evaluate('ev_volume_integrate_mat.2.Omega(m.density, u)',
                          mode='el_avg', copy_materials=False, verbose=False)

    lam, mu = mc.lame_from_youngpoisson(young, poisson,
                                        plane=options.plane)
    alam = nm.average(lam)
    amu = nm.average(mu)
    adensity = nm.average(density)

    cp = nm.sqrt((alam + 2.0 * amu) / adensity)
    cs = nm.sqrt(amu / adensity)
    output('average p-wave speed:', cp)
    output('average shear wave speed:', cs)

    log_names = [r'$\omega_p$', r'$\omega_s$']
    log_plot_kwargs = [{'ls' : '--', 'color' : 'k'},
                       {'ls' : '--', 'color' : 'gray'}]

    if options.mode == 'omega':
        fun = lambda wmag, wdir: (cp * wmag, cs * wmag)

    else:
        fun = lambda wmag, wdir: (wmag / cp, wmag / cs)

    return fun, log_names, log_plot_kwargs