def assemble(mtx_d):
    m = Material('m', D=mtx_d, rho=density)

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

    t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=m, v=v, u=u)
    eq1 = Equation('stiffness', t1)
    eq2 = Equation('mass', t2)
    lhs_eqs = Equations([eq1, eq2])

    pb = Problem('modal', equations=lhs_eqs)

    pb.time_update()
    n_rbm = dim * (dim + 1) / 2

    pb.update_materials()

    tmp = time.time()
    # Assemble stiffness and mass matrices.
    mtx_k = eq1.evaluate(mode='weak', dw_mode='matrix', asm_obj=pb.mtx_a)
    mtx_m = mtx_k.copy()
    mtx_m.data[:] = 0.0
    mtx_m = eq2.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx_m)

    return mtx_k, mtx_m

def evalValAndDeriv(D):
    m = Material('m', D = D, rho = 2700.0)

    integral = Integral('i', order=2)

    t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_volume_dot(m.rho, v, u)', integral, omega, m=m, v=v, u=u)
    eq1 = Equation('stiffness', t1)
    eq2 = Equation('mass', t2)
    lhs_eqs = Equations([eq1, eq2])

    pb = Problem('modal', equations = lhs_eqs)

    pb.time_update()
    n_rbm = dim * (dim + 1) / 2

    pb.update_materials()

    # Assemble stiffness and mass matrices.
    mtx_k = eq1.evaluate(mode='weak', dw_mode='matrix', asm_obj=pb.mtx_a)
    mtx_m = mtx_k.copy()
    mtx_m.data[:] = 0.0
    mtx_m = eq2.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx_m)

    eigs0, evecs0 = scipy.sparse.linalg.eigsh(mtx_k, k = 10, M = mtx_m, which = 'SM')

    eigs = eigs0[3:]
    evecs = evecs0[:, 3:]

    dydmu = numpy.array([evecs[:, i].T.dot(dKdmu.dot(evecs[:, i])) for i in range(evecs.shape[1])])
    dydlambda = numpy.array([evecs[:, i].T.dot(dKdlambda.dot(evecs[:, i])) for i in range(evecs.shape[1])])

    return eigs, dydmu, dydlambda

def solve_problem(shape, dims, young, poisson, force, transform=None):
    domain = make_domain(dims[:2], shape, transform=transform)

    omega = domain.regions['Omega']
    gamma1 = domain.regions['Gamma1']
    gamma2 = domain.regions['Gamma2']

    field = Field.from_args('fu', nm.float64, 6, omega, approx_order=1,
                            poly_space_base='shell10x')
    u = FieldVariable('u', 'unknown', field)
    v = FieldVariable('v', 'test', field, primary_var_name='u')

    thickness = dims[2]
    if transform is None:
        pload = [[0.0, 0.0, force / shape[1], 0.0, 0.0, 0.0]] * shape[1]

    elif transform == 'bend':
        pload = [[force / shape[1], 0.0, 0.0, 0.0, 0.0, 0.0]] * shape[1]

    elif transform == 'twist':
        pload = [[0.0, force / shape[1], 0.0, 0.0, 0.0, 0.0]] * shape[1]

    m = Material('m', D=sh.create_elastic_tensor(young=young, poisson=poisson),
                 values={'.drill' : 1e-7})
    load = Material('load', values={'.val' : pload})

    aux = Integral('i', order=3)
    qp_coors, qp_weights = aux.get_qp('3_8')
    qp_coors[:, 2] = thickness * (qp_coors[:, 2] - 0.5)
    qp_weights *= thickness

    integral = Integral('i', coors=qp_coors, weights=qp_weights, order='custom')

    t1 = Term.new('dw_shell10x(m.D, m.drill, v, u)',
                  integral, omega, m=m, v=v, u=u)
    t2 = Term.new('dw_point_load(load.val, v)',
                  integral, gamma2, load=load, v=v)
    eq = Equation('balance', t1 - t2)
    eqs = Equations([eq])

    fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem('elasticity with shell10x', equations=eqs, nls=nls, ls=ls)
    pb.time_update(ebcs=Conditions([fix_u]))

    state = pb.solve()

    return pb, state, u, gamma2

    def test_solving(self):
        from sfepy.base.base import IndexedStruct
        from sfepy.discrete import (FieldVariable, Material, Problem, Function,
                                    Equation, Equations, Integral)
        from sfepy.discrete.conditions import Conditions, EssentialBC
        from sfepy.terms import Term
        from sfepy.solvers.ls import ScipyDirect
        from sfepy.solvers.nls import Newton
        from sfepy.mechanics.matcoefs import stiffness_from_lame

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

        m = Material('m', D=stiffness_from_lame(self.dim, 1.0, 1.0))
        f = Material('f', val=[[0.02], [0.01]])

        bc_fun = Function('fix_u_fun', fix_u_fun,
                          extra_args={'extra_arg' : 'hello'})

        fix_u = EssentialBC('fix_u', self.gamma1, {'u.all' : bc_fun})
        shift_u = EssentialBC('shift_u', self.gamma2, {'u.0' : 0.1})

        integral = Integral('i', order=3)

        t1 = Term.new('dw_lin_elastic(m.D, v, u)',
                      integral, self.omega, m=m, v=v, u=u)

        t2 = Term.new('dw_volume_lvf(f.val, v)', integral, self.omega, f=f, v=v)

        eq = Equation('balance', t1 + t2)
        eqs = Equations([eq])

        ls = ScipyDirect({})

        nls_status = IndexedStruct()
        nls = Newton({}, lin_solver=ls, status=nls_status)

        pb = Problem('elasticity', equations=eqs)
        ## pb.save_regions_as_groups('regions')

        pb.set_bcs(ebcs=Conditions([fix_u, shift_u]))
        pb.set_solver(nls)

        state = pb.solve()

        name = op.join(self.options.out_dir, 'test_high_level_solving.vtk')
        pb.save_state(name, state)

        ok = nls_status.condition == 0
        if not ok:
            self.report('solver did not converge!')

        _ok = state.has_ebc()
        if not _ok:
            self.report('EBCs violated!')

        ok = ok and _ok

        return ok

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 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']})
    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())
    problem.solve()

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

    def from_conf(conf, options):
        from sfepy.discrete import Problem

        problem = Problem.from_conf(conf, init_equations=False)
        test = Test(problem=problem,
                    conf=conf, options=options)
        return test

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

    parser = OptionParser(usage=usage, version='%prog')
    parser.add_option('-n', '--no-plot',
                      action="store_true", dest='no_plot',
                      default=False, help=helps['no_plot'])
    options, args = parser.parse_args()

    required, other = get_standard_keywords()
    # Use this file as the input file.
    conf = ProblemConf.from_file(__file__, required, other)

    # Create problem instance, but do not set equations.
    problem = Problem.from_conf(conf, init_equations=False)

    # Solve the problem. Output is ignored, results stored by using the
    # step_hook.
    u_t = solve_branch(problem, linear_tension)
    u_c = solve_branch(problem, linear_compression)

    # Get pressure load by calling linear_*() for each time step.
    ts = problem.get_timestepper()
    load_t = nm.array([linear_tension(ts, nm.array([[0.0]]), 'qp')['val']
                       for aux in ts.iter_from(0)],
                      dtype=nm.float64).squeeze()
    load_c = nm.array([linear_compression(ts, nm.array([[0.0]]), 'qp')['val']
                       for aux in ts.iter_from(0)],
                      dtype=nm.float64).squeeze()

    # Join the branches.
    displacements = {}
    for key in u_t.keys():
        displacements[key] = nm.r_[u_c[key][::-1], u_t[key]]
    load = nm.r_[load_c[::-1], load_t]


    if plt is None:
        output('matplotlib cannot be imported, printing raw data!')
        output(displacements)
        output(load)
    else:
        legend = []
        for key, val in six.iteritems(displacements):
            plt.plot(load, val)
            legend.append(key)

        plt.legend(legend, loc = 2)
        plt.xlabel('tension [kPa]')
        plt.ylabel('displacement [mm]')
        plt.grid(True)

        plt.gcf().savefig('pressure_displacement.png')

        if not options.no_plot:
            plt.show()

def vary_omega1_size( problem ):
    """Vary size of \Omega1. Saves also the regions into options['output_dir'].

    Input:
      problem: Problem instance
    Return:
      a generator object:
      1. creates new (modified) problem
      2. yields the new (modified) problem and output container
      3. use the output container for some logging
      4. yields None (to signal next iteration to Application)
    """
    from sfepy.discrete import Problem
    from sfepy.solvers.ts import get_print_info
    
    output.prefix = 'vary_omega1_size:'

    diameters = nm.linspace( 0.1, 0.6, 7 ) + 0.001
    ofn_trunk, output_format = problem.ofn_trunk, problem.output_format
    output_dir = problem.output_dir
    join = os.path.join

    conf = problem.conf
    cf = conf.get_raw( 'functions' )
    n_digit, aux, d_format = get_print_info( len( diameters ) + 1 )
    for ii, diameter in enumerate( diameters ):
        output( 'iteration %d: diameter %3.2f' % (ii, diameter) )

        cf['select_circ'] = (lambda coors, domain=None: 
                             select_circ(coors[:,0], coors[:,1], 0, diameter),)
        conf.edit('functions', cf)
        problem = Problem.from_conf(conf)

        problem.save_regions( join( output_dir, ('regions_' + d_format) % ii ),
                              ['Omega_1'] )
        region = problem.domain.regions['Omega_1']
        if not region.has_cells():
            raise ValueError('region %s has no cells!' % region.name)

        ofn_trunk = ofn_trunk + '_' + (d_format % ii)
        problem.setup_output(output_filename_trunk=ofn_trunk,
                             output_dir=output_dir,
                             output_format=output_format)

        out = []
        yield problem, out

        out_problem, state = out[-1]

        filename = join( output_dir,
                         ('log_%s.txt' % d_format) % ii )
        fd = open( filename, 'w' )
        log_item = '$r(\Omega_1)$: %f\n' % diameter
        fd.write( log_item )
        fd.write( 'solution:\n' )
        nm.savetxt(fd, state())
        fd.close()

        yield None

    def from_conf(conf, options):
        from sfepy.discrete import Problem

        problem = Problem.from_conf(conf)
        problem.time_update()

        test = Test(problem=problem, conf=conf, options=options)
        return test

def assemble_matrices(define, mod, pars, set_wave_dir, options):
    """
    Assemble the blocks of dispersion eigenvalue problem matrices.
    """
    define_problem = functools.partial(define,
                                       filename_mesh=options.mesh_filename,
                                       pars=pars,
                                       approx_order=options.order,
                                       refinement_level=options.refine,
                                       solver_conf=options.solver_conf,
                                       plane=options.plane,
                                       post_process=options.post_process)

    conf = ProblemConf.from_dict(define_problem(), mod)

    pb = Problem.from_conf(conf)
    pb.dispersion_options = options
    pb.set_output_dir(options.output_dir)
    dim = pb.domain.shape.dim

    # Set the normalized wave vector direction to the material(s).
    wdir = nm.asarray(options.wave_dir[:dim], dtype=nm.float64)
    wdir = wdir / nm.linalg.norm(wdir)
    set_wave_dir(pb, wdir)

    bbox = pb.domain.mesh.get_bounding_box()
    size = (bbox[1] - bbox[0]).max()
    scaling0 = apply_unit_multipliers([1.0], ['length'],
                                      options.unit_multipliers)[0]
    scaling = scaling0
    if options.mesh_size is not None:
        scaling *= options.mesh_size / size
    output('scaling factor of periodic cell mesh coordinates:', scaling)
    output('new mesh size with applied unit multipliers:', scaling * size)
    pb.domain.mesh.coors[:] *= scaling
    pb.set_mesh_coors(pb.domain.mesh.coors, update_fields=True)

    bzone = 2.0 * nm.pi / (scaling * size)
    output('1. Brillouin zone size:', bzone * scaling0)
    output('1. Brillouin zone size with applied unit multipliers:', bzone)

    pb.time_update()
    pb.update_materials()

    # Assemble the matrices.
    mtxs = {}
    for key, eq in pb.equations.iteritems():
        mtxs[key] = mtx = pb.mtx_a.copy()
        mtx = eq.evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx)
        mtx.eliminate_zeros()
        output_array_stats(mtx.data, 'nonzeros in %s' % key)

        output('symmetry checks:')
        output('%s - %s^T:' % (key, key), max_diff_csr(mtx, mtx.T))
        output('%s - %s^H:' % (key, key), max_diff_csr(mtx, mtx.H))

    return pb, wdir, bzone, mtxs

    def test_stokes_slip_bc(self):
        import scipy.sparse as sp

        from sfepy.base.conf import ProblemConf
        from sfepy.discrete import Problem
        import examples.navier_stokes.stokes_slip_bc as ssb

        conf = ProblemConf.from_module(ssb)
        pb = Problem.from_conf(conf, init_solvers=False)
        pb.time_update()
        variables = pb.get_variables()

        adi = variables.adi
        lcdi = variables.lcdi
        mtx = variables.mtx_lcbc

        ok = adi.var_names == lcdi.var_names
        self.report('same adi-lcdi ordering:', ok)

        ublock = mtx[adi.indx['u']]
        ir, ic = ublock.nonzero()
        ir += adi.indx['u'].start

        i0, i1 = adi.indx['u'].start, adi.indx['u'].stop
        _ok0 = (i0 <= ir).all() and  (ir < i1).all()
        self.report('u block rows in [%d %d[: %s' % (i0, i1, _ok0))

        i0, i1 = lcdi.indx['u'].start, lcdi.indx['u'].stop
        _ok1 = (i0 <= ic).all() and  (ic < i1).all()
        self.report('u block cols in [%d %d[: %s' % (i0, i1, _ok1))

        ok = ok and _ok0 and _ok1

        pblock = mtx[adi.indx['p']]
        ir, ic, iv = sp.find(pblock)
        ir += adi.indx['p'].start

        i0, i1 = adi.indx['p'].start, adi.indx['p'].stop
        _ok0 = (i0 <= ir).all() and  (ir < i1).all()
        self.report('p block rows in [%d %d[: %s' % (i0, i1, _ok0))

        i0, i1 = lcdi.indx['p'].start, lcdi.indx['p'].stop
        _ok1 = (i0 <= ic).all() and  (ic < i1).all()
        self.report('p block cols in [%d %d[: %s' % (i0, i1, _ok1))

        ok = ok and _ok0 and _ok1

        _ok0 = (len(ir) == adi.n_dof['p'])
        self.report('p block size correct:', _ok0)
        _ok1 = ((ir - adi.indx['p'].start) == (ic - lcdi.indx['p'].start)).all()
        self.report('p block diagonal:', _ok1)
        _ok2 = (iv == 1.0).all()
        self.report('p block identity:', _ok2)

        ok = ok and _ok0 and _ok1 and _ok2

        return ok

def make_h1_projection_data(target, eval_data):
    """
    Project scalar data given by a material-like `eval_data()` function to a
    scalar `target` field variable using the :math:`H^1` dot product.
    """
    order = target.field.approx_order * 2
    integral = Integral('i', order=order)

    un = target.name
    v = FieldVariable('v', 'test', target.field, primary_var_name=un)
    lhs1 = Term.new('dw_volume_dot(v, %s)' % un, integral,
                    target.field.region, v=v, **{un : target})
    lhs2 = Term.new('dw_laplace(v, %s)' % un, integral,
                    target.field.region, v=v, **{un : target})

    def _eval_data(ts, coors, mode, **kwargs):
        if mode == 'qp':
            val = eval_data(ts, coors, mode, 'val', **kwargs)
            gval = eval_data(ts, coors, mode, 'grad', **kwargs)
            return {'val' : val, 'gval' : gval}

    m = Material('m', function=_eval_data)
    rhs1 = Term.new('dw_volume_lvf(m.val, v)', integral, target.field.region,
                    m=m, v=v)
    rhs2 = Term.new('dw_diffusion_r(m.gval, v)', integral, target.field.region,
                    m=m, v=v)

    eq = Equation('projection', lhs1 + lhs2 - rhs1 - rhs2)
    eqs = Equations([eq])

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

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

    pb.time_update()

    # This sets the target variable with the projection solution.
    pb.solve()

    if nls_status.condition != 0:
        output('H1 projection: solver did not converge!')

def main():
    from sfepy import data_dir

    parser = OptionParser(usage=usage, version="%prog")
    parser.add_option("-s", "--show", action="store_true", dest="show", default=False, help=help["show"])
    options, args = parser.parse_args()

    mesh = Mesh.from_file(data_dir + "/meshes/2d/rectangle_tri.mesh")
    domain = Domain("domain", mesh)

    min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
    eps = 1e-8 * (max_x - min_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, "vector", omega, approx_order=2)

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

    m = Material("m", lam=1.0, mu=1.0)
    f = Material("f", val=[[0.02], [0.01]])

    integral = Integral("i", order=3)

    t1 = Term.new("dw_lin_elastic_iso(m.lam, m.mu, v, u)", integral, omega, m=m, v=v, u=u)
    t2 = Term.new("dw_volume_lvf(f.val, v)", integral, omega, f=f, v=v)
    eq = Equation("balance", t1 + t2)
    eqs = Equations([eq])

    fix_u = EssentialBC("fix_u", gamma1, {"u.all": 0.0})

    bc_fun = Function("shift_u_fun", shift_u_fun, extra_args={"shift": 0.01})
    shift_u = EssentialBC("shift_u", gamma2, {"u.0": bc_fun})

    ls = ScipyDirect({})

    nls_status = IndexedStruct()
    nls = Newton({}, lin_solver=ls, status=nls_status)

    pb = Problem("elasticity", equations=eqs, nls=nls, ls=ls)
    pb.save_regions_as_groups("regions")

    pb.time_update(ebcs=Conditions([fix_u, shift_u]))

    vec = pb.solve()
    print nls_status

    pb.save_state("linear_elasticity.vtk", vec)

    if options.show:
        view = Viewer("linear_elasticity.vtk")
        view(vector_mode="warp_norm", rel_scaling=2, is_scalar_bar=True, is_wireframe=True)

def run(domain, order):
    omega = domain.create_region('Omega', 'all')
    bbox = domain.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    min_y, max_y = bbox[:, 1]
    eps = 1e-8 * (max_x - min_x)
    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')
    gamma3 = domain.create_region('Gamma3',
                                  'vertices in y < %.10f' % (min_y + eps),
                                  'facet')
    gamma4 = domain.create_region('Gamma4',
                                  'vertices in y > %.10f' % (max_y - eps),
                                  'facet')

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

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

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

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

    fix1 = EssentialBC('fix1', gamma1, {'u.0' : 0.4})
    fix2 = EssentialBC('fix2', gamma2, {'u.0' : 0.0})

    def get_shift(ts, coors, region):
        return nm.ones_like(coors[:, 0])

    dof_map_fun = Function('dof_map_fun', per.match_x_line)
    shift_fun = Function('shift_fun', get_shift)

    sper = LinearCombinationBC('sper', [gamma3, gamma4], {'u.0' : 'u.0'},
                               dof_map_fun, 'shifted_periodic',
                               arguments=(shift_fun,))

    ls = ScipyDirect({})

    pb = Problem('laplace', equations=eqs, auto_solvers=None)

    pb.time_update(ebcs=Conditions([fix1, fix2]), lcbcs=Conditions([sper]))

    ev = pb.get_evaluator()
    nls = Newton({}, lin_solver=ls,
                 fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

    pb.set_solver(nls)

    state = pb.solve()

    return pb, state

    def from_conf(conf, options):
        from sfepy.discrete import Problem

        pb = Problem.from_conf(conf, init_solvers=False)
        pb.time_update()
        mtx = pb.equations.evaluate(mode='weak', dw_mode='matrix',
                                    asm_obj=pb.mtx_a)

        test = Test(mtx=mtx, conf=conf, options=options)
        return test

    def __init__(self, conf, options, output_prefix, init_equations=True, **kwargs):
        """`kwargs` are passed to  Problem.from_conf()

        Command-line options have precedence over conf.options."""
        Application.__init__(self, conf, options, output_prefix)
        self.setup_options()

        is_eqs = init_equations
        if hasattr(options, "solve_not") and options.solve_not:
            is_eqs = False
        self.problem = Problem.from_conf(conf, init_equations=is_eqs, **kwargs)

        self.setup_output_info(self.problem, self.options)

def run(domain, order):
    omega = domain.create_region("Omega", "all")
    bbox = domain.get_mesh_bounding_box()
    min_x, max_x = bbox[:, 0]
    min_y, max_y = bbox[:, 1]
    eps = 1e-8 * (max_x - min_x)
    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")
    gamma3 = domain.create_region("Gamma3", "vertices in y < %.10f" % (min_y + eps), "facet")
    gamma4 = domain.create_region("Gamma4", "vertices in y > %.10f" % (max_y - eps), "facet")

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

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

    integral = Integral("i", order=2 * order)

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

    fix1 = EssentialBC("fix1", gamma1, {"u.0": 0.4})
    fix2 = EssentialBC("fix2", gamma2, {"u.0": 0.0})

    def get_shift(ts, coors, region):
        return nm.ones_like(coors[:, 0])

    dof_map_fun = Function("dof_map_fun", per.match_x_line)
    shift_fun = Function("shift_fun", get_shift)

    sper = LinearCombinationBC(
        "sper", [gamma3, gamma4], {"u.0": "u.0"}, dof_map_fun, "shifted_periodic", arguments=(shift_fun,)
    )

    ls = ScipyDirect({})

    pb = Problem("laplace", equations=eqs, auto_solvers=None)

    pb.time_update(ebcs=Conditions([fix1, fix2]), lcbcs=Conditions([sper]))

    ev = pb.get_evaluator()
    nls = Newton({}, lin_solver=ls, fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix)

    pb.set_solver(nls)

    state = pb.solve()

    return pb, state

def recover_micro_hook( micro_filename, region, macro,
                        naming_scheme = 'step_iel',
                        recovery_file_tag='' ):

    # Create a micro-problem instance.
    required, other = get_standard_keywords()
    required.remove( 'equations' )
    pb = Problem.from_conf_file(micro_filename, required=required, other=other,
                                init_equations=False, init_solvers=False)

    coefs_filename = pb.conf.options.get('coefs_filename', 'coefs')
    output_dir = pb.conf.options.get('output_dir', '.')
    coefs_filename = op.join(output_dir, coefs_filename) + '.h5'

    # Coefficients and correctors
    coefs = Coefficients.from_file_hdf5( coefs_filename )
    corrs = get_correctors_from_file( dump_names = coefs.dump_names ) 

    recovery_hook = pb.conf.options.get('recovery_hook', None)

    if recovery_hook is not None:
        recovery_hook = pb.conf.get_function(recovery_hook)

        aux = max(pb.domain.shape.n_gr, 2)
        format = get_print_info( aux, fill = '0' )[1] \
            + '_' + get_print_info( pb.domain.mesh.n_el, fill = '0' )[1]

        for ig, ii, iel in region.iter_cells():
            print 'ig: %d, ii: %d, iel: %d' % (ig, ii, iel)

            local_macro = {}
            for k, v in macro.iteritems():
                local_macro[k] = v[ii,0]

            out = recovery_hook( pb, corrs, local_macro )

            # save data
            suffix = format % (ig, iel)
            micro_name = pb.get_output_name(extra='recovered_'\
                                            + recovery_file_tag + suffix)
            filename = op.join(output_dir, op.basename(micro_name))
            fpv = pb.conf.options.get('file_per_var', False)
            pb.save_state(filename, out=out,
                          file_per_var=fpv)

    def test_save_ebc(self):
        from sfepy.discrete import (FieldVariable, Integral,
                                    Equation, Equations, Problem)
        from sfepy.discrete.conditions import Conditions, EssentialBC
        from sfepy.terms import Term

        name = op.join(self.options.out_dir,
                       op.splitext(op.basename(__file__))[0])

        integral = Integral('i', order=1)

        u = self.variables['u']
        v = FieldVariable('v', 'test', u.field, primary_var_name='u')

        p = self.variables['p']
        q = FieldVariable('q', 'test', p.field, primary_var_name='p')

        regions = self.problem.domain.regions
        omega = regions['Omega']

        # Problem.save_ebc() requires to have equations defined.
        t1 = Term.new('dw_lin_elastic(v, u)',
                      integral, omega, v=v, u=u)
        t2 = Term.new('dw_laplace(q, p)', integral, omega, q=q, p=p)
        eq = Equation('aux', t1 + t2)
        eqs = Equations([eq])

        pb = Problem('test', equations=eqs, auto_solvers=False)

        all_ebcs = []
        all_ebcs.append(EssentialBC('fix_u1', regions['RightFix'],
                                    {'u.all' : nm.array([0.0, 1.0])}))
        all_ebcs.append(EssentialBC('fix_u2', regions['LeftStrip'],
                                    {'u.0' : 0.0, 'u.1' : 1.0}))
        all_ebcs.append(EssentialBC('fix_p1', regions['LeftFix'],
                                    {'p.all' : 0.0}))
        all_ebcs.append(EssentialBC('fix_p2', regions['RightStrip'],
                                    {'p.0' : 0.0}))
        ebcs = Conditions(all_ebcs)

        pb.time_update(ebcs=ebcs)

        pb.save_ebc(name + '_ebcs_f.vtk', ebcs=ebcs, force=True)
        pb.save_ebc(name + '_ebcs.vtk', ebcs=ebcs, default=-1, force=False)

        return True