def nodal_stress(out, pb, state, extend=False, integrals=None):
    '''
    Calculate stresses at nodal points.
    '''

    # Point load.
    mat = pb.get_materials()['Load']
    P = 2.0 * mat.get_data('special', 'val')[1]

    # Calculate nodal stress.
    pb.time_update()

    if integrals is None: integrals = pb.get_integrals()

    stress = pb.evaluate('ev_cauchy_stress.ivn.Omega(Asphalt.D, u)', mode='qp',
                         integrals=integrals, copy_materials=False)
    sfield = Field.from_args('stress', nm.float64, (3,),
                             pb.domain.regions['Omega'])
    svar = FieldVariable('sigma', 'parameter', sfield,
                         primary_var_name='(set-to-None)')
    svar.set_from_qp(stress, integrals['ivn'])

    print('\n==================================================================')
    print('Given load = %.2f N' % -P)
    print('\nAnalytical solution')
    print('===================')
    print('Horizontal tensile stress = %.5e MPa/mm' % (-2.*P/(nm.pi*150.)))
    print('Vertical compressive stress = %.5e MPa/mm' % (-6.*P/(nm.pi*150.)))
    print('\nFEM solution')
    print('============')
    print('Horizontal tensile stress = %.5e MPa/mm' % (svar()[0]))
    print('Vertical compressive stress = %.5e MPa/mm' % (-svar()[1]))
    print('==================================================================')
    return out

def nodal_stress(out, pb, state, extend=False, integrals=None):
    """
    Calculate stresses at nodal points.
    """

    # Point load.
    mat = pb.get_materials()["Load"]
    P = 2.0 * mat.get_data("special", "val")[1]

    # Calculate nodal stress.
    pb.time_update()

    if integrals is None:
        integrals = pb.get_integrals()

    stress = pb.evaluate("ev_cauchy_stress.ivn.Omega(Asphalt.D, u)", mode="qp", integrals=integrals)
    sfield = Field.from_args("stress", nm.float64, (3,), pb.domain.regions["Omega"])
    svar = FieldVariable("sigma", "parameter", sfield, primary_var_name="(set-to-None)")
    svar.set_data_from_qp(stress, integrals["ivn"])

    print "\n=================================================================="
    print "Given load = %.2f N" % -P
    print "\nAnalytical solution"
    print "==================="
    print "Horizontal tensile stress = %.5e MPa/mm" % (-2.0 * P / (nm.pi * 150.0))
    print "Vertical compressive stress = %.5e MPa/mm" % (-6.0 * P / (nm.pi * 150.0))
    print "\nFEM solution"
    print "============"
    print "Horizontal tensile stress = %.5e MPa/mm" % (svar()[0][0])
    print "Vertical compressive stress = %.5e MPa/mm" % (-svar()[0][1])
    print "=================================================================="
    return out

    def test_volume_tl(self):
        from sfepy.discrete import FieldVariable

        fu = self.problem.fields['vector']
        fq = self.problem.fields['scalar']

        var_u = FieldVariable('u', 'parameter', fu,
                              primary_var_name='(set-to-None)')
        var_q = FieldVariable('q', 'test', fq,
                              primary_var_name='(set-to-None)')

        var_u.set_data(nm.linspace(0, 0.004, var_u.n_dof))

        vval = self.problem.evaluate('dw_tl_volume.i.Omega( q, u )',
                                     term_mode='volume', q=var_q, u=var_u)

        sval = self.problem.evaluate('d_tl_volume_surface.i.Gamma( u )',
                                     u=var_u)

        ok = abs(vval - sval) < 1e-14

        self.report('TL: by volume: %e == by surface: %e -> %s' %
                    (vval, sval, ok))

        return ok

    def test_surface_evaluate(self):
        from sfepy.discrete import FieldVariable
        problem = self.problem

        us = problem.get_variables()['us']
        vec = nm.empty(us.n_dof, dtype=us.dtype)
        vec[:] = 1.0
        us.set_data(vec)

        expr = 'ev_surface_integrate.i.Left( us )'
        val = problem.evaluate(expr, us=us)
        ok1 = nm.abs(val - 1.0) < 1e-15
        self.report('with unknown: %s, value: %s, ok: %s'
                    % (expr, val, ok1))

        ps1 = FieldVariable('ps1', 'parameter', us.get_field(),
                            primary_var_name='(set-to-None)')
        ps1.set_data(vec)

        expr = 'ev_surface_integrate.i.Left( ps1 )'
        val = problem.evaluate(expr, ps1=ps1)
        ok2 = nm.abs(val - 1.0) < 1e-15
        self.report('with parameter: %s, value: %s, ok: %s'
                    % (expr, val, ok2))
        ok2 = True

        return ok1 and ok2

    def test_projection_iga_fem(self):
        from sfepy.discrete import FieldVariable
        from sfepy.discrete.fem import FEDomain, Field
        from sfepy.discrete.iga.domain import IGDomain
        from sfepy.mesh.mesh_generators import gen_block_mesh
        from sfepy.discrete.iga.domain_generators import gen_patch_block_domain
        from sfepy.discrete.projections import (make_l2_projection,
                                                make_l2_projection_data)

        shape = [10, 12, 12]
        dims = [5, 6, 6]
        centre = [0, 0, 0]
        degrees = [2, 2, 2]

        nurbs, bmesh, regions = gen_patch_block_domain(dims, shape, centre,
                                                       degrees,
                                                       cp_mode='greville',
                                                       name='iga')
        ig_domain = IGDomain('iga', nurbs, bmesh, regions=regions)

        ig_omega = ig_domain.create_region('Omega', 'all')
        ig_field = Field.from_args('iga', nm.float64, 1, ig_omega,
                                   approx_order='iga', poly_space_base='iga')
        ig_u = FieldVariable('ig_u', 'parameter', ig_field,
                             primary_var_name='(set-to-None)')

        mesh = gen_block_mesh(dims, shape, centre, name='fem')
        fe_domain = FEDomain('fem', mesh)

        fe_omega = fe_domain.create_region('Omega', 'all')
        fe_field = Field.from_args('fem', nm.float64, 1, fe_omega,
                                   approx_order=2)
        fe_u = FieldVariable('fe_u', 'parameter', fe_field,
                             primary_var_name='(set-to-None)')

        def _eval_data(ts, coors, mode, **kwargs):
            return nm.prod(coors**2, axis=1)[:, None, None]

        make_l2_projection_data(ig_u, _eval_data)

        make_l2_projection(fe_u, ig_u) # This calls ig_u.evaluate_at().

        coors = 0.5 * nm.random.rand(20, 3) * dims

        ig_vals = ig_u.evaluate_at(coors)
        fe_vals = fe_u.evaluate_at(coors)

        ok = nm.allclose(ig_vals, fe_vals, rtol=0.0, atol=1e-12)
        if not ok:
            self.report('iga-fem projection failed!')
            self.report('coors:')
            self.report(coors)
            self.report('iga fem diff:')
            self.report(nm.c_[ig_vals, fe_vals, nm.abs(ig_vals - fe_vals)])

        return ok

    def create_subequations(self, var_names, known_var_names=None):
        """
        Create sub-equations containing only terms with the given virtual
        variables.

        Parameters
        ----------
        var_names : list
            The list of names of virtual variables.
        known_var_names : list
            The list of  names of (already) known state variables.

        Returns
        -------
        subequations : Equations instance
            The sub-equations.
        """
        from sfepy.discrete import FieldVariable

        known_var_names = get_default(known_var_names, [])

        objs = []
        for iv, var_name in enumerate(var_names):
            terms = [term.copy(name=term.name)
                     for eq in self for term in eq.terms
                     if term.get_virtual_name() == var_name]

            # Make parameter variables from known state variables in terms
            # arguments.
            for known_name in known_var_names:
                for term in terms:
                    if known_name in term.arg_names:
                        ii = term.arg_names.index(known_name)
                        state = self.variables[known_name]
                        par = FieldVariable(known_name, 'parameter',
                                            state.field,
                                            primary_var_name='(set-to-None)')
                        term.args[ii] = par
                        term._kwargs[known_name] = par
                        par.set_data(state())

            new_terms = Terms(terms)
            objs.append(Equation('eq_%d' % iv, new_terms))

        subequations = Equations(objs)

        return subequations

    def test_variables(self):
        from sfepy.discrete import FieldVariable

        u = FieldVariable('u', 'parameter', self.field,
                          primary_var_name='(set-to-None)')

        u.set_constant(1.0)

        vec = u() # Nodal values.

        ok = nm.allclose(vec, 1.0)

        ## print u() 
        ## print u.get_vector() # Coefficient vector w.r.t. the field space basis.
        ## print u(gamma1)
        ## print u.get_vector(gamma2)

        return ok

def prepare_variable(filename, n_components):
    from sfepy.discrete import FieldVariable
    from sfepy.discrete.fem import Mesh, FEDomain, Field

    mesh = Mesh.from_file(filename)

    bbox = mesh.get_bounding_box()
    dd = bbox[1, :] - bbox[0, :]
    data = nm.sin(4.0 * nm.pi * mesh.coors[:, 0:1] / dd[0]) * nm.cos(4.0 * nm.pi * mesh.coors[:, 1:2] / dd[1])

    domain = FEDomain("domain", mesh)
    omega = domain.create_region("Omega", "all")
    field = Field.from_args("field", nm.float64, n_components, omega, approx_order=2)

    u = FieldVariable("u", "parameter", field, primary_var_name="(set-to-None)")
    u.set_from_mesh_vertices(data * nm.arange(1, n_components + 1)[None, :])

    return u

    def test_projection_tri_quad(self):
        from sfepy.discrete.projections import make_l2_projection

        source = FieldVariable('us', 'unknown', self.field)

        coors = self.field.get_coor()
        vals = nm.sin(2.0 * nm.pi * coors[:,0] * coors[:,1])
        source.set_data(vals)

        name = op.join(self.options.out_dir,
                       'test_projection_tri_quad_source.vtk')
        source.save_as_mesh(name)

        mesh = Mesh.from_file('meshes/2d/square_quad.mesh',
                              prefix_dir=sfepy.data_dir)
        domain = FEDomain('domain', mesh)

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


        field = Field.from_args('bilinear', nm.float64, 'scalar', omega,
                                approx_order=1)

        target = FieldVariable('ut', 'unknown', field)

        make_l2_projection(target, source)

        name = op.join(self.options.out_dir,
                       'test_projection_tri_quad_target.vtk')
        target.save_as_mesh(name)

        bbox = self.field.domain.get_mesh_bounding_box()
        x = nm.linspace(bbox[0, 0] + 0.001, bbox[1, 0] - 0.001, 20)
        y = nm.linspace(bbox[0, 1] + 0.001, bbox[1, 1] - 0.001, 20)

        xx, yy = nm.meshgrid(x, y)
        test_coors = nm.c_[xx.ravel(), yy.ravel()].copy()

        vec1 = source.evaluate_at(test_coors)
        vec2 = target.evaluate_at(test_coors)

        ok = (nm.abs(vec1 - vec2) < 0.01).all()

        return ok

    def test_variables(self):
        from sfepy.discrete import FieldVariable, Integral

        ok = True

        u = FieldVariable('u', 'parameter', self.field,
                          primary_var_name='(set-to-None)')
        u.set_constant(1.0)
        vec = u() # Nodal values.

        _ok = nm.allclose(vec, 1.0)
        self.report('set constant:', _ok)
        ok = _ok and ok

        def fun(coors):
            val = nm.empty_like(coors)
            val[:, 0] = 2 * coors[:, 1] - coors[:, 0]
            val[:, 1] = coors[:, 0] + 3 * coors[:, 1]
            return val
        u.set_from_function(fun)

        coors = u.field.get_coor()
        eu = u.evaluate_at(coors)

        _ok = nm.allclose(eu, fun(coors), rtol=0.0, atol=1e-13)
        self.report('set from function:', _ok)
        ok = _ok and ok

        integral = Integral('i', order=2)
        gu_qp = u.evaluate(mode='grad', integral=integral)

        # du_i/dx_j, i = column, j = row.
        gu = nm.array([[-1.,  1.],
                       [ 2.,  3.]])
        _ok = nm.allclose(gu_qp, gu[None, None, ...], rtol=0.0, atol=1e-13)
        self.report('set from function - gradient:', _ok)
        ok = _ok and ok

        u_qp = gu_qp[..., :, :1]
        u.set_from_qp(u_qp, integral)
        vu = u()

        _ok = (nm.allclose(vu[::2], -1, rtol=0.0, atol=1e-13) and
               nm.allclose(vu[1::2], 2, rtol=0.0, atol=1e-13))
        self.report('set from qp:', _ok)
        ok = _ok and ok

        return ok

def prepare_variable(filename, n_components):
    from sfepy.discrete import FieldVariable
    from sfepy.discrete.fem import Mesh, FEDomain, Field

    mesh = Mesh.from_file(filename)

    bbox = mesh.get_bounding_box()
    dd = bbox[1,:] - bbox[0,:]
    data = (nm.sin(4.0 * nm.pi * mesh.coors[:,0:1] / dd[0])
            * nm.cos(4.0 * nm.pi * mesh.coors[:,1:2] / dd[1]))

    domain = FEDomain('domain', mesh)
    omega = domain.create_region('Omega', 'all')
    field = Field.from_args('field', nm.float64, n_components, omega,
                            approx_order=2)

    u = FieldVariable('u', 'parameter', field,
                      primary_var_name='(set-to-None)')
    u.set_from_mesh_vertices(nm.c_[tuple([data] * n_components)])

    return u

    def test_project_tensors(self):
        from sfepy.discrete import FieldVariable
        from sfepy.discrete.projections import project_by_component

        ok = True

        u = FieldVariable('u', 'parameter', self.field,
                          primary_var_name='(set-to-None)')
        u.set_constant(1.0)

        component = FieldVariable('component', 'parameter', self.field,
                                  primary_var_name='(set-to-None)')

        nls_options = {'eps_a' : 1e-16, 'i_max' : 1}

        u_qp = u.evaluate()
        u2 = FieldVariable('u2', 'parameter', self.field,
                           primary_var_name='(set-to-None)')
        project_by_component(u2, u_qp, component, self.field.approx_order,
                             nls_options=nls_options)

        _ok = self.compare_vectors(u(), u2())
        ok = ok and _ok

        gu_qp = u.evaluate(mode='grad')

        gfield = Field.from_args('gu', nm.float64, 2, self.field.region,
                                 approx_order=self.field.approx_order)
        gu = FieldVariable('gu', 'parameter', gfield,
                           primary_var_name='(set-to-None)')

        project_by_component(gu, gu_qp, component, gfield.approx_order,
                             nls_options=nls_options)

        _ok = self.compare_vectors(gu(), nm.zeros_like(gu()))
        ok = ok and _ok

        return ok

def save_basis_on_mesh(mesh, options, output_dir, lin,
                       permutations=None, suffix=''):
    if permutations is not None:
        mesh = mesh.copy()
        gel = GeometryElement(mesh.descs[0])
        perms = gel.get_conn_permutations()[permutations]
        conn = mesh.cmesh.get_cell_conn()
        n_el, n_ep = conn.num, gel.n_vertex
        offsets = nm.arange(n_el) * n_ep

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

    domain = FEDomain('domain', mesh)

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

    if options.plot_dofs:
        import sfepy.postprocess.plot_dofs as pd
        import sfepy.postprocess.plot_cmesh as pc
        ax = pc.plot_wireframe(None, mesh.cmesh)
        ax = pd.plot_global_dofs(ax, field.get_coor(), field.econn)
        ax = pd.plot_local_dofs(ax, field.get_coor(), field.econn)
        if options.dofs is not None:
            ax = pd.plot_nodes(ax, field.get_coor(), field.econn,
                               field.poly_space.nodes,
                               get_dofs(options.dofs, var.n_dof))
        pd.plt.show()

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

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

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

        var.set_data(vec)

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

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

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

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

def make_term_args(arg_shapes, arg_kinds, arg_types, ats_mode, domain,
                   material_value=None, poly_space_base=None):
    from sfepy.base.base import basestr
    from sfepy.discrete import FieldVariable, Material, Variables, Materials
    from sfepy.discrete.fem import Field
    from sfepy.solvers.ts import TimeStepper
    from sfepy.mechanics.tensors import dim2sym

    omega = domain.regions['Omega']
    dim = domain.shape.dim
    sym = dim2sym(dim)

    def _parse_scalar_shape(sh):
        if isinstance(sh, basestr):
            if sh == 'D':
                return dim

            elif sh == 'D2':
                return dim**2

            elif sh == 'S':
                return sym

            elif sh == 'N': # General number ;)
                return 1

            else:
                return int(sh)

        else:
            return sh

    def _parse_tuple_shape(sh):
        if isinstance(sh, basestr):
            return [_parse_scalar_shape(ii.strip()) for ii in sh.split(',')]

        else:
            return (int(sh),)

    args = {}
    str_args = []
    materials = []
    variables = []
    for ii, arg_kind in enumerate(arg_kinds):
        if arg_kind != 'ts':
            if ats_mode is not None:
                extended_ats = arg_types[ii] + ('/%s' % ats_mode)

            else:
                extended_ats = arg_types[ii]

            try:
                sh = arg_shapes[arg_types[ii]]

            except KeyError:
                sh = arg_shapes[extended_ats]

        if arg_kind.endswith('variable'):
            shape = _parse_scalar_shape(sh[0] if isinstance(sh, tuple) else sh)
            field = Field.from_args('f%d' % ii, nm.float64, shape, omega,
                                    approx_order=1,
                                    poly_space_base=poly_space_base)

            if arg_kind == 'virtual_variable':
                if sh[1] is not None:
                    istate = arg_types.index(sh[1])

                else:
                    # Only virtual variable in arguments.
                    istate = -1
                    # -> Make fake variable.
                    var = FieldVariable('u-1', 'unknown', field)
                    var.set_constant(0.0)
                    variables.append(var)

                var = FieldVariable('v', 'test', field,
                                    primary_var_name='u%d' % istate)

            elif arg_kind == 'state_variable':
                var = FieldVariable('u%d' % ii, 'unknown', field)
                var.set_constant(0.0)

            elif arg_kind == 'parameter_variable':
                var = FieldVariable('p%d' % ii, 'parameter', field,
                                    primary_var_name='(set-to-None)')
                var.set_constant(0.0)

            variables.append(var)
            str_args.append(var.name)
            args[var.name] = var

        elif arg_kind.endswith('material'):
            if sh is None: # Switched-off opt_material.
                continue

            prefix = ''
            if isinstance(sh, basestr):
                aux = sh.split(':')
                if len(aux) == 2:
                    prefix, sh = aux
#.........这里部分代码省略.........

#.........这里部分代码省略.........
                                              verbose=True)

    output('...done in', time.clock() - tt)

    output('creating solver...')
    tt = time.clock()

    conf = Struct(method='bcgsl', precond='jacobi', sub_precond=None,
                  i_max=10000, eps_a=1e-50, eps_r=1e-6, eps_d=1e4,
                  verbose=True)
    status = {}
    ls = PETScKrylovSolver(conf, comm=comm, mtx=pmtx, status=status)

    field_ranges = {}
    for ii, variable in enumerate(variables.iter_state(ordered=True)):
        field_ranges[variable.name] = lfds[ii].petsc_dofs_range

    ls.set_field_split(field_ranges, comm=comm)

    ev = PETScParallelEvaluator(pb, pdofs, drange, True,
                                psol, comm, verbose=True)

    nls_status = {}
    conf = Struct(method='newtonls',
                  i_max=5, eps_a=0, eps_r=1e-5, eps_s=0.0,
                  verbose=True)
    nls = PETScNonlinearSolver(conf, pmtx=pmtx, prhs=prhs, comm=comm,
                               fun=ev.eval_residual,
                               fun_grad=ev.eval_tangent_matrix,
                               lin_solver=ls, status=nls_status)

    output('...done in', time.clock() - tt)

    output('solving...')
    tt = time.clock()

    state = pb.create_state()
    state.apply_ebc()

    ev.psol_i[...] = state()
    ev.gather(psol, ev.psol_i)

    psol = nls(psol)

    ev.scatter(ev.psol_i, psol)
    sol0_i = ev.psol_i[...]

    output('...done in', time.clock() - tt)

    output('saving solution...')
    tt = time.clock()

    state.set_full(sol0_i)
    out = state.create_output_dict()

    filename = os.path.join(options.output_dir, 'sol_%02d.h5' % comm.rank)
    pb.domain.mesh.write(filename, io='auto', out=out)

    gather_to_zero = pl.create_gather_to_zero(psol)

    psol_full = gather_to_zero(psol)

    if comm.rank == 0:
        sol = psol_full[...].copy()

        u = FieldVariable('u', 'parameter', field1,
                          primary_var_name='(set-to-None)')
        remap = gfds[0].id_map
        ug = sol[remap]

        p = FieldVariable('p', 'parameter', field2,
                          primary_var_name='(set-to-None)')
        remap = gfds[1].id_map
        pg = sol[remap]

        if (((order_u == 1) and (order_p == 1))
            or (options.linearization == 'strip')):
            out = u.create_output(ug)
            out.update(p.create_output(pg))
            filename = os.path.join(options.output_dir, 'sol.h5')
            mesh.write(filename, io='auto', out=out)

        else:
            out = u.create_output(ug, linearization=Struct(kind='adaptive',
                                                           min_level=0,
                                                           max_level=order_u,
                                                           eps=1e-3))

            filename = os.path.join(options.output_dir, 'sol_u.h5')
            out['u'].mesh.write(filename, io='auto', out=out)

            out = p.create_output(pg, linearization=Struct(kind='adaptive',
                                                           min_level=0,
                                                           max_level=order_p,
                                                           eps=1e-3))

            filename = os.path.join(options.output_dir, 'sol_p.h5')
            out['p'].mesh.write(filename, io='auto', out=out)

    output('...done in', time.clock() - tt)

    def test_continuity(self):
        from sfepy.base.base import Struct
        from sfepy.mesh.mesh_generators import gen_block_mesh
        from sfepy.discrete import FieldVariable
        from sfepy.discrete.fem import FEDomain, Field
        from sfepy.discrete.projections import make_l2_projection_data
        import sfepy.discrete.fem.refine_hanging as rh

        dims = [1.5, 2.0, 1.3]
        shape = [3, 3, 3]
        centre = [0.0, 0.0, 0.0]
        probe_gens = {"2_4": _gen_lines_2_4, "3_8": _gen_grid_3_8}

        ok = True
        for key in self.gel_names:
            gel = self.gels[key]
            probe_gen = probe_gens[key]

            perms = gel.get_conn_permutations()

            dim = gel.dim
            for io, order in enumerate(range(1, 4)):
                mesh00 = gen_block_mesh(dims[:dim], shape[:dim], centre[:dim], name="block")

                for ip, perm in enumerate(perms):
                    self.report("geometry: %s, order: %d, permutation: %d: %s" % (key, order, ip, perm))
                    mesh0 = mesh00.copy()
                    conn = mesh0.cmesh.get_conn(dim, 0).indices
                    conn = conn.reshape((mesh0.n_el, -1))
                    conn[-1, :] = conn[-1, perm]

                    domain0 = FEDomain("d", mesh0)

                    refine = nm.zeros(mesh0.n_el, dtype=nm.uint8)
                    refine[:-1] = 1

                    subs = None
                    domain, subs = rh.refine(domain0, refine, subs=subs)

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

                    field.substitute_dofs(subs)

                    uvar = FieldVariable("u", "parameter", field, primary_var_name="(set-to-None)")

                    field.restore_dofs(store=True)
                    field.substitute_dofs(subs=None, restore=True)

                    make_l2_projection_data(uvar, eval_fun)

                    field.restore_dofs()

                    bbox = domain.get_mesh_bounding_box()
                    eps = 1e-7

                    save = False
                    for ii, (probe0, probe1) in enumerate(probe_gen(bbox, eps)):
                        probe0.set_options(close_limit=0.0)
                        probe1.set_options(close_limit=0.0)

                        pars0, vals0 = probe0(uvar)
                        pars1, vals1 = probe1(uvar)

                        assert_(nm.allclose(pars0, pars1, atol=1e-14, rtol=0.0))

                        _ok = nm.allclose(vals0, vals1, atol=20.0 * eps, rtol=0.0)
                        if not _ok:
                            save = True
                            self.report("probe %d failed! (max. error: %e)" % (ii, nm.abs(vals0 - vals1).max()))

                        ok = ok and _ok

                    if (ip == 0) or save:
                        out = uvar.create_output()
                        filenames = _build_filenames(self.options.out_dir, key, order, ip)

                        domain.mesh.write(filenames[0], out=out)

                        linearization = Struct(kind="adaptive", min_level=0, max_level=4, eps=1e-2)

                        out = uvar.create_output(linearization=linearization)
                        val = out["u"]
                        mesh = val.get("mesh", domain.mesh)
                        mesh.write(filenames[1], out=out)

        return ok