def test_interpolation(self):
        from sfepy import data_dir
        from sfepy.discrete.fem import Mesh
        from sfepy.linalg import make_axis_rotation_matrix

        fname = in_dir(self.options.out_dir)

        meshes = {
            'tp' : Mesh.from_file(data_dir + '/meshes/3d/block.mesh'),
            'si' : Mesh.from_file(data_dir + '/meshes/3d/cylinder.mesh'),
        }

        datas = gen_datas(meshes)

        for field_name in ['scalar_si', 'vector_si', 'scalar_tp', 'vector_tp']:
            m1 = meshes[field_name[-2:]]
            for ia, angle in enumerate(nm.linspace(0.0, nm.pi, 11)):
                self.report('%s: %d. angle: %f' % (field_name, ia, angle))
                shift = [0.0, 0.0, 0.0]
                mtx = make_axis_rotation_matrix([0, 1, 0], angle)

                m2 = m1.copy('rotated mesh')
                m2.transform_coors(mtx)

                data = datas[field_name]
                u1, u2 = do_interpolation(m2, m1, data, field_name)

                if ia == 0:
                    u1.save_as_mesh(fname('test_mesh_interp_%s_u1.vtk'
                                          % field_name))

                u2.save_as_mesh(fname('test_mesh_interp_%s_u2.%03d.vtk'
                                      % (field_name, ia)))

        return True

    def test_write_read_meshes(self):
        """
        Try to write and then read all supported formats.
        """
        from sfepy.discrete.fem import Mesh
        from sfepy.discrete.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 six.iteritems(supported_formats):
            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 main():
    parser = ArgumentParser(description=__doc__.rstrip(),
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('filename', help=helps['filename'])
    parser.add_argument('-d', '--detailed',
                        action='store_true', dest='detailed',
                        default=False, help=helps['detailed'])
    options = parser.parse_args()

    mesh = Mesh.from_file(options.filename)

    output(mesh.cmesh)
    output('element types:', mesh.descs)
    output('nodal BCs:', sorted(mesh.nodal_bcs.keys()))

    bbox = mesh.get_bounding_box()
    output('bounding box:\n%s'
           % '\n'.join('%s: [%14.7e, %14.7e]' % (name, bbox[0, ii], bbox[1, ii])
                       for ii, name in enumerate('xyz'[:mesh.dim])))

    output('centre:           [%s]'
           % ', '.join('%14.7e' % ii for ii in 0.5 * (bbox[0] + bbox[1])))
    output('coordinates mean: [%s]'
           % ', '.join('%14.7e' % ii for ii in mesh.coors.mean(0)))

    if not options.detailed: return

    domain = FEDomain(mesh.name, mesh)

    for dim in range(1, mesh.cmesh.tdim + 1):
        volumes = mesh.cmesh.get_volumes(dim)
        output('volumes of %d %dD entities:\nmin: %.7e mean: %.7e median:'
               ' %.7e max: %.7e'
               % (mesh.cmesh.num[dim], dim, volumes.min(), volumes.mean(),
                  nm.median(volumes), volumes.max()))

    euler = lambda mesh: nm.dot(mesh.cmesh.num, [1, -1, 1, -1])
    ec = euler(mesh)
    output('Euler characteristic:', ec)

    graph = mesh.create_conn_graph(verbose=False)
    n_comp, _ = graph_components(graph.shape[0], graph.indptr, graph.indices)
    output('number of connected components:', n_comp)

    if mesh.dim > 1:
        region = domain.create_region('surf', 'vertices of surface', 'facet')
        surf_mesh = Mesh.from_region(region, mesh,
                                     localize=True, is_surface=True)
        FEDomain(surf_mesh.name, surf_mesh) # Calls CMesh.setup_entities().

        sec = euler(surf_mesh)
        output('surface Euler characteristic:', sec)
        if mesh.dim == 3:
            output('surface genus:', (2.0 - sec) / 2.0)

        surf_graph = surf_mesh.create_conn_graph(verbose=False)
        n_comp, _ = graph_components(surf_graph.shape[0],
                                     surf_graph.indptr, surf_graph.indices)
        output('number of connected surface components:', n_comp)

def refine_region(domain0, region0, region1):
    """
    Coarse cell sub_cells[ii, 0] in mesh0 is split into sub_cells[ii, 1:] in
    mesh1.

    The new fine cells are interleaved among the original coarse cells so that
    the indices of the coarse cells do not change.

    The cell groups are preserved. The vertex groups are preserved only in the
    coarse (non-refined) cells.
    """
    if region1 is None:
        return domain0, None

    mesh0 = domain0.mesh
    mesh1 = Mesh.from_region(region1, mesh0)
    domain1 = FEDomain('d', mesh1)
    domain1r = domain1.refine()
    mesh1r = domain1r.mesh

    n_cell = region1.shape.n_cell
    n_sub = 4 if mesh0.cmesh.tdim == 2 else 8

    sub_cells = nm.empty((n_cell, n_sub + 1), dtype=nm.uint32)
    sub_cells[:, 0] = region1.cells
    sub_cells[:, 1] = region1.cells
    aux = nm.arange((n_sub - 1) * n_cell, dtype=nm.uint32)
    sub_cells[:, 2:] = mesh0.n_el + aux.reshape((n_cell, -1))

    coors0, vgs0, conns0, mat_ids0, descs0 = mesh0._get_io_data()
    coors, vgs, _conns, _mat_ids, descs = mesh1r._get_io_data()

    # Preserve vertex groups of non-refined cells.
    vgs[:len(vgs0)] = vgs0

    def _interleave_refined(c0, c1):
        if c1.ndim == 1:
            c0 = c0[:, None]
            c1 = c1[:, None]

        n_row, n_col = c1.shape
        n_new = region0.shape.n_cell + n_row

        out = nm.empty((n_new, n_col), dtype=c0.dtype)
        out[region0.cells] = c0[region0.cells]
        out[region1.cells] = c1[::n_sub]
        aux = c1.reshape((-1, n_col * n_sub))
        out[mesh0.n_el:] = aux[:, n_col:].reshape((-1, n_col))

        return out

    conn = _interleave_refined(conns0[0], _conns[0])
    mat_id = _interleave_refined(mat_ids0[0], _mat_ids[0]).squeeze()

    mesh = Mesh.from_data('a', coors, vgs, [conn], [mat_id], descs)
    domain = FEDomain('d', mesh)

    return domain, sub_cells

    def test_interpolation_two_meshes(self):
        from sfepy import data_dir
        from sfepy.discrete import Variables
        from sfepy.discrete.fem import Mesh, FEDomain, Field

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

        m2 = Mesh.from_file(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 = FEDomain('d1', m1)
        omega1 = d1.create_region('Omega', 'all')
        field1 = Field.from_args('scalar_tp', nm.float64, (1,1), omega1,
                                 approx_order=1)
        ff1 = {field1.name : field1}

        d2 = FEDomain('d2', m2)
        omega2 = d2.create_region('Omega', 'all')
        field2 = Field.from_args('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 __init__(self, name, nurbs, bmesh, regions=None, **kwargs):
        """
        Create an IGA domain.

        Parameters
        ----------
        name : str
            The domain name.
        """
        Domain.__init__(self, name, nurbs=nurbs, bmesh=bmesh, regions=regions,
                        **kwargs)
        from sfepy.discrete.fem.geometry_element import create_geometry_elements
        from sfepy.discrete.fem import Mesh
        from sfepy.discrete.fem.utils import prepare_remap

        tconn = iga.get_bezier_topology(bmesh.conn, nurbs.degrees)
        itc = nm.unique(tconn)

        remap = prepare_remap(itc, bmesh.conn.max() + 1)

        ltcoors = bmesh.cps[itc]
        ltconn = remap[tconn]

        n_nod, dim = ltcoors.shape
        n_el = ltconn.shape[0]
        self.shape = Struct(n_nod=n_nod, dim=dim, tdim=0, n_el=n_el)

        desc = '%d_%d' % (dim, bmesh.conn.shape[1])
        mat_id = nm.zeros(bmesh.conn.shape[0], dtype=nm.int32)
        eval_mesh = Mesh.from_data(self.name + '_eval', nurbs.cps, None,
                                   [nurbs.conn], [mat_id], [desc])
        self.eval_mesh = eval_mesh

        desc = '%d_%d' % (dim, 2**dim)
        mat_id = nm.zeros(ltconn.shape[0], dtype=nm.int32)
        self.mesh = Mesh.from_data(self.name + '_topo', ltcoors, None, [ltconn],
                                   [mat_id], [desc])

        self.cmesh = self.mesh.cmesh
        gels = create_geometry_elements()
        self.cmesh.set_local_entities(gels)
        self.cmesh.setup_entities()

        self.shape.tdim = self.cmesh.tdim

        self.gel = gels[desc]

        if regions is not None:
            self.vertex_set_bcs = {}
            for key, val in self.regions.iteritems():
                self.vertex_set_bcs[key] = remap[val]

        self.reset_regions()

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

        mesh = Mesh("test mesh", data_dir + "/meshes/various_formats/abaqus_tet.inp")
        mesh.nodal_bcs["set0"] = [0, 7]
        domain = FEDomain("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 test_rcm(self):
        from sfepy import data_dir
        from sfepy.linalg import rcm, permute_in_place, save_sparse_txt
        from sfepy.discrete.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 _get_bqp(geometry, order):
    from sfepy.discrete import Integral
    from sfepy.discrete.fem.geometry_element import GeometryElement
    from sfepy.discrete.fem import Mesh, FEDomain, Field

    gel = GeometryElement(geometry)

    mesh = Mesh.from_data('aux', gel.coors, None,
                          [gel.conn[None, :]], [[0]], [geometry])
    domain = FEDomain('domain', mesh)
    omega = domain.create_region('Omega', 'all')
    surf =  domain.create_region('Surf', 'vertices of surface', 'facet')
    field = Field.from_args('f', nm.float64, shape=1,
                            region=omega, approx_order=1)
    field.setup_surface_data(surf)

    integral = Integral('aux', order=order)
    field.create_bqp('Surf', integral)

    sd = field.surface_data['Surf']
    qp = field.qp_coors[(integral.order, sd.bkey)]

    output('geometry:', geometry, 'order:', order, 'num. points:',
           qp.vals.shape[1], 'true_order:',
           integral.qps[gel.surface_facet_name].order)
    output('min. weight:', qp.weights.min())
    output('max. weight:', qp.weights.max())

    return (gel, qp.vals.reshape((-1, mesh.dim)),
            nm.tile(qp.weights, qp.vals.shape[0]))

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 main():
    parser = ArgumentParser(description=__doc__,
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('mesh_dir')
    options = parser.parse_args()

    mesh_dir = options.mesh_dir

    mesh_files = []
    for (dirpath, dirnames, filenames) in os.walk(mesh_dir):
        for ii in filenames:
            _, ext = os.path.splitext(ii)
            if ext.lower() in ['.mesh', '.vtk']:
                mesh_files.append(dirpath + os.path.sep + ii)

    for ii in mesh_files:
        base, ext = os.path.splitext(ii)
        fname_out = base + '.png'
        if ext == '.mesh':
            fname_in = 'aux.vtk'
            mesh = Mesh.from_file(ii)
            mesh.write(fname_in, io='auto')

        else:
            fname_in = ii

        print(('writing %s...' % fname_out))
        gen_shot(fname_in, fname_out)

    def test_read_meshes(self):
        """Try to read all listed meshes."""
        from sfepy.discrete.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 main():
    parser = OptionParser(usage=usage)
    (options, args) = parser.parse_args()

    if len(args) != 1:
        parser.print_help()
        sys.exit(1)

    mesh_dir = args[0]

    mesh_files = []
    for (dirpath, dirnames, filenames) in os.walk(mesh_dir):
        for ii in filenames:
            _, ext = os.path.splitext(ii)
            if ext.lower() in ['.mesh', '.vtk']:
                mesh_files.append(dirpath + os.path.sep + ii)

    for ii in mesh_files:
        base, ext = os.path.splitext(ii)
        fname_out = base + '.png'
        if ext == '.mesh':
            fname_in = 'aux.vtk'
            mesh = Mesh.from_file(ii)
            mesh.write(fname_in, io='auto')

        else:
            fname_in = ii

        print('writing %s...' % fname_out)
        gen_shot(fname_in, fname_out)

def main():
    parser = ArgumentParser(description=__doc__)
    parser.add_argument('--version', action='version', version='%(prog)s')
    parser.add_argument('filename')
    options = parser.parse_args()

    filename = options.filename

    mesh = Mesh.from_file(filename)
    output('Mesh:')
    output('  dimension: %d, vertices: %d, elements: %d'
           % (mesh.dim, mesh.n_nod, mesh.n_el))

    domain = FEDomain('domain', mesh)
    output(domain.cmesh)
    domain.cmesh.cprint(1)
    dim = domain.cmesh.dim

    entities_opts = [
        {'color' : 'k', 'label_global' : 12, 'label_local' : 8},
        {'color' : 'b', 'label_global' : 12, 'label_local' : 8},
        {'color' : 'g', 'label_global' : 12, 'label_local' : 8},
        {'color' : 'r', 'label_global' : 12},
    ]
    if dim == 2: entities_opts.pop(2)

    pc.plot_cmesh(None, domain.cmesh,
                  wireframe_opts = {'color' : 'k'},
                  entities_opts=entities_opts)

    plt.show()

def gen_two_bodies(dims0, shape0, centre0, dims1, shape1, centre1, shift1):
    from sfepy.discrete.fem import Mesh
    from sfepy.mesh.mesh_generators import gen_block_mesh

    m0 = gen_block_mesh(dims0, shape0, centre0)
    m1 = gen_block_mesh(dims1, shape1, centre1)

    coors = nm.concatenate((m0.coors, m1.coors + shift1), axis=0)

    desc = m0.descs[0]
    c0 = m0.get_conn(desc)
    c1 = m1.get_conn(desc)
    conn = nm.concatenate((c0, c1 + m0.n_nod), axis=0)

    ngroups = nm.zeros(coors.shape[0], dtype=nm.int32)
    ngroups[m0.n_nod:] = 1

    mat_id = nm.zeros(conn.shape[0], dtype=nm.int32)
    mat_id[m0.n_el:] = 1

    name = 'two_bodies.mesh'

    mesh = Mesh.from_data(name, coors, ngroups, [conn], [mat_id], m0.descs)

    mesh.write(name, io='auto')

    return mesh

    def test_refine_3_8(self):
        mesh = Mesh.from_file(data_dir + '/meshes/elements/3_8_1.mesh')
        domain = refine(FEDomain('domain', mesh), self.options.out_dir, 1)

        ok = compare_mesh('3_8', domain.mesh.coors, domain.mesh.get_conn('3_8'))

        return ok

def save_basis(nurbs, pars):
    """
    Save a NURBS object basis on a FE mesh corresponding to the given
    parametrization in VTK files.

    Parameters
    ----------
    nurbs : igakit.nurbs.NURBS instance
        The NURBS object.
    pars : sequence of array, optional
        The values of parameters in each parametric dimension.
    """
    coors, conn, desc = create_linear_fe_mesh(nurbs, pars)
    mat_id = nm.zeros(conn.shape[0], dtype=nm.int32)
    mesh = Mesh.from_data('nurbs', coors, None, [conn], [mat_id], [desc])

    n_dof = nurbs.weights.ravel().shape[0]
    variable = nm.zeros(n_dof, dtype=nm.float64)
    field = variable.reshape(nurbs.weights.shape)
    for ic in xrange(n_dof):
        variable[ic - 1] = 0.0
        variable[ic] = 1.0

        vals = nurbs.evaluate(field, *pars).reshape((-1))
        out = {}
        out['bf'] = Struct(name='output_data', mode='vertex',
                           data=vals[:, None])
        mesh.write('iga_basis_%03d.vtk' % ic, io='auto', out=out)

def main():
    parser = ArgumentParser(description=__doc__.rstrip(),
                            formatter_class=RawDescriptionHelpFormatter)
    parser.add_argument('filename', help=helps['filename'])
    parser.add_argument('-d', '--detailed',
                        action='store_true', dest='detailed',
                        default=False, help=helps['detailed'])
    options = parser.parse_args()

    mesh = Mesh.from_file(options.filename)

    output(mesh.cmesh)
    output('element types:', mesh.descs)
    output('nodal BCs:', sorted(mesh.nodal_bcs.keys()))

    bbox = mesh.get_bounding_box()
    output('bounding box: %s'
           % ', '.join('%s: [%s, %s]' % (name, bbox[0, ii], bbox[1, ii])
                       for ii, name in enumerate('xyz'[:mesh.dim])))

    output('centre:', mesh.coors.mean(0))

    if not options.detailed: return

    from sfepy.discrete.fem.geometry_element import create_geometry_elements
    gels = create_geometry_elements()
    mesh.cmesh.set_local_entities(gels)
    mesh.cmesh.setup_entities()

    for dim in range(1, mesh.cmesh.tdim + 1):
        volumes = mesh.cmesh.get_volumes(dim)
        output('volumes of %d %dD entities: min: %s mean: %s max: %s'
               % (mesh.cmesh.num[dim],
                  dim, volumes.min(), volumes.mean(), volumes.max()))

def triangulate(mesh, verbose=False):
    """
    Triangulate a 2D or 3D tensor product mesh: quadrilaterals->triangles,
    hexahedrons->tetrahedrons.

    Parameters
    ----------
    mesh : Mesh
        The input mesh.

    Returns
    -------
    mesh : Mesh
        The triangulated mesh.
    """
    conns = None
    for k, new_desc in [('3_8', '3_4'), ('2_4', '2_3')]:
        if k in mesh.descs:
            conns = mesh.get_conn(k)
            break

    if conns is not None:
        nelo = conns.shape[0]
        output('initial mesh: %d elements' % nelo, verbose=verbose)

        new_conns = elems_q2t(conns)
        nn = new_conns.shape[0] // nelo
        new_cgroups = nm.repeat(mesh.cmesh.cell_groups, nn)

        output('new mesh: %d elements' % new_conns.shape[0], verbose=verbose)
        mesh = Mesh.from_data(mesh.name, mesh.coors,
                              mesh.cmesh.vertex_groups,
                              [new_conns], [new_cgroups], [new_desc])

    return mesh

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.discrete.fem import Mesh, FEDomain

    if level > 0:
        mesh = Mesh.from_file(filename)
        domain = FEDomain(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