Python Mesh.from_file方法代码示例

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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.')

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
    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

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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.")

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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()
#.........这里部分代码省略.........

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
def extract_time_history(filename, extract, verbose=True):
    """Extract time history of a variable from a multi-time-step results file.

    Parameters
    ----------
    filename : str
        The name of file to extract from.
    extract : str
        The description of what to extract in a string of comma-separated
        description items. A description item consists of: name of the variable
        to extract, mode ('e' for elements, 'n' for nodes), ids of the nodes or
        elements (given by the mode). Example: 'u n 10 15, p e 0' means
        variable 'u' in nodes 10, 15 and variable 'p' in element 0.
    verbose : bool
        Verbosity control.

    Returns
    -------
    ths : dict
        The time histories in a dict with variable names as keys. If a
        nodal variable is requested in elements, its value is a dict of histories
        in the element nodes.
    ts : TimeStepper instance
        The time stepping information.
    """
    output('extracting selected data...', verbose=verbose)

    output('selection:', extract, verbose=verbose)

    ##
    # Parse extractions.
    pes = OneTypeList(Struct)
    for chunk in extract.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 ValueError('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 ValueError('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 = MeshIO.any_from_filename(filename)
    ths = {}
    for pe in pes:
        mode, nname = io.read_data_header(pe.var)
        output(mode, nname, verbose=verbose)

        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 ValueError('cannot extract cell data %s in nodes!' % pe.var)
            
        ths[pe.var] = th

    output('...done', verbose=verbose)

    ts = TimeStepper(*io.read_time_stepper())

    return ths, ts

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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...'
#.........这里部分代码省略.........

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
def main():
    version = open( op.join( init_sfepy.install_dir,
                             'VERSION' ) ).readlines()[0][:-1]

    parser = OptionParser( usage = usage, version = "%prog " + version )
    parser.add_option( "--mesh",
                       action = "store_true", dest = "mesh",
                       default = False, help = help['mesh'] )
    parser.add_option( "--2d",
                       action = "store_true", dest = "dim2",
                       default = False, help = help['dim'] )
    parser.add_option( "-o", "", metavar = 'filename',
                       action = "store", dest = "output_filename_trunk",
                       default = "mesh", help = help['filename'] )
    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'] )
    parser.add_option( "--dft",
                       action = "store_true", dest = "dft",
                       default = False, help = help['dft'] )
    parser.add_option( "-p", "--plot",
                       action = "store_true", dest = "plot",
                       default = False, help = help['plot'] )

    options, args = parser.parse_args()

    if len( args ) == 1:
        filename_in = args[0];
    elif len( args ) == 0:
        if options.oscillator:
            dim = MeshIO.any_from_filename("tmp/mesh.vtk").read_dimension()
            if dim == 2:
                filename_in = "input/quantum/oscillator2d.py"
            else:
                assert_( dim == 3 )
                filename_in = "input/quantum/oscillator3d.py"
            options.dim = dim
            print "Dimension:", dim
        elif options.well:
            dim = MeshIO.any_from_filename("tmp/mesh.vtk").read_dimension()
            if dim == 2:
                filename_in = "input/quantum/well2d.py"
            else:
                assert_( dim == 3 )
                filename_in = "input/quantum/well3d.py"
            options.dim = dim
            print "Dimension:", dim
        elif options.hydrogen:
            dim = MeshIO.any_from_filename("tmp/mesh.vtk").read_dimension()
            if dim == 2:
                filename_in = "input/quantum/hydrogen2d.py"
            else:
                assert_( dim == 3 )
                filename_in = "input/quantum/hydrogen3d.py"
            options.dim = dim
            print "Dimension:", dim
        elif options.boron:
            dim = MeshIO.any_from_filename("tmp/mesh.vtk").read_dimension()
            if dim == 2:
                filename_in = "input/quantum/boron2d.py"
            else:
                assert_( dim == 3 )
                filename_in = "input/quantum/boron3d.py"
            options.dim = dim
            print "Dimension:", dim
        elif options.mesh:
            try:
                os.makedirs("tmp")
            except OSError, e:
                if e.errno != 17: # [Errno 17] File exists
                    raise
            if options.dim2:
                print "Dimension: 2"
                os.system("cp database/quantum/square.geo tmp/mesh.geo")
                os.system("gmsh -2 tmp/mesh.geo -format mesh")
                os.system("script/mesh_to_vtk.py tmp/mesh.mesh tmp/mesh.vtk")
            else:
                print "Dimension: 3"
                import geom
                from sfepy.fem.mesh import Mesh
                try:
                    from site_cfg import tetgen_path
                except ImportError:
                    tetgen_path = '/usr/bin/tetgen'
                os.system("gmsh -0 database/box.geo -o tmp/x.geo")
                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")
#.........这里部分代码省略.........

# 需要导入模块: from sfepy.fem.mesh import Mesh [as 别名]
# 或者: from sfepy.fem.mesh.Mesh import from_file [as 别名]
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]

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