def main():
    parser = OptionParser(usage=usage, version="%prog")

    options, args = parser.parse_args()

    if (len(args) == 1):
        mesh_filename = args[0];
    else:
        parser.print_help(),
        return

    mesh = Mesh('mesh', mesh_filename)
    print mesh

    domain = Domain('domain', mesh)
    print domain

    reg = domain.create_region('Surface',
                               'nodes of surface',
                               {'can_cells' : True})

    dual_mesh = DualMesh(reg)
    dual_mesh.save('dual_mesh.mesh',)
    dual_mesh.save_axes('axes.vtk',)

    print dual_mesh

    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 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):
        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 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 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 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 do_interpolation(m2, m1, data, field_name, force=False):
    """Interpolate data from m1 to m2. """
    from sfepy.fem import Domain, H1NodalVolumeField, Variables

    fields = {
        'scalar_si' : ((1,1), 'Omega', 2),
        'vector_si' : ((3,1), 'Omega', 2),
        'scalar_tp' : ((1,1), 'Omega', 1),
        'vector_tp' : ((3,1), 'Omega', 1),
    }

    d1 = Domain('d1', m1)

    omega1 = d1.create_region('Omega', 'all')

    f = fields[field_name]

    field1 = H1NodalVolumeField('f', nm.float64, f[0], d1.regions[f[1]],
                                approx_order=f[2])
    ff = {field1.name : field1}

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

    d2 = Domain('d2', m2)
    omega2 = d2.create_region('Omega', 'all')

    field2 = H1NodalVolumeField('f', nm.float64, f[0], d2.regions[f[1]],
                                approx_order=f[2])
    ff2 = {field2.name : field2}

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

    if not force:
        # 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.5)

    else:
        coors = u2.field.get_coor()
        vals = u1.evaluate_at(coors, close_limit=0.5)
        u2.set_data(vals)

    return u1, u2

    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.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

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

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

        coors = self.field.get_coor()
        vals = nm.sin(2.0 * nm.pi * coors[:,0] * coors[:,1])
        source.data_from_any(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 = Domain('domain', mesh)

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


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

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

        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 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', 'nodes of surface')

            domains.append(domain)

        integrals = {'Omega' : Integral('iv', kind='v', order=3),
                     'Gamma' : Integral('is', kind='s', order=3)}

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

def mesh_hook(mesh, mode):
    """
    Load and refine a mesh here.
    """
    if mode == 'read':
        mesh = Mesh.from_file(base_mesh)
        domain = Domain(mesh.name, mesh)
        for ii in range(3):
            output('refine %d...' % ii)
            domain = domain.refine()
            output('... %d nodes %d elements'
                   % (domain.shape.n_nod, domain.shape.n_el))

        domain.mesh.name = '2_4_2_refined'

        return domain.mesh

    elif mode == 'write':
        pass

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

        source = FieldVariable("us", "unknown", self.field, 1)

        coors = self.field.get_coor()
        vals = nm.sin(2.0 * nm.pi * coors[:, 0] * coors[:, 1])
        source.data_from_any(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 = Domain("domain", mesh)

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

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

        target = FieldVariable("ut", "unknown", field, 1)

        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

import numpy as np

from sfepy.fem import (Mesh, Domain, Field,
                       FieldVariable,
                       Material, Integral,
                       Equation, Equations,
                       ProblemDefinition)
from sfepy.terms import Term
from sfepy.fem.conditions import Conditions, EssentialBC
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.postprocess import Viewer

mesh = Mesh.from_file('meshes/2d/square_tri2.mesh')
domain = Domain('domain', mesh)

omega = domain.create_region('Omega', 'all')
left = domain.create_region('Left',
                            'vertices in x < -0.999',
                            'facet')
right = domain.create_region('Right',
                             'vertices in x > 0.999',
                             'facet')
bottom = domain.create_region('Bottom',
                              'vertices in y < -0.999',
                              'facet')
top = domain.create_region('Top',
                           'vertices in y > 0.999',
                           'facet')

domain.save_regions_as_groups('regions.vtk')

#.........这里部分代码省略.........
        for ip in get_dofs(options.dofs, ps.n_nod):
            output('shape function %d...' % ip)

            def eval_dofs(iels, rx):
                if options.derivative == 0:
                    bf = ps.eval_base(rx).squeeze()
                    rvals = bf[None, :, ip:ip+1]

                else:
                    bfg = ps.eval_base(rx, diff=True)
                    rvals = bfg[None, ..., ip]

                return rvals

            def eval_coors(iels, rx):
                bf = gps.eval_base(rx).squeeze()
                coors = nm.dot(bf, gel.coors)[None, ...]
                return coors

            (level, coors, conn,
             vdofs, mat_ids) = create_output(eval_dofs, eval_coors, 1,
                                             ps, min_level=lin.min_level,
                                             max_level=lin.max_level,
                                             eps=lin.eps)
            out = {
                'bf' : Struct(name='output_data',
                              mode='vertex', data=vdofs,
                              var_name='bf', dofs=None)
            }

            mesh = Mesh.from_data('bf_mesh', coors, None, [conn], [mat_ids],
                                  [options.geometry])

            name = name_template % ip
            mesh.write(name, out=out)

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

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

        domain = Domain('domain', mesh)

        if options.permutations:
            permutations = [int(ii) for ii in options.permutations.split(',')]
            output('using connectivity permutations:', permutations)
            for group in domain.iter_groups():
                perms = group.gel.get_conn_permutations()[permutations]
                offsets = nm.arange(group.shape.n_el) * group.shape.n_ep

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

            domain.setup_facets()

        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, 1)

        if options.plot_dofs:
            import sfepy.postprocess.plot_dofs as pd
            group = domain.groups[0]
            ax = pd.plot_mesh(None, mesh.coors, mesh.conns[0], group.gel.edges)
            ax = pd.plot_global_dofs(ax, field.get_coor(), field.aps[0].econn)
            ax = pd.plot_local_dofs(ax, field.get_coor(), field.aps[0].econn)
            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.vtk' % _format)
        for ip in get_dofs(options.dofs, var.n_dof):
            output('dof %d...' % ip)

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

            var.data_from_any(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
            out['u'].mesh.write(name, out=out)

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

def _gen_common_data(order, gels, report):
    import sfepy
    from sfepy.base.base import Struct
    from sfepy.linalg import combine
    from sfepy.fem import Mesh, Domain, Field, FieldVariable, Integral
    from sfepy.fem.global_interp import get_ref_coors

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

    for geom, poly_space_base in combine([['2_4', '3_8'],
                                          ['lagrange', 'lobatto']]):
        report('geometry: %s, base: %s' % (geom, poly_space_base))

        mesh0 = Mesh.from_file('meshes/elements/%s_2.mesh' % geom,
                               prefix_dir=sfepy.data_dir)
        gel = gels[geom]

        perms = gel.get_conn_permutations()

        qps, qp_weights = integral.get_qp(gel.surface_facet.name)
        zz = nm.zeros_like(qps[:, :1])
        qps = nm.hstack(([qps] + [zz]))

        rot = rots[geom]
        if rot is not None:
            pass

        shift = shifts[geom]
        rcoors = nm.ascontiguousarray(qps
                                      + shift[:1, :] - shift[1:, :])
        ccoors = nm.ascontiguousarray(qps
                                      + shift[:1, :] + shift[1:, :])

        for ir, pr in enumerate(perms):
            for ic, pc in enumerate(perms):
                report('ir: %d, ic: %d' % (ir, ic))

                mesh = mesh0.copy()
                conn = mesh.conns[0]
                conn[0, :] = conn[0, pr]
                conn[1, :] = conn[1, pc]

                cache = Struct(mesh=mesh)

                domain = Domain('domain', mesh)
                omega = domain.create_region('Omega', 'all')
                region = domain.create_region('Facet', rsels[geom])
                field = Field.from_args('f', nm.float64, shape=1,
                                        region=omega, approx_order=order,
                                        poly_space_base=poly_space_base)
                var = FieldVariable('u', 'unknown', field, 1)
                report('# dofs: %d' % var.n_dof)

                vec = nm.empty(var.n_dof, dtype=var.dtype)

                ap = field.aps[0]
                ps = ap.interp.poly_spaces['v']

                dofs = field.get_dofs_in_region_group(region, 0,
                                                      merge=False)
                edofs, fdofs = nm.unique(dofs[1]), nm.unique(dofs[2])

                rrc, rcells, rstatus = get_ref_coors(field, rcoors,
                                                     cache=cache)
                crc, ccells, cstatus = get_ref_coors(field, ccoors,
                                                     cache=cache)

                yield (geom, poly_space_base, qp_weights, mesh, ir, ic,
                       ap, ps, rrc, crc, vec, edofs, fdofs)

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('-c', '--center', metavar='center',
                      action='store', dest='center',
                      default=None, help=help['center'])
    parser.add_option('-r', '--refine', metavar='level',
                      action='store', type=int, dest='refine',
                      default=0, help=help['refine'])
    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 = _parse_val_or_vec(options.scale, 'scale', parser)
    center = _parse_val_or_vec(options.center, 'center', parser)

    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)

    if center is not None:
        cc = 0.5 * mesh.get_bounding_box().sum(0)
        shift = center - cc
        tr = nm.c_[nm.eye(mesh.dim, dtype=nm.float64), shift[:, None]]
        mesh.transform_coors(tr)

    if options.refine > 0:
        domain = Domain(mesh.name, mesh)
        output('initial mesh: %d nodes %d elements'
               % (domain.shape.n_nod, domain.shape.n_el))

        for ii in range(options.refine):
            output('refine %d...' % ii)
            domain = domain.refine()
            output('... %d nodes %d elements'
                   % (domain.shape.n_nod, domain.shape.n_el))

        mesh = domain.mesh

    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 main():
    parser = OptionParser( usage = usage, version = "%prog " + version )
    parser.add_option( "-m", "--mesh",
                       action = "store_true", dest = "save_mesh",
                       default = True,
                       help = "save surface mesh [default: %default]" )
    parser.add_option( "-n", "--no-surface",
                       action = "store_true", dest = "no_surface",
                       default = False,
                       help = "do not output surface [default: %default]" )
    (options, args) = parser.parse_args()

    if (len( args ) == 2):
        filename_in = args[0];
        filename_out = args[1];
    else:
        parser.print_help(),
        return

    if (filename_in == '-'):
        file_in = sys.stdin
    else:
        file_in = open( filename_in, "r" ); 

    mesh = Mesh.from_file( filename_in )

    if (filename_in != '-'):
        file_in.close()

    domain = Domain.from_mesh( mesh, op.join( init_sfepy.install_dir, 'eldesc' ) )
    domain.setup_groups()

    if domain.has_faces():
        domain.fix_element_orientation()
        domain.setup_neighbour_lists( create_edge_list = False )

        lst, surf_faces = domain.surface_faces()

        surf_mesh = Mesh.from_surface( surf_faces, mesh )

        if options.save_mesh:
            base, ext = op.splitext( op.basename( filename_in ) )
            surf_mesh.write( "surf_" + base + '.mesh', io = 'auto' )

        if options.no_surface:
            return

        n_nod = mesh.nod0.shape[0]
        gr_s = surface_graph( surf_faces, n_nod )
##         import sfepy.base.plotutils as plu
##         plu.spy( gr_s )
##         plu.pylab.show()

        n_comp, comps = surface_components( gr_s, surf_faces )
#        print 'components:', n_comp

        ccs, comps = comps, nm.zeros( (0,1), nm.int32 )
        for cc in ccs:
            comps = nm.concatenate( (comps, cc[:,nm.newaxis]), 0 )

        out = nm.concatenate( (lst, comps), 1 )

        if (filename_out == '-'):
            file_out = sys.stdout
        else:
            file_out = open( filename_out, "w" ); 
        for row in out:
            file_out.write( '%d %d %d %d\n' % (row[0], row[1], row[2], row[3]) )
        if (filename_out != '-'):
            file_out.close()