def from_conf(conf, dpb, apb):

        opts = conf.options
        regions = dpb.domain.regions
        if opts.use_mesh_velocity:
            import tables as pt

            fd = pt.open_file(opts.mesh_velocity_filename, mode="r")
            aux = fd.get_node("/u").read()
            nu = nm.asarray(aux, dtype=nm.float64)
            fd.close()

        else:
            nu = None

        sp_boxes = ffd.read_spline_box_hdf5(opts.ffd_spline_data)
        dsg_vars = ffd.read_dsg_vars_hdf5(opts.ffd_spline_data)
        dsg_vars.renumber_by_boxes(sp_boxes)
        dsg_vars.normalize_null_space_base()
        print dsg_vars.indx.shape
        print dsg_vars.null_space_b.shape

        control_region = regions[opts.control_domain]
        design_region = regions[opts.design_domain]

        from sfepy.fem.mesh import Mesh

        cmm = Mesh.from_region(control_region, dpb.domain.mesh)
        dmm = Mesh.from_region(design_region, dpb.domain.mesh)
        cmm.write("control.mesh", io="auto")
        dmm.write("design.mesh", io="auto")

        SOFC = ShapeOptimFlowCase
        obj = SOFC(
            dpb=dpb,
            apb=apb,
            sp_boxes=sp_boxes,
            dsg_vars=dsg_vars,
            problem_type=opts.problem,
            objective_function_type=opts.objective_function,
            var_map=opts.var_map,
            nu=nu,
            use_mesh_velocity=opts.use_mesh_velocity,
            save_dir=opts.save_dir,
            save_iter_sols=opts.save_iter_sols,
            save_control_points=opts.save_control_points,
            save_dsg_vars=opts.save_dsg_vars,
            test_terms_if_test=opts.test_terms_if_test,
        )

        equations = getattr(conf, "_".join(("equations_sensitivity", opts.problem, opts.objective_function)))

        obj.obj_fun_term = equations["objective"]
        obj.sens_terms = equations["sensitivity"]

        obj.n_var = dsg_vars.val.shape[0]

        obj.create_evaluables()

        return obj

    def create_mesh(self, extra_nodes=True):
        """
        Create a mesh from the field region, optionally including the field
        extra nodes.
        """
        mesh = self.domain.mesh

        if self.approx_order != 0:
            conns, mat_ids, descs = [], [], []
            for ig, ap in self.aps.iteritems():
                group = self.domain.groups[ig]
                if extra_nodes:
                    conn = ap.econn
                else:
                    offset = group.shape.n_ep
                    conn = ap.econn[:,:offset]
                conns.append(conn)
                mat_ids.append(mesh.mat_ids[ig])
                descs.append(mesh.descs[ig])

            if extra_nodes:
                coors = self.coors

            else:
                coors = self.coors[:self.n_vertex_dof]

            mesh = Mesh.from_data(self.name, coors, None, conns,
                                  mat_ids, descs)

        return mesh

def dump_to_vtk(filename, output_filename_trunk=None, step0=0, steps=None,
                fields=None, linearization=None):
    """Dump a multi-time-step results file into a sequence of VTK files."""
    def _save_step(suffix, out, mesh):
        if linearization is not None:
            output('linearizing...')
            out = _linearize(out, fields, linearization)
            output('...done')
            for key, val in out.iteritems():
                lmesh = val.get('mesh', mesh)
                lmesh.write(output_filename_trunk + '_' + key + suffix,
                            io='auto', out={key : val})
                if hasattr(val, 'levels'):
                    output('max. refinement per group:', val.levels)

        else:
            mesh.write(output_filename_trunk + suffix, io='auto', out=out)

    output('dumping to VTK...')

    io = MeshIO.any_from_filename(filename)
    mesh = Mesh.from_file(filename, io=io)

    if output_filename_trunk is None:
        output_filename_trunk = get_trunk(filename)

    try:
        ts = TimeStepper(*io.read_time_stepper())
        times, nts, dts = extract_times(filename)

    except ValueError:
        output('no time stepping info found, assuming single step')

        out = io.read_data(0)
        if out is not None:
            _save_step('.vtk', out, mesh)

        ret = None

    else:
        ts.times = times
        ts.n_step = times.shape[0]

        if steps is None:
            iterator = ts.iter_from(step0)

        else:
            iterator = [(step, ts.times[step]) for step in steps]

        for step, time in iterator:
            output(ts.format % (step, ts.n_step - 1))
            out = io.read_data(step)
            if out is None: break

            _save_step('.' + ts.suffix % step + '.vtk', out, mesh)

        ret = ts.suffix

    output('...done')
    return ret

def dump_to_vtk( filename, options, steps = None ):
    output( 'dumping to VTK...' )
    
    mesh = Mesh.from_file( filename )

    io = HDF5MeshIO( filename )
    ts = TimeStepper( *io.read_time_stepper() )

    if options.output_filename_trunk:
        ofn_trunk = options.output_filename_trunk
    else:
        ofn_trunk = get_trunk( filename )

    if steps is None:
        iterator = ts.iter_from( options.step0 )
    else:
        iterator = [(step, ts.times[step]) for step in steps]

    for step, time in iterator:
        output( ts.format % (step, ts.n_step - 1) )
        out = io.read_data( step )
        if out is None: break
        mesh.write( ofn_trunk + ts.suffix % step + '.vtk',
                    io = 'auto', out = out )

    output( '...done' )
    return ts.suffix

def main():
    parser = OptionParser(usage=usage, version="%prog 42")
    parser.add_option("-s", "--scale", type=int, metavar='scale',
                      action="store", dest="scale",
                      default=2, help=help['scale'])
    parser.add_option("-r", "--repeat", type='str', metavar='nx,ny[,nz]',
                      action="callback", dest="repeat",
                      callback=parse_repeat, default=None, help=help['repeat'])
    parser.add_option("-e", "--eps", type=float, metavar='eps',
                      action="store", dest="eps",
                      default=1e-8, help=help['eps'])
    (options, args) = parser.parse_args()

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

    output = Output('genPerMesh:')
    output('scale:', options.scale)
    output('repeat:', options.repeat)
    output('eps:', options.eps)

    mesh_in = Mesh.from_file(filename_in)
    mesh_out = gen_tiled_mesh(mesh_in, options.repeat, 1./options.scale,
                              options.eps)
    mesh_out.write(filename_out, io='auto')
    output('done.')

def gen_block_mesh(dims, shape, centre, mat_id=0, name='block',
                   coors=None, verbose=True):
    """
    Generate a 2D or 3D block mesh. The dimension is determined by the
    lenght of the shape argument.

    Parameters
    ----------
    dims : array of 2 or 3 floats
        Dimensions of the block.
    shape : array of 2 or 3 ints
        Shape (counts of nodes in x, y, z) of the block mesh.
    centre : array of 2 or 3 floats
        Centre of the block.
    mat_id : int, optional
        The material id of all elements.
    name : string
        Mesh name.
    verbose : bool
        If True, show progress of the mesh generation.

    Returns
    -------
    mesh : Mesh instance
    """
    dims = nm.asarray(dims, dtype=nm.float64)
    shape = nm.asarray(shape, dtype=nm.int32)
    centre = nm.asarray(centre, dtype=nm.float64)

    dim = shape.shape[0]

    centre = centre[:dim]
    dims = dims[:dim]

    n_nod = nm.prod(shape)
    output('generating %d vertices...' % n_nod, verbose=verbose)

    x0 = centre - 0.5 * dims
    dd = dims / (shape - 1)

    ngrid = nm.mgrid[[slice(ii) for ii in shape]]
    ngrid.shape = (dim, n_nod)

    coors = x0 + ngrid.T * dd
    output('...done', verbose=verbose)

    n_el = nm.prod(shape - 1)
    output('generating %d cells...' % n_el, verbose=verbose)

    mat_ids = nm.empty((n_el,), dtype=nm.int32)
    mat_ids.fill(mat_id)

    conn, desc = get_tensor_product_conn(shape)
    output('...done', verbose=verbose)

    mesh = Mesh.from_data(name, coors, None, [conn], [mat_ids], [desc])
    return mesh

    def test_mesh_smoothing(self):
        from sfepy.mesh.mesh_tools import smooth_mesh
        from sfepy.fem.mesh import Mesh
        from sfepy import data_dir

        mesh = Mesh.from_file(data_dir + '/meshes/3d/cylinder.vtk')
        vol0 = get_volume(mesh.conns[0], mesh.coors)
        mesh.coors = smooth_mesh(mesh, n_iter=10)
        vol1 = get_volume(mesh.conns[0], mesh.coors)
        filename = op.join(self.options.out_dir, 'smoothed_cylinder.vtk')
        mesh.write(filename)
        frac = vol1 / vol0

        if (frac < 0.967) and (frac > 0.966):
            self.report('mesh smoothed')
            return True

        else:
            self.report('mesh smoothed, volume mismatch!')
            return False

def dump_to_vtk(filename, output_filename_trunk=None, step0=0, steps=None):
    """Dump a multi-time-step results file into a sequence of VTK files."""
    output('dumping to VTK...')
    
    io = MeshIO.any_from_filename(filename)
    mesh = Mesh.from_file(filename, io=io)

    if output_filename_trunk is None:
        output_filename_trunk = get_trunk(filename)

    try:
        ts = TimeStepper(*io.read_time_stepper())

    except:
        output('no time stepping info found, assuming single step')

        out = io.read_data(0)
        if out is not None:
            mesh.write(output_filename_trunk + '.vtk', io='auto', out=out)

        ret = None

    else:
        if steps is None:
            iterator = ts.iter_from(step0)

        else:
            iterator = [(step, ts.times[step]) for step in steps]

        for step, time in iterator:
            output(ts.format % (step, ts.n_step - 1))
            out = io.read_data(step)
            if out is None: break
            mesh.write('.'.join((output_filename_trunk,
                                 ts.suffix % step, 'vtk')),
                       io='auto', out=out)

        ret = ts.suffix

    output('...done')
    return ret

    def create_mesh_from_control_points( self ):
        offset = 0
        dim = self.spbs[0].cxyz.shape[1]
        coors = nm.empty((0, dim), dtype=nm.float64)
        conns = []
        mat_ids = []
        descs = []
        for ib, spb in enumerate( self.spbs ):
            n_nod = spb.cxyz.shape[0]
            coors = nm.concatenate( (coors, spb.cxyz), 0 )
            descs.append( '3_2' )

            conn = []
            for ij in xrange( spb.cpi.shape[1] ):
                for ik in xrange( spb.cpi.shape[2] ):
                    inx = spb.cpi[:,ij,ik]
                    row = [[p1, p2] for p1, p2 in zip( inx[:-1], inx[1:] )]
                    conn.extend( row )
            for ij in xrange( spb.cpi.shape[0] ):
                for ik in xrange( spb.cpi.shape[2] ):
                    inx = spb.cpi[ij,:,ik]
                    row = [[p1, p2] for p1, p2 in zip( inx[:-1], inx[1:] )]
                    conn.extend( row )
            for ij in xrange( spb.cpi.shape[0] ):
                for ik in xrange( spb.cpi.shape[1] ):
                    inx = spb.cpi[ij,ik,:]
                    row = [[p1, p2] for p1, p2 in zip( inx[:-1], inx[1:] )]
                    conn.extend( row )

            aux = nm.empty(len(conn), dtype=nm.int32)
            aux.fill(ib)
            mat_ids.append(aux)

            conns.append( offset + nm.array( conn, dtype = nm.int32 ) )
            offset += n_nod

        mesh = Mesh.from_data('control_points', coors, None, conns,
                              mat_ids, descs)
        return mesh

def main():
    parser = OptionParser(usage=usage, version="%prog 42")
    parser.add_option(
        "-s", "--scale", type=int, metavar="scale", action="store", dest="scale", default=2, help=help["scale"]
    )
    parser.add_option(
        "-r",
        "--repeat",
        type="str",
        metavar="nx,ny[,nz]",
        action="callback",
        dest="repeat",
        callback=parse_repeat,
        default=None,
        help=help["repeat"],
    )
    parser.add_option(
        "-e", "--eps", type=float, metavar="eps", action="store", dest="eps", default=1e-8, help=help["eps"]
    )
    parser.add_option("-n", "--no-mvd", action="store_true", dest="nomvd", default=False, help=help["nomvd"])
    (options, args) = parser.parse_args()

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

    output = Output("genPerMesh:")
    output("scale:", options.scale)
    output("repeat:", options.repeat)
    output("eps:", options.eps)

    mesh_in = Mesh.from_file(filename_in)
    mesh_out = compose_periodic_mesh(mesh_in, options.scale, options.repeat, options.eps, check_mvd=not options.nomvd)
    mesh_out.write(filename_out, io="auto")
    output("done.")

def solve_eigen_problem1( conf, options ):


    pb = ProblemDefinition.from_conf( conf )
    dim = pb.domain.mesh.dim

    pb.time_update()

    dummy = pb.create_state_vector()

    output( 'assembling lhs...' )
    tt = time.clock()
    mtx_a = eval_term_op( dummy, conf.equations['lhs'], pb,
                       dw_mode = 'matrix', tangent_matrix = pb.mtx_a )
    output( '...done in %.2f s' % (time.clock() - tt) )

    output( 'assembling rhs...' )
    tt = time.clock()
    mtx_b = eval_term_op( dummy, conf.equations['rhs'], pb,
                       dw_mode = 'matrix', tangent_matrix = pb.mtx_a.copy() )
    output( '...done in %.2f s' % (time.clock() - tt) )

    #mtxA.save( 'tmp/a.txt', format='%d %d %.12f\n' )
    #mtxB.save( 'tmp/b.txt', format='%d %d %.12f\n' )
    try:
        n_eigs = conf.options.n_eigs
    except AttributeError:
        n_eigs = mtx_a.shape[0]

    if n_eigs is None:
        n_eigs = mtx_a.shape[0]

##     mtx_a.save( 'a.txt', format='%d %d %.12f\n' )
##     mtx_b.save( 'b.txt', format='%d %d %.12f\n' )
    print 'computing resonance frequencies...'
    eig = Solver.any_from_conf( pb.get_solver_conf( conf.options.eigen_solver ) )
    eigs, mtx_s_phi = eig( mtx_a, mtx_b, conf.options.n_eigs )
    from sfepy.fem.mesh import Mesh
    bounding_box = Mesh.from_file("tmp/mesh.vtk").get_bounding_box()
    # this assumes a box (3D), or a square (2D):
    a = bounding_box[1][0] - bounding_box[0][0]
    E_exact = None
    if options.hydrogen or options.boron:
        if options.hydrogen:
            Z = 1
        elif options.boron:
            Z = 5
        if options.dim == 2:
            E_exact = [-float(Z)**2/2/(n-0.5)**2/4 for n in [1]+[2]*3+[3]*5 +\
                    [4]*8 + [5]*15]
        elif options.dim == 3:
            E_exact = [-float(Z)**2/2/n**2 for n in [1]+[2]*2**2+[3]*3**2 ]
    if options.well:
        if options.dim == 2:
            E_exact = [pi**2/(2*a**2)*x for x in [2, 5, 5, 8, 10, 10, 13, 13,
                17, 17, 18, 20, 20 ] ]
        elif options.dim == 3:
            E_exact = [pi**2/(2*a**2)*x for x in [3, 6, 6, 6, 9, 9, 9, 11, 11,
                11, 12, 14, 14, 14, 14, 14, 14, 17, 17, 17] ]
    if options.oscillator:
        if options.dim == 2:
            E_exact = [1] + [2]*2 + [3]*3 + [4]*4 + [5]*5 + [6]*6
        elif options.dim == 3:
            E_exact = [float(1)/2+x for x in [1]+[2]*3+[3]*6+[4]*10 ]
    if E_exact is not None:
        print "a=%f" % a
        print "Energies:"
        print     "n      exact         FEM      error"

        for i, e in enumerate(eigs):
            from numpy import NaN
            if i < len(E_exact):
                exact = E_exact[i]
                err = 100*abs((exact - e)/exact)
            else:
                exact = NaN
                err = NaN
            print "%d:  %.8f   %.8f  %5.2f%%" % (i, exact, e, err)
    else:
        print eigs
##     import sfepy.base.plotutils as plu
##     plu.spy( mtx_b, eps = 1e-12 )
##     plu.pylab.show()
##     pause()
    n_eigs = eigs.shape[0]

    mtx_phi = nm.empty( (pb.variables.di.ptr[-1], mtx_s_phi.shape[1]),
                       dtype = nm.float64 )
    for ii in xrange( n_eigs ):
        mtx_phi[:,ii] = pb.variables.make_full_vec( mtx_s_phi[:,ii] )

    save = get_default_attr( conf.options, 'save_eig_vectors', None )
    out = {}
    for ii in xrange( n_eigs ):
        if save is not None:
            if (ii > save[0]) and (ii < (n_eigs - save[1])): continue
        aux = pb.state_to_output( mtx_phi[:,ii] )
        key = aux.keys()[0]
        out[key+'%03d' % ii] = aux[key]

#.........这里部分代码省略.........

def extract_time_history( filename, options ):
    output( 'extracting selected data...' )

    el = options.extract_list
    output( 'extraction list:', el )

    ##
    # Parse extractions.
    pes = OneTypeList( Struct )
    for chunk in el.split( ',' ):
        aux =  chunk.strip().split()
        pes.append( Struct( var = aux[0],
                            mode = aux[1],
                            indx = map( int, aux[2:] ),
                            igs = None ) )

    ##
    # Verify array limits, set igs for element data, shift indx.
    mesh = Mesh.from_file( filename )
    n_el, n_els, offs = mesh.n_el, mesh.n_els, mesh.el_offsets
    for pe in pes:
        if pe.mode == 'n':
            for ii in pe.indx:
                if (ii < 0) or (ii >= mesh.n_nod):
                    raise IndexError, 'node index 0 <= %d < %d'\
                          % (ii, mesh.n_nod)

        if pe.mode == 'e':
            pe.igs = []
            for ii, ie in enumerate( pe.indx[:] ):
                if (ie < 0) or (ie >= n_el):
                    raise IndexError, 'element index 0 <= %d < %d'\
                          % (ie, n_el)
                ig = (ie < n_els).argmax()
                pe.igs.append( ig )
                pe.indx[ii] = ie - offs[ig]

##     print pes

    ##
    # Extract data.
    # Assumes only one element group (ignores igs)!
    io = HDF5MeshIO( filename )
    ths = {}
    for pe in pes:
        mode, nname = io.read_data_header( pe.var )
        print mode, nname
        if ((pe.mode == 'n' and mode == 'vertex') or
            (pe.mode == 'e' and mode == 'cell')):
            th = io.read_time_history( nname, pe.indx )

        elif pe.mode == 'e' and mode == 'vertex':
            conn = mesh.conns[0]
            th = {}
            for iel in pe.indx:
                ips = conn[iel]
                th[iel] = io.read_time_history( nname, ips )
        else:
            raise RuntimeError, 'cannot extract cell data %s in nodes' % pe.var
            
        ths[pe.var] = th
    return ths

def main():
    parser = OptionParser(usage=usage, version='%prog ' + sfepy.__version__)
    parser.add_option('-c', '--conf', metavar='"key : value, ..."',
                      action='store', dest='conf', type='string',
                      default=None, help= help['conf'])
    parser.add_option('-O', '--options', metavar='"key : value, ..."',
                      action='store', dest='app_options', type='string',
                      default=None, help=help['options'])
    parser.add_option('-o', '', metavar='filename',
                      action='store', dest='output_filename_trunk',
                      default=None, help=help['filename'])
    parser.add_option('--create-mesh',
                      action='store_true', dest='create_mesh',
                      default=False, help=help['create_mesh'])
    parser.add_option('--2d',
                      action='store_true', dest='dim2',
                      default=False, help=help['dim'])
    parser.add_option('-m', '--mesh', metavar='filename',
                      action='store', dest='mesh',
                      default=None, help=help['mesh'])
    parser.add_option('--mesh-dir', metavar='dirname',
                      action='store', dest='mesh_dir',
                      default='tmp', help=help['mesh_dir'])
    parser.add_option('--oscillator',
                      action='store_true', dest='oscillator',
                      default=False, help=help['oscillator'])
    parser.add_option('--well',
                      action='store_true', dest='well',
                      default=False, help=help['well'])
    parser.add_option('--hydrogen',
                      action='store_true', dest='hydrogen',
                      default=False, help=help['hydrogen'])
    parser.add_option('--boron',
                      action='store_true', dest='boron',
                      default=False, help=help['boron'])

    options, args = parser.parse_args()

    if options.create_mesh and options.mesh:
        output('--create-mesh and --mesh options are mutually exclusive!')
        return

    if len(args) == 1:
        filename_in = args[0];
        auto_mesh_name = False

    elif len(args) == 0:
        auto_mesh_name = True

        mesh_filename = os.path.join(options.mesh_dir, 'mesh.vtk')
        ensure_path(mesh_filename)

        if options.oscillator:
            filename_in = fix_path("examples/quantum/oscillator.py")

        elif options.well:
            filename_in = fix_path("examples/quantum/well.py")

        elif options.hydrogen:
            filename_in = fix_path("examples/quantum/hydrogen.py")

        elif options.boron:
            filename_in = fix_path("examples/quantum/boron.py")

        elif options.create_mesh:
            output('generating mesh...')
            try:
                os.makedirs("tmp")
            except OSError, e:
                if e.errno != 17: # [Errno 17] File exists
                    raise
            if options.dim2:
                output("dimension: 2")
                gp = fix_path('meshes/quantum/square.geo')
                os.system("cp %s tmp/mesh.geo" % gp)
                os.system("gmsh -2 tmp/mesh.geo -format mesh")
                mtv = fix_path('script/mesh_to_vtk.py')
                os.system("%s tmp/mesh.mesh %s" % (mtv, mesh_filename))
            else:
                output("dimension: 3")
                import sfepy.geom as geom
                from sfepy.fem.mesh import Mesh
                try:
                    from site_cfg import tetgen_path
                except ImportError:
                    tetgen_path = '/usr/bin/tetgen'
                gp = fix_path('meshes/quantum/box.geo')
                os.system("gmsh -0 %s -o tmp/x.geo" % gp)
                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(mesh_filename, io="auto")
            output("...mesh written to %s" % mesh_filename)
            return

        else:
            parser.print_help()
#.........这里部分代码省略.........

def gen_block_mesh(dims, shape, centre, name='block'):
    """
    Generate a 2D or 3D block mesh. The dimension is determined by the
    lenght of the shape argument.

    Parameters
    ----------
    dims : array of 2 or 3 floats
        Dimensions of the block.
    shape : array of 2 or 3 ints
        Shape (counts of nodes in x, y, z) of the block mesh.
    centre : array of 2 or 3 floats
        Centre of the block.

    name : string
        Mesh name.

    Returns
    -------
    mesh : Mesh instance
    """
    dims = nm.asarray(dims, dtype=nm.float64)
    shape = nm.asarray(shape, dtype=nm.int32)
    centre = nm.asarray(centre, dtype=nm.float64)

    dim = shape.shape[0]

    centre = centre[:dim]
    dims = dims[:dim]

    x0 = centre - 0.5 * dims
    dd = dims / (shape - 1)

    grid = nm.zeros(shape, dtype = nm.int32)
    n_nod = nm.prod(shape)
    coors = nm.zeros((n_nod, dim), dtype = nm.float64)

    bar = MyBar("       nodes:")
    bar.init(n_nod)
    for ii, ic in enumerate(cycle(shape)):
        grid[tuple(ic)] = ii
        coors[ii] = x0 + ic * dd
        if not (ii % 100):
            bar.update(ii)
    bar.update(ii + 1)

    n_el = nm.prod(shape - 1)
    mat_id = nm.zeros((n_el,), dtype = nm.int32)

    if (dim == 2):
        conn = nm.zeros((n_el, 4), dtype = nm.int32)
        bar = MyBar("       elements:")
        bar.init(n_el)
        for ii, (ix, iy) in enumerate(cycle(shape - 1)):
            conn[ii,:] = [grid[ix  ,iy], grid[ix+1,iy  ],
                          grid[ix+1,iy+1], grid[ix  ,iy+1]]
            if not (ii % 100):
                bar.update(ii)
        bar.update(ii + 1)
        desc = '2_4'

    else:
        conn = nm.zeros((n_el, 8), dtype = nm.int32)
        bar = MyBar("       elements:")
        bar.init(n_el)
        for ii, (ix, iy, iz) in enumerate(cycle(shape - 1)):
            conn[ii,:] = [grid[ix  ,iy  ,iz  ], grid[ix+1,iy  ,iz  ],
                          grid[ix+1,iy+1,iz  ], grid[ix  ,iy+1,iz  ],
                          grid[ix  ,iy  ,iz+1], grid[ix+1,iy  ,iz+1],
                          grid[ix+1,iy+1,iz+1], grid[ix  ,iy+1,iz+1]]
            if not (ii % 100):
                bar.update(ii)
        bar.update(ii + 1)
        desc = '3_8'

    mesh = Mesh.from_data(name, coors, None, [conn], [mat_id], [desc])
    return mesh

def main():

    parser = OptionParser( usage = usage, version = "%prog 42" )
    parser.add_option( "-s", "--scale", type = int, metavar = 'scale',
                       action = "store", dest = "scale",
                       default = 2, help = help['scale'] )
    parser.add_option( "-e", "--eps", type = float, metavar = 'eps',
                       action = "store", dest = "eps",
                       default = 1e-8, help = help['eps'] )
    parser.add_option( "-t", "--test",
                       action = "store_true", dest = "test",
                       default = False, help = help['test'] )
    parser.add_option( "-n", "--no-mvd",
                       action = "store_true", dest = "nomvd",
                       default = False, help = help['nomvd'] )
    (options, args) = parser.parse_args()

    if options.test:
        test()
        return

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

    print 'scale:', options.scale
    print 'eps:', options.eps

    mesh_in = Mesh.from_file( filename_in )
    bbox = mesh_in.get_bounding_box()
    print 'bbox:\n', bbox
    mscale = bbox[1] - bbox[0]
    centre0 = 0.5 * (bbox[1] + bbox[0])
    print 'centre:\n', centre0

    scale = nm.array( options.scale, dtype = nm.float64 )

    # Normalize original coordinates.
    coor0 = (mesh_in.nod0[:,:-1] - centre0) / (mscale)
    dim = mesh_in.dim

    coor0, mesh_in.conns = fix_double_nodes( coor0, mesh_in.conns, options.eps )
    if not options.nomvd:
        mes0 = get_min_edge_size( coor0, mesh_in.conns )
        mvd0 = get_min_vertex_distance( coor0, mes0 )
        if mes0 > (mvd0 + options.eps):
            print '          original min. "edge" length: %.5e' % mes0
            print 'original approx. min. vertex distance: %.5e' % mvd0
            print '-> still double nodes in input mesh!'
            print 'try increasing eps...'
            raise ValueError

    for indx in cycle( [options.scale] * dim ):
        aindx = nm.array( indx, dtype = nm.float64 )
        centre = 0.5 * (2.0 * aindx - scale + 1.0)
        print indx, centre

        if aindx.sum() == 0:
            coor = coor0 + centre
            conns = mesh_in.conns
        else:
            coor1 = coor0 + centre
            conns1 = mesh_in.conns

            cmap = find_map( coor, coor1, eps = options.eps )
            if not cmap.size:
                print 'non-periodic mesh!'
#                raise ValueError
            else:
                print cmap.size / 2
            coor, conns = merge_mesh( coor, conns, coor1, conns1, cmap,
                                     eps = options.eps )

    if not options.nomvd:
        mes = get_min_edge_size( coor, conns )
        mvd = get_min_vertex_distance( coor, mes0 )

        print '          original min. "edge" length: %.5e' % mes0
        print '             final min. "edge" length: %.5e' % mes
        print 'original approx. min. vertex distance: %.5e' % mvd0
        print '   final approx. min. vertex distance: %.5e' % mvd
        if mvd < 0.99999 * mvd0:
            if mvd0 < (mes0 - options.eps):
                print '-> probably non-periodic input mesh!'
                print '   ... adjacent sides were not connected!'
                print '   try increasing eps...'
            else:
                print '-> input mesh might be periodic'
                print '   try increasing eps...'
        else:
            print '-> input mesh looks periodic'
    else:
        print 'non-periodic input mesh detection skipped!'

    print 'renormalizing...'
    coor = (coor * mscale) / scale
    print 'saving...'
#.........这里部分代码省略.........

def gen_mesh_from_poly(filename, verbose=True):
    """
    Import mesh generated by tetgen or triangle.

    Parameters
    ----------
    filename : string
        file name

    Returns
    -------
    mesh : Mesh instance
        triangular or tetrahedral mesh
    """

    def getnodes(fnods,up):
        f=file(fnods)
        l=[int(x) for x in f.readline().split()]
        npoints,dim,nattrib,nbound=l
        if verbose: up.init(npoints)
        nodes=[]
        for line in f:
            if line[0]=="#": continue
            l=[float(x) for x in line.split()]
            l = l[:(dim + 1)]
            l[0]=int(l[0])
            nodes.append(tuple(l))
            assert l[0]==len(nodes)
        assert npoints==len(nodes)
        return nodes

    def getele(fele,up):
        f=file(fele)
        l=[int(x) for x in f.readline().split()]
        nele,nnod,nattrib=l
        #we have either linear or quadratic tetrahedra:
        if nnod in [4,10]:
            elem = 'tetra'
            linear = (nnod == 4)
        if nnod in [3, 7]:
            elem = 'tri'
            linear = (nnod == 3)

        # if nattrib!=1:
        #     raise "tetgen didn't assign an entity number to each element (option -A)"
        els=[]
        regions={}
        for line in f:
            if line[0]=="#": continue
            l=[int(x) for x in line.split()]
            if elem == 'tri':
                if linear:
                    assert (len(l) - 1 - nattrib) == 3
                    els.append((l[0],l[1],l[2],l[3]))
                    regionnum=l[5]
                else:
                    assert len(l)-2 == 10
                    els.append((l[0],54,l[1],l[2],l[3],l[4],
                                l[5],l[6],l[7],l[8],l[9],l[10]))
                    regionnum=l[11]
            if elem == 'tetra':
                if linear:
                    assert len(l)-2 == 4
                    els.append((l[0],54,l[1],l[2],l[3],l[4]))
                    regionnum=l[5]
                else:
                    assert len(l)-2 == 10
                    els.append((l[0],54,l[1],l[2],l[3],l[4],
                                l[5],l[6],l[7],l[8],l[9],l[10]))
                    regionnum=l[11]
            if regionnum==0:
                print "see %s, element # %d"%(fele,l[0])
                raise "there are elements not belonging to any physical entity"
            if regions.has_key(regionnum):
                regions[regionnum].append(l[0])
            else:
                regions[regionnum]=[l[0]]
            assert l[0]==len(els)
            if verbose: up.update(l[0])
        return els,regions,linear

    def getBCfaces(ffaces,up):
        f=file(ffaces)
        l=[int(x) for x in f.readline().split()]
        nfaces,nattrib=l
        if nattrib!=1:
            raise "tetgen didn't assign an entity number to each face \
(option -A)"
        if verbose: up.init(nfaces)
        faces={}
        for line in f:
            if line[0]=="#": continue
            l=[int(x) for x in line.split()]
            assert len(l)==5
            regionnum=l[4]
            if regionnum==0: continue
            if faces.has_key(regionnum):
                faces[regionnum].append((l[1],l[2],l[3]))
            else:
                faces[regionnum]=[(l[1],l[2],l[3])]
#.........这里部分代码省略.........

def gen_cylinder_mesh(dims, shape, centre, axis='x', force_hollow=False,
                      is_open=False, open_angle=0.0, non_uniform=False,
                      name='cylinder'):
    """
    Generate a cylindrical mesh along an axis. Its cross-section can be
    ellipsoidal.

    Parameters
    ----------
    axis: one of 'x', 'y', 'z'
        The axis of the cylinder.
    dims : array of 5 floats
        Dimensions of the cylinder: inner surface semi-axes a1, b1, outer
        surface semi-axes a2, b2, length.
    shape : array of 3 ints
        Shape (counts of nodes in radial, circumferential and longitudinal
        directions) of the cylinder mesh.
    centre : array of 3 floats
        Centre of the cylinder.

    force_hollow : boolean
        Force hollow mesh even if inner radii a1 = b1 = 0.

    is_open : boolean
        Generate an open cylinder segment.
    open_angle : float
        Opening angle in radians.

    non_uniform : boolean
        If True, space the mesh nodes in radial direction so that the element
        volumes are (approximately) the same, making thus the elements towards
        the outer surface thinner.

    name : string
        Mesh name.

    Returns
    -------
    mesh : Mesh instance
    """
    dims = nm.asarray(dims, dtype=nm.float64)
    shape = nm.asarray(shape, dtype=nm.int32)
    centre = nm.asarray(centre, dtype=nm.float64)

    a1, b1, a2, b2, length = dims
    nr, nfi, nl = shape
    origin = centre - nm.array([0.5 * length, 0.0, 0.0])

    dfi = 2.0 * (nm.pi - open_angle) / nfi
    if is_open:
        nnfi = nfi + 1
    else:
        nnfi = nfi

    is_hollow = force_hollow or not (max(abs(a1), abs(b1)) < 1e-15)

    if is_hollow:
        mr = 0
    else:
        mr = (nnfi - 1) * nl

    grid = nm.zeros((nr, nnfi, nl), dtype=nm.int32)

    n_nod = nr * nnfi * nl - mr
    coors = nm.zeros((n_nod, 3), dtype=nm.float64)

    angles = nm.linspace(open_angle, open_angle+(nfi)*dfi, nfi+1)
    xs = nm.linspace(0.0, length, nl)
    if non_uniform:
        ras = nm.zeros((nr,), dtype=nm.float64)
        rbs = nm.zeros_like(ras)
        advol = (a2**2 - a1**2) / (nr - 1)
        bdvol = (b2**2 - b1**2) / (nr - 1)
        ras[0], rbs[0] = a1, b1
        for ii in range(1, nr):
            ras[ii] = nm.sqrt(advol + ras[ii-1]**2)
            rbs[ii] = nm.sqrt(bdvol + rbs[ii-1]**2)
    else:
        ras = nm.linspace(a1, a2, nr)
        rbs = nm.linspace(b1, b2, nr)

    # This is 3D only...
    bar = MyBar("       nodes:")
    bar.init(n_nod)
    ii = 0
    for ix in range(nr):
        a, b = ras[ix], rbs[ix]
        for iy, fi in enumerate(angles[:nnfi]):
            for iz, x in enumerate(xs):
                grid[ix,iy,iz] = ii
                coors[ii] = origin + [x, a * nm.cos(fi), b * nm.sin(fi)]
                if not (ii % 100):
                    bar.update(ii)
                ii += 1

                if not is_hollow and (ix == 0):
                    if iy > 0:
                        grid[ix,iy,iz] = grid[ix,0,iz]
                        ii -= 1
    print
#.........这里部分代码省略.........

def gen_mesh_from_goem(geo, a=None, quadratic=False, verbose=True,
                       refine=False, polyfilename='./meshgen.poly',
                       out='mesh', **kwargs):
    """
    Runs mesh generator - tetgen for 3D or triangle for 2D meshes.

    Parameters
    ----------
    geo : geometry
        geometry description
    a : int, optional
        a maximum area/volume constraint
    quadratic : bool, optional
        set True for quadratic elements
    verbose : bool, optional
        detailed information
    refine : bool, optional
        refines mesh

    Returns
    -------
    mesh : Mesh instance
        triangular or tetrahedral mesh
    """

    import os.path as op
    import pexpect

    # write geometry to poly file
    geo.to_poly_file(polyfilename)

    if not refine:
        params = "-Apq"
    else:
        params = "-Arq"
    if verbose:
        params = params + " -Q"
    if a != None and not refine:
        params = params + " -a%f" % (a)
    if refine:
        params = params + " -a"
    if quadratic:
        params = params + " -o2"
    params = params + " %s" % (polyfilename)

    meshgen_call = {2: 'triangle', 3: 'tetgen'}
    cmd = "%s %s" % (meshgen_call[geo.dim], params)
    if verbose: print "Generating mesh using", cmd
    if geo.dim == 2:
        p=pexpect.run(cmd, timeout=None)
        bname, ext = op.splitext(polyfilename)
        mesh = Mesh.from_file(bname + '.1.node')
        mesh.write(bname + '.' + out)
    if geo.dim == 3:
        p=pexpect.spawn(cmd, timeout=None)
        if not refine:
            p.expect("Opening %s." % (polyfilename))
        else:
            p.expect("Opening %s.node.\r\n" % (polyfilename))
            p.expect("Opening %s.ele.\r\n" % (polyfilename))
            p.expect("Opening %s.face.\r\n" % (polyfilename))
            p.expect("Opening %s.vol." % (polyfilename))
        assert p.before == ""
        p.expect(pexpect.EOF)
        if p.before != "\r\n":
            print p.before
            raise "Error when running mesh generator (see above for output): %s" % cmd