def test_write_read_meshes(self):
        """
        Try to write and then read all supported formats.
        """
        from sfepy.fem import Mesh
        from sfepy.fem.meshio import supported_formats, supported_capabilities

        conf_dir = op.dirname(__file__)
        mesh0 = Mesh.from_file(data_dir
                               + '/meshes/various_formats/small3d.mesh',
                               prefix_dir=conf_dir)

        oks = []
        for suffix, format_ in supported_formats.iteritems():
            if isinstance(format_, tuple):
                continue
            if 'w' not in supported_capabilities[format_]: continue

            filename = op.join(self.options.out_dir, 'test_mesh_wr' + suffix)
            self.report('%s format: %s' % (suffix, filename))

            mesh0.write(filename, io='auto')
            mesh1 = Mesh.from_file(filename)

            oks.extend(self._compare_meshes(mesh0, mesh1))

        return sum(oks) == len(oks)

    def test_interpolation_two_meshes(self):
        from sfepy import data_dir
        from sfepy.fem import Mesh, Domain, H1NodalVolumeField, Variables

        m1 = Mesh('source mesh', data_dir + '/meshes/3d/block.mesh')

        m2 = Mesh('target mesh', data_dir + '/meshes/3d/cube_medium_tetra.mesh')
        m2.coors *= 2.0

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

        variables1 = {
            'u'       : ('unknown field', 'scalar_tp', 0),
            'v'       : ('test field',    'scalar_tp', 'u'),
        }

        variables2 = {
            'u'       : ('unknown field', 'scalar_si', 0),
            'v'       : ('test field',    'scalar_si', 'u'),
        }

        d1 = Domain('d1', m1)
        omega1 = d1.create_region('Omega', 'all')
        field1 = H1NodalVolumeField('scalar_tp', nm.float64, (1,1), omega1,
                                    approx_order=1)
        ff1 = {field1.name : field1}

        d2 = Domain('d2', m2)
        omega2 = d2.create_region('Omega', 'all')
        field2 = H1NodalVolumeField('scalar_si', nm.float64, (1,1), omega2,
                                    approx_order=0)
        ff2 = {field2.name : field2}

        vv1 = Variables.from_conf(transform_variables(variables1), ff1)
        u1 = vv1['u']
        u1.set_from_mesh_vertices(data)

        vv2 = Variables.from_conf(transform_variables(variables2), ff2)
        u2 = vv2['u']

        # Performs interpolation, if other field differs from self.field
        # or, in particular, is defined on a different mesh.
        u2.set_from_other(u1, strategy='interpolation', close_limit=0.1)

        fname = in_dir(self.options.out_dir)
        u1.save_as_mesh(fname('test_mesh_interp_block_scalar.vtk'))
        u2.save_as_mesh(fname('test_mesh_interp_cube_scalar.vtk'))

        return True

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 from_conf(conf, options):
        import sfepy
        from sfepy.fem import Mesh, Domain, H1NodalVolumeField
        mesh = Mesh.from_file('meshes/2d/rectangle_tri.mesh',
                              prefix_dir=sfepy.data_dir)
        domain = Domain('domain', mesh)
        dim = domain.shape.dim

        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 = H1NodalVolumeField('fu', nm.float64, 'vector', omega,
                                   approx_order=2)

        test = Test(conf=conf, options=options, dim=dim,
                    omega=omega, gamma1=gamma1, gamma2=gamma2,
                    field=field)
        return test

    def save_axes(self, filename):
        coors = []
        conns = []
        mat_ids = []
        offset = 0
        for ig, dual_surface in self.dual_surfaces.iteritems():
            cc = nm.r_[dual_surface.edge_centre_coors,
                       dual_surface.dual_coors]
            coors.append(cc)

            conn = dual_surface.conn.copy() + offset
            conn[:,1:] += dual_surface.edge_centre_coors.shape[0]
            conns.append(conn)
            
            mat_id = nm.empty((conn.shape[0],), dtype=nm.int32)
            mat_id[:] = ig
            mat_ids.append(mat_id)

            offset += cc.shape[0]

        coors = nm.concatenate(coors, axis=0)

        out = {}
        for ig, dual_surface in self.dual_surfaces.iteritems():
            eto = edge_data_to_output
            out['en_%d' % ig] = eto(coors, conns[ig], dual_surface.e_sort,
                                    dual_surface.edge_normals)
            out['ed_%d' % ig] = eto(coors, conns[ig], dual_surface.e_sort,
                                    dual_surface.edge_dirs)
            out['eo_%d' % ig] = eto(coors, conns[ig], dual_surface.e_sort,
                                    dual_surface.edge_ortho)

        dual_mesh = Mesh.from_data('dual_mesh_vectors', coors, None, conns,
                                   mat_ids, ['2_4'] * len(conns))
        dual_mesh.write(filename, io='auto', out=out)

    def test_read_meshes( self ):
        """Try to read all listed meshes."""
        from sfepy.fem import Mesh

        conf_dir = op.dirname(__file__)
        meshes = {}
        for ii, filename in enumerate( filename_meshes ):
            self.report( '%d. mesh: %s' % (ii + 1, filename) )
            mesh = Mesh.from_file(filename, prefix_dir=conf_dir)

            assert_(mesh.dim == (mesh.coors.shape[1]))
            assert_(mesh.n_nod == (mesh.coors.shape[0]))
            assert_(mesh.n_nod == (mesh.ngroups.shape[0]))
            assert_(mesh.n_el == sum(mesh.n_els))
            for ig, conn in enumerate( mesh.conns ):
                assert_(conn.shape[0] == len(mesh.mat_ids[ig]))
                assert_(conn.shape[0] == mesh.n_els[ig])
                assert_(conn.shape[1] == mesh.n_e_ps[ig])
                
            self.report( 'read ok' )
            meshes[filename] = mesh

        self.meshes = meshes

        return True

    def test_rcm(self):
        from sfepy import data_dir
        from sfepy.linalg import rcm, permute_in_place, save_sparse_txt
        from sfepy.fem import Mesh

        filename = data_dir + '/meshes/2d/special/square_triquad.mesh'

        self.report('testing reversed Cuthill-McKee permutation')

        conf_dir = op.dirname(__file__)
        mesh = Mesh.from_file(filename, prefix_dir=conf_dir)

        graph = mesh.create_conn_graph()
        graph0 = graph.copy()

        save_sparse_txt(op.join(self.options.out_dir, 'test_rcm_graph_orig'),
                        graph, fmt='%d %d %d\n')

        perm = rcm(graph)

        permute_in_place(graph, perm)
        save_sparse_txt(op.join(self.options.out_dir, 'test_rcm_graph_rcm'),
                        graph, fmt='%d %d %d\n')
        
        assert_((graph0.indptr != graph.indptr).any())
        assert_((graph0.indices != graph.indices).any())

        permute_in_place(graph, perm, inverse=True)
        save_sparse_txt(op.join(self.options.out_dir, 'test_rcm_graph_rcm_inv'),
                        graph, fmt='%d %d %d\n')

        assert_((graph0.indptr == graph.indptr).all())
        assert_((graph0.indices == graph.indices).all())

        return True

def refine_3_8(mesh_in, ed, fa):
    """
    Refines hexahedral mesh by cutting cutting each edge in half and
    making 8 new finer hexahedrons out of one coarser one.
    """
    # Unique edge centres.
    e_coors, e_uid = ed.get_coors()
    e_centres = 0.5 * nm.sum(e_coors, axis=1)

    # Unique face centres.
    f_coors, f_uid = fa.get_coors()
    f_centres = 0.25 * nm.sum(f_coors, axis=1)

    # Unique element centres.
    coors = mesh_in.get_element_coors()
    centres = 0.125 * nm.sum(coors, axis=1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres, f_centres, centres]

    o1 = mesh_in.n_nod
    o2 = o1 + e_centres.shape[0]
    o3 = o2 + f_centres.shape[0]

    st = nm.vstack

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        e_indx = ed.indx[ig]
        f_indx = fa.indx[ig]
        off = mesh_in.el_offsets[ig]
        n_el  = conn.shape[0]

        e_nodes = ed.uid_i[e_indx].reshape((n_el, 12)) + o1
        f_nodes = fa.uid_i[f_indx].reshape((n_el, 6)) + o2
        nodes = nm.arange(n_el) + off + o3

        c = nm.c_[conn, e_nodes, f_nodes, nodes].T

        new_conn = st([c[0], c[8], c[20], c[11], c[16], c[22], c[26], c[21],
                       c[1], c[9], c[20], c[8], c[17], c[24], c[26], c[22],
                       c[2], c[10], c[20], c[9], c[18], c[25], c[26], c[24],
                       c[3], c[11], c[20], c[10], c[19], c[21], c[26], c[25],
                       c[4], c[15], c[23], c[12], c[16], c[21], c[26], c[22],
                       c[5], c[12], c[23], c[13], c[17], c[22], c[26], c[24],
                       c[6], c[13], c[23], c[14], c[18], c[24], c[26], c[25],
                       c[7], c[14], c[23], c[15], c[19], c[25], c[26], c[21]]).T
        new_conn = new_conn.reshape((8 * n_el, 8))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(8)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, conns,
                          mat_ids, mesh_in.descs )

    return mesh

    def test_invariance_qp(self):
        from sfepy import data_dir
        from sfepy.fem import (Mesh, Domain, H1NodalVolumeField,
                               Variables, Integral)
        from sfepy.terms import Term
        from sfepy.fem.mappings import get_physical_qps

        mesh = Mesh('source mesh', data_dir + '/meshes/3d/block.mesh')

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

        variables = {
            'u'       : ('unknown field', 'scalar_tp', 0),
            'v'       : ('test field',    'scalar_tp', 'u'),
        }

        domain = Domain('domain', mesh)
        omega = domain.create_region('Omega', 'all')
        field = H1NodalVolumeField('scalar_tp', nm.float64, 1, omega,
                                   approx_order=1)
        ff = {field.name : field}

        vv = Variables.from_conf(transform_variables(variables), ff)
        u = vv['u']
        u.set_from_mesh_vertices(data)

        integral = Integral('i', order=2)
        term = Term.new('ev_volume_integrate(u)', integral, omega, u=u)
        term.setup()
        val1, _ = term.evaluate(mode='qp')
        val1 = val1.ravel()

        qps = get_physical_qps(omega, integral)
        coors = qps.get_merged_values()

        val2 = u.evaluate_at(coors).ravel()

        self.report('max. difference:', nm.abs(val1 - val2).max())
        ok = nm.allclose(val1, val2, rtol=0.0, atol=1e-12)
        self.report('invariance in qp: %s' % ok)

        return ok

    def from_conf( conf, options ):
        from sfepy import data_dir
        from sfepy.fem import Mesh, Domain, Functions

        mesh = Mesh('test mesh',
                    data_dir + '/meshes/various_formats/abaqus_tet.inp')
        mesh.nodal_bcs['set0'] = [0, 7]
        domain = Domain('test domain', mesh)

        conf_functions = {
            'get_vertices' : (get_vertices,),
            'get_cells' : (get_cells,),
        }
        functions = Functions.from_conf(transform_functions(conf_functions))

        test = Test(conf=conf, options=options,
                    domain=domain, functions=functions)
        return test

    def from_conf(conf, options):
        mesh = Mesh.from_file("meshes/2d/square_unit_tri.mesh", prefix_dir=sfepy.data_dir)
        domain = Domain("domain", mesh)

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

        field = H1NodalVolumeField("linear", nm.float64, "scalar", omega, approx_order=1)

        test = Test(conf=conf, options=options, omega=omega, field=field)
        return test

def refine_3_4(mesh_in, ed):
    """
    Refines tetrahedra by cutting each edge in half and making 8 new
    finer tetrahedra out of one coarser one. Old nodal coordinates come
    first in `coors`, then the new ones. The new tetrahedra are similar
    to the old one, no degeneration is supposed to occur as at most 3
    congruence classes of tetrahedra appear, even when re-applied
    iteratively (provided that `conns` are not modified between two
    applications - ordering of vertices in tetrahedra matters not only
    for positivity of volumes).

    References:

    - Juergen Bey: Simplicial grid refinement: on Freudenthal s algorithm and 
      the optimal number of congruence classes, Numer.Math. 85 (2000), 
      no. 1, 1--29, or
    - Juergen Bey: Tetrahedral grid refinement, Computing 55 (1995), 
      no. 4, 355--378, or
      http://citeseer.ist.psu.edu/bey95tetrahedral.html
    """
    # Unique edge centres.
    e_coors, e_uid = ed.get_coors()
    e_centres = 0.5 * nm.sum(e_coors, axis=1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres]

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        indx = ed.indx[ig]
        n_el  = conn.shape[0]

        e_nodes = ed.uid_i[indx].reshape((n_el, 6)) + mesh_in.n_nod

        c = nm.c_[conn, e_nodes].T

        new_conn = nm.vstack([c[0], c[4], c[6], c[7],
                              c[4], c[1], c[5], c[8],
                              c[6], c[5], c[2], c[9],
                              c[7], c[8], c[9], c[3],
                              c[4], c[6], c[7], c[8],
                              c[4], c[6], c[8], c[5],
                              c[6], c[7], c[8], c[9],
                              c[6], c[5], c[9], c[8]]).T
        new_conn = new_conn.reshape((8 * n_el, 4))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(8)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, conns,
                          mat_ids, mesh_in.descs )

    return mesh

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", "nodes in x < %.10f" % (min_x + eps))
    gamma2 = domain.create_region("Gamma2", "nodes in x > %.10f" % (max_x - eps))

    field = Field("fu", nm.float64, "vector", omega, space="H1", poly_space_base="lagrange", approx_order=2)

    u = FieldVariable("u", "unknown", field, mesh.dim)
    v = FieldVariable("v", "test", field, mesh.dim, 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 = ProblemDefinition("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 from_conf(conf, options):
        mesh = Mesh.from_file('meshes/2d/square_unit_tri.mesh',
                              prefix_dir=sfepy.data_dir)
        domain = Domain('domain', mesh)

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

        field = Field('linear', nm.float64, 'scalar', omega,
                      space='H1', poly_space_base='lagrange', approx_order=1)

        test = Test(conf=conf, options=options, omega=omega, field=field)
        return test

def main():
    parser = OptionParser(usage=usage)
    parser.add_option("-s", "--scale", metavar='scale',
                      action="store", dest="scale",
                      default=None, help=help['scale'])
    parser.add_option("-f", "--format", metavar='format',
                      action="store", type='string', dest="format",
                      default=None, help=help['format'])
    parser.add_option("-l", "--list", action="store_true", 
                      dest="list", help=help['list'])
    (options, args) = parser.parse_args()

    if options.list:
        output_writable_meshes()
        sys.exit(0)

    if len(args) != 2:
        parser.print_help()
        sys.exit(1)
    
    scale = options.scale
    if scale is not None:
        try:
            try:
                scale = float(scale)
            except ValueError:
                scale = [float(ii) for ii in scale.split(',')]
            scale = nm.array(scale, dtype=nm.float64, ndmin=1)
        except:
            output('bad scale! (%s)' % scale)
            parser.print_help()
            sys.exit(1)
        
    filename_in, filename_out = args
    
    mesh = Mesh.from_file(filename_in)

    if scale is not None:
        if len(scale) == 1:
            tr = nm.eye(mesh.dim, dtype=nm.float64) * scale
        elif len(scale) == mesh.dim:
            tr = nm.diag(scale)
        else:
            raise ValueError('bad scale! (%s)' % scale)
        mesh.transform_coors(tr)

    io = MeshIO.for_format(filename_out, format=options.format,
                           writable=True)

    output('writing %s...' % filename_out)
    mesh.write(filename_out, io=io)
    output('...done')

    def test_read_meshes( self ):
        """Try to read all listed meshes."""
        from sfepy.fem import Mesh

        meshes = {}
        for ii, filename in enumerate( filename_meshes ):
            self.report( '%d. mesh: %s' % (ii + 1, filename) )
            mesh = Mesh.from_file( filename )
            self.report( 'read ok' )
            meshes[filename] = mesh

        self.meshes = meshes

        return True

    def test_normals(self):
        """
        Check orientations of surface normals on the reference elements.
        """
        import sfepy
        from sfepy.fem import Mesh, Domain, Integral
        from sfepy.fem.poly_spaces import PolySpace
        from sfepy.fem.mappings import SurfaceMapping
        from sfepy.linalg import normalize_vectors

        ok = True

        for geom in ['2_3', '2_4', '3_4', '3_8']:
            mesh = Mesh.from_file('meshes/elements/%s_1.mesh' % geom,
                                  prefix_dir=sfepy.data_dir)
            domain = Domain('domain', mesh)
            surface = domain.create_region('Surface', 'nodes of surface')
            domain.create_surface_group(surface)

            sd = domain.surface_groups[0][surface.name]

            coors = domain.get_mesh_coors()
            gel = domain.geom_els[geom].surface_facet
            ps = PolySpace.any_from_args('aux', gel, 1)

            mapping = SurfaceMapping(coors, sd.get_connectivity(), ps)

            integral = Integral('i', order=1)
            vals, weights = integral.get_qp(gel.name)

            # Evaluate just in the first quadrature point...
            geo = mapping.get_mapping(vals[:1], weights[:1])

            expected = expected_normals[geom].copy()
            normalize_vectors(expected)

            _ok = nm.allclose(expected, geo.normal[:, 0, :, 0],
                              rtol=0.0, atol=1e-14)
            self.report('%s: %s' % (geom, _ok))

            if not _ok:
                self.report('expected:')
                self.report(expected)
                self.report('actual:')
                self.report(geo.normal[:, 0, :, 0])

            ok = ok and _ok

        return ok

def mesh():
    if len(sys.argv) == 3:
        geom_filename = sys.argv[1]
        vtk_filename = sys.argv[2]
    else:
        print "Usage: %s <gmsh_filename> <mesh_filename>" % sys.argv[0]
        return

    os.system("gmsh -0 %s -o tmp/x.geo" % geom_filename)
    g = geom.read_gmsh("tmp/x.geo")
    g.printinfo()
    geom.write_tetgen(g, "tmp/t.poly")
    geom.runtetgen("tmp/t.poly", a=0.03, Q=1.0, quadratic=False, tetgenpath=tetgen_path)

    m = Mesh.from_file("tmp/t.1.node")
    m.write(vtk_filename, io="auto")

def refine_2_4(mesh_in, ed):
    """
    Refines mesh out of quadrilaterals by cutting cutting each edge in
    half and making 4 new finer quadrilaterals out of one coarser one.
    """
    # Unique edge centres.
    e_coors, e_uid = ed.get_coors()
    e_centres = 0.5 * nm.sum(e_coors, axis=1)

    # Unique element centres.
    coors = mesh_in.get_element_coors()
    centres = 0.25 * nm.sum(coors, axis=1)

    # New coordinates after the original ones.
    coors = nm.r_[mesh_in.coors, e_centres, centres]

    o1 = mesh_in.n_nod
    o2 = o1 + e_centres.shape[0]

    conns = []
    mat_ids = []
    for ig, conn in enumerate(mesh_in.conns):
        e_indx = ed.indx[ig]
        off = mesh_in.el_offsets[ig]
        n_el  = conn.shape[0]

        e_nodes = ed.uid_i[e_indx].reshape((n_el, 4)) + o1
        nodes = nm.arange(n_el) + off + o2

        c = nm.c_[conn, e_nodes, nodes].T

        new_conn = nm.vstack([c[0], c[4], c[8], c[7],
                              c[1], c[5], c[8], c[4],
                              c[2], c[6], c[8], c[5],
                              c[3], c[7], c[8], c[6]]).T
        new_conn = new_conn.reshape((4 * n_el, 4))
        conns.append(new_conn)

        new_mat_id = mesh_in.mat_ids[ig].repeat(4)
        mat_ids.append(new_mat_id)

    mesh = Mesh.from_data(mesh_in.name + '_r', coors, None, conns,
                          mat_ids, mesh_in.descs )

    return mesh

    def from_conf(conf, options):
        from sfepy.fem import Mesh, Domain, Integral

        domains = []
        for filename in filename_meshes:
            mesh = Mesh.from_file(filename)
            domain = Domain('domain_%s' % mesh.name.replace(data_dir, ''),
                            mesh)
            domain.create_region('Omega', 'all')
            domain.create_region('Gamma', 'vertices of surface', 'facet')

            domains.append(domain)

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

        test = Test(domains=domains, integral=integral,
                    conf=conf, options=options)
        return test