def runimport(geometryfile, iproc, filename=None):
     import start as start
     cubit = start.start_cubit()
     cfg = start.start_cfg(filename=filename)
     file1 = cfg.output_dir + '/' + geometryfile
     cubitcommand = 'open "' + file1 + '"  '
     cubit.cmd(cubitcommand)

def importgeometry(geometryfile, iproc=0, filename=None):
    import start as start
    cfg = start.start_cfg(filename=filename)
    mpiflag, iproc, numproc, mpi = start.start_mpi()

    if iproc == 0:
        print 'importing geometry....'
    a = ['ok from ' + str(iproc)]

    mpi.barrier()
    total_a = mpi.allgather(a)
    if iproc == 0:
        print total_a

    def runimport(geometryfile, iproc, filename=None):
        import start as start
        cubit = start.start_cubit()
        cfg = start.start_cfg(filename=filename)
        file1 = cfg.output_dir + '/' + geometryfile
        cubitcommand = 'open "' + file1 + '"  '
        cubit.cmd(cubitcommand)

    if cfg.parallel_import:
        runimport(geometryfile, iproc, filename=filename)
    else:
        if iproc == 0:
            runimport(geometryfile, iproc, filename=filename)
            for i in range(1, mpi.size):
                mpi.send('import', i)
                msg, status = mpi.recv(i)
        else:
            msg, status = mpi.recv(0)
            runimport(geometryfile, iproc, filename=filename)
            mpi.send('ok' + str(iproc), 0)

def savegeometry(iproc=0,surf=False,filename=None):
    import start as start
    cfg                         = start.start_cfg(filename=filename)
    mpiflag,iproc,numproc,mpi   = start.start_mpi()
    
    def runsave(geometryfile,iproc,filename=None):
        import start as start
        cubit                   = start.start_cubit()
        cfg                         = start.start_cfg(filename=filename)
        flag=[0]
        ner=cubit.get_error_count()
        cubitcommand= 'save as "'+ cfg.output_dir+'/'+geometryfile+ '"  overwrite' 
        cubit.cmd(cubitcommand)                                                    
        ner2=cubit.get_error_count()                                             
        if ner == ner2:
            flag=[1]
        return flag
        
    if surf:
        geometryfile='surf_vol_'+str(iproc)+'.cub'
    else:
        geometryfile='geometry_vol_'+str(iproc)+'.cub'
        
    flagsaved=[0]
    infosave=(iproc,flagsaved)
    
    mpi.barrier()
    total_saved=mpi.allgather(flagsaved)
    if isinstance(total_saved,int): total_saved=[total_saved]
    
    ind=0
    saving=True
    while saving:
        if len(total_saved) != sum(total_saved):
            #
            if not flagsaved[0]: 
                flagsaved=runsave(geometryfile,iproc,filename=filename)
                if flagsaved[0]:
                    infosave=(iproc,flagsaved[0])        
                    if numproc > 1:
                        f=open('geometry_saved'+str(iproc),'w')
                        f.close()
            mpi.barrier()
            total_saved=mpi.allgather(flagsaved)
            if isinstance(total_saved,int): total_saved=[total_saved]
            ind=ind+1
        else:
            saving=False
        if ind > len(total_saved)+10: saving=False
        print sum(total_saved),'/',len(total_saved),' saved'
    
    info_total_saved=mpi.allgather(infosave)
    if isinstance(info_total_saved,int): info_total_saved=[info_total_saved]
    
    if iproc==0:
        f=open('geometry_saving.log','w')
        f.write('\n'.join(str(x) for x in info_total_saved))                
        f.close()

def refinement(nvol, vol, filename=None):
    import start as start
    cfg = start.start_cfg(filename=filename)
    from utilities import get_v_h_list
    #
    # vertical refinement
    # for nvol = 3
    #
    # ___________________________ interface 4
    #
    # vol 2
    # ___________________________ interface 3
    #
    # vol 1
    # ___________________________ interface 2
    #
    # vol 0
    # ___________________________ interface 1
    #
    #
    if cfg.ntripl != 0:
        if len(cfg.refinement_depth) != 0:
            _, _, _, _, _, tsurf = get_v_h_list([vol[nvol - 1].ID])
            tsurf = ' '.join(str(x) for x in tsurf)
            for idepth in cfg.refinement_depth:
                cubitcommand = 'refine node in surf  ' + \
                    str(tsurf) + ' numsplit 1 bias 1.0 depth ' + str(idepth)
                cubit.cmd(cubitcommand)
        else:
            for ir in cfg.tripl:
                if ir == 1:
                    print ('interface = 1 means that the refinement'
                           'interface is at the bottom of the volume')
                    txt = ' all '
                    idepth = 1
                    cubitcommand = 'refine hex in vol  ' + txt
                elif ir != nvol + 1:
                    txt = ''
                    for id_vol_ref in range(ir - 1, nvol):
                        txt = txt + str(vol[id_vol_ref].ID) + ' '
                    cubitcommand = 'refine hex in vol  ' + txt
                else:
                    # refinement on the top surface
                    _, _, _, _, _, tsurf = get_v_h_list([vol[ir - 2].ID])
                    tsurf = ' '.join(str(x) for x in tsurf)
                    idepth = 1
                    cubitcommand = 'refine node in surf ' + str(tsurf) + \
                                   ' numsplit 1 bias 1.0 depth ' + str(idepth)
                cubit.cmd(cubitcommand)

    sandwich = 'verticalsandwich_volume_ascii_regulargrid_mpiregularmap'
    if not nvol and cfg.volume_type == sandwich:
            # AAA
            # Volume 2 is always in between the 2nd and 3rd vertical
            # surfaces from the left
        cubitcommand = "refine node in volume 2 numsplit 1 depth 0"
        cubit.cmd(cubitcommand)

def mesh(filename=None):
    """create the mesh"""
    import start as start
    cfg                     = start.start_cfg(filename=filename)
    mpiflag,iproc,numproc,mpi   = start.start_mpi()
    #
    if cfg.map_meshing_type == 'regularmap':
            mesh_layercake_regularmap(filename=filename)
    else:
        print 'error: map_meshing_type ', cfg.map_meshing_type,' not implemented'

def surfaces(filename=None):
    """create the volumes"""
    import start as start
    print'volume'
    cfg                     = start.start_cfg(filename=filename)
    print cfg
    if cfg.volume_type == 'layercake_volume_ascii_regulargrid_regularmap':
            layercake_volume_ascii_regulargrid_mpiregularmap(filename=filename,onlysurface=True)
    elif cfg.volume_type == 'layercake_volume_fromacis_mpiregularmap':
            layercake_volume_fromacis_mpiregularmap(filename=filename,onlysurface=True)
    elif cfg.volume_type == 'verticalsandwich_volume_ascii_regulargrid_mpiregularmap':
            layercake_volume_ascii_regulargrid_mpiregularmap(filename=filename,verticalsandwich=True,onlysurface=True)

 def runsave(geometryfile,iproc,filename=None):
     import start as start
     cubit                   = start.start_cubit()
     cfg                         = start.start_cfg(filename=filename)
     flag=[0]
     ner=cubit.get_error_count()
     cubitcommand= 'save as "'+ cfg.output_dir+'/'+geometryfile+ '"  overwrite' 
     cubit.cmd(cubitcommand)                                                    
     ner2=cubit.get_error_count()                                             
     if ner == ner2:
         flag=[1]
     return flag

def surfaces(filename=None):
    """creating the surfaces defined in the parameter files
       #
       nsurf = number of surfaces
       surf_type = list of stype of method for the creation of the surface
                   -- regulare_grid (u and v lines)
                   -- skin (u lines)       
    """
    #
    import start as start
    cfg                     = start.start_cfg(filename=filename)
    #
    #
    for isurface in range(0,cfg.nsurf):
        surf_type=cfg.surf_type[isurface]
        if surf_type == 'regular_grid':
            surface_regular_grid(isurface,cfgname=filename)
        elif surf_type == 'skin':
            surface_skin(isurface,cfgname=filename)

def read_grid(filename=None):
    import sys
    import start as start
    mpiflag,iproc,numproc,mpi   = start.start_mpi()
    #
    numpy                       = start.start_numpy()
    cfg                         = start.start_cfg(filename=filename)
    from utilities import geo2utm
    
    #if cfg.irregulargridded_surf==True then cfg.nx and cfg.ny are the desired number of point along the axis....
    if cfg.nx and cfg.ny:
        nx=cfg.nx
        ny=cfg.ny
        if cfg.nstep:
            nx=min(cfg.nx,int(cfg.nx/cfg.nstep)+1)
            ny=min(cfg.ny,int(cfg.ny/cfg.nstep)+1)
            nstep=cfg.nstep
        else:
            nstep=1
    else:
        try:
            xstep=cfg.step
            ystep=cfg.step
        except:
            xstep=cfg.xstep
            ystep=cfg.ystep
        nx= int((cfg.longitude_max-cfg.longitude_min)/xstep)+1
        ny= int((cfg.latitude_max-cfg.latitude_min)/ystep)+1
        nstep=1
    #
    
    if cfg.irregulargridded_surf:
        xt,xstep=numpy.linspace(cfg.xmin, cfg.xmax, num=nx, retstep=True)
        yt,ystep=numpy.linspace(cfg.ymin, cfg.ymax, num=ny, retstep=True)
    
    elev=numpy.zeros([nx,ny,cfg.nz],float)
    #
    if  cfg.bottomflat: 
        elev[:,:,0] = cfg.depth_bottom
        bottomsurface=1
    else:
        bottomsurface=0
        
    for inz in range(bottomsurface,cfg.nz-1):
        grdfilename=cfg.filename[inz-bottomsurface]

        if cfg.irregulargridded_surf:
            coordx,coordy,elev_1=process_irregular_surfacefiles(iproc,nx,ny,cfg.xmin,cfg.xmax,cfg.ymin,cfg.ymax,xstep,ystep,grdfilename)
        else:
            coordx,coordy,elev_1=process_surfacefiles(iproc,nx,ny,nstep,grdfilename,cfg.unit,cfg.lat_orientation)
        elev[:,:,inz]=elev_1[:,:]
        #
    
    inz=cfg.nz-1 #last surface
    if cfg.sea:
        elev[:,:,inz]=elev[:,:,inz-1]
    else:
        #try:
        grdfile = cfg.filename[inz-bottomsurface]
        print 'reading ',cfg.filename[inz-bottomsurface]
        if cfg.irregulargridded_surf:
            coordx,coordy,elev_1=process_irregular_surfacefiles(iproc,nx,ny,cfg.xmin,cfg.xmax,cfg.ymin,cfg.ymax,xstep,ystep,grdfile)
        else:
            coordx,coordy,elev_1=process_surfacefiles(iproc,nx,ny,nstep,grdfile,cfg.unit,cfg.lat_orientation)
        elev[:,:,inz]=elev_1[:,:]
        #except:
        #     txt='error reading: '+  str( cfg.filename[inz-bottomsurface] )
        #    raise NameError, txt
        
        
        if cfg.subduction:
          print 'subduction'
          top=elev[:,:,inz]
          slab=elev[:,:,inz-1]
          subcrit=numpy.abs(top-slab)<cfg.subduction_thres
          top[subcrit]=slab[subcrit]+cfg.subduction_thres
          print len(top[subcrit])
          elev[:,:,inz]=top
    return coordx,coordy,elev,nx,ny

        cpuxmax = 1
    print cpuxmax, cpuymax

    if cpml:
        if not cpml_size:
            print 'specify the size of the cpml boundaries'
            import sys
            sys.exit()
        elif cpml_size <= 0:
            print 'cpml size negative, please check the parameters'
            import sys
            sys.exit()

    if chkcfg is True:
        import start as start
        cfg = start.start_cfg()
        d = cfg.__dict__
        ks = d.keys()
        ks.sort()
        for k in ks:
            if '__' not in k and '<' not in str(d[k]) and d[k] is not None:
                txt = str(k) + ' -----> ' + \
                    str(d[k])
                txt = txt.replace("'", "").replace('"', '')
                print txt
    else:
        try:
            import start as start
            cfg = start.start_cfg()
            f = open('cfg.log', 'w')
            print>>f, 'CFG FILE: ', cfg_name

def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandwich=False,onlysurface=False):
    import sys
    import start as start
    #
    mpiflag,iproc,numproc,mpi   = start.start_mpi()
    #
    numpy                       = start.start_numpy()
    cfg                         = start.start_cfg(filename=filename)                       
    
    from utilities import geo2utm, savegeometry,savesurf,cubit_command_check
    
    from math import sqrt
    #
    try:
        mpi.barrier()
    except:
        pass
    #
    #
    command = "comment '"+"PROC: "+str(iproc)+"/"+str(numproc)+" '"
    cubit_command_check(iproc,command,stop=True)
    if verticalsandwich: cubit.cmd("comment 'Starting Vertical Sandwich'")
    #
    #get icpuy,icpux values
    if mpiflag:
        icpux = iproc % cfg.nproc_xi
        icpuy = int(iproc / cfg.nproc_xi)
    else:
        icpuy=int(cfg.id_proc/cfg.nproc_xi)
        icpux=cfg.id_proc%cfg.nproc_xi
    
    #
    if  cfg.geometry_format == 'ascii':
        #for the original surfaces
        #number of points in the files that describe the topography
        import local_volume
        if cfg.localdir_is_globaldir:
            if iproc == 0 or not mpiflag:
                coordx_0,coordy_0,elev_0,nx_0,ny_0=local_volume.read_grid(filename)
                print 'end of reading grd files '+str(nx_0*ny_0)+ ' points'
            else:
                pass
            if iproc == 0 or not mpiflag:
                coordx=mpi.bcast(coordx_0)
            else:
                coordx=mpi.bcast()
            if iproc == 0 or not mpiflag:
                coordy=mpi.bcast(coordy_0)
            else:
                coordy=mpi.bcast()
            if iproc == 0 or not mpiflag:
                elev=mpi.bcast(elev_0)
            else:
                elev=mpi.bcast()
            if iproc == 0 or not mpiflag:
                nx=mpi.bcast(nx_0)
            else:
                nx=mpi.bcast()       
            if iproc == 0 or not mpiflag:
                ny=mpi.bcast(ny_0)
            else:
                ny=mpi.bcast()
        else:
            coordx,coordy,elev,nx,ny=local_volume.read_grid(filename)
        print str(iproc)+ ' end of receving grd files '
        nx_segment=int(nx/cfg.nproc_xi)+1
        ny_segment=int(ny/cfg.nproc_eta)+1
        
    elif cfg.geometry_format=='regmesh': # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        if cfg.depth_bottom != cfg.zdepth[0]:
            if iproc == 0: print 'the bottom of the block is at different depth than depth[0] in the configuration file'
        nx= cfg.nproc_xi+1
        ny= cfg.nproc_eta+1
        nx_segment=2
        ny_segment=2
        #if iproc == 0: print nx,ny,cfg.cpux,cfg.cpuy
        xp=(cfg.xmax-cfg.xmin)/float((nx-1))
        yp=(cfg.ymax-cfg.ymin)/float((ny-1))
        #
        elev=numpy.zeros([nx,ny,cfg.nz],float)
        coordx=numpy.zeros([nx,ny],float)
        coordy=numpy.zeros([nx,ny],float)
        #
        #
        xlength=(cfg.xmax-cfg.xmin)/float(cfg.nproc_xi) #length of x slide for chunk
        ylength=(cfg.ymax-cfg.ymin)/float(cfg.nproc_eta) #length of y slide for chunk
        nelem_chunk_x=1    
        nelem_chunk_y=1
        ivxtot=nelem_chunk_x+1
        ivytot=nelem_chunk_y+1 
        xstep=xlength #distance between vertex on x
        ystep=ylength
        for i in range(0,cfg.nz):
            elev[:,:,i] = cfg.zdepth[i]
        
        icoord=0
        for iy in range(0,ny):
            for ix in range(0,nx):
                icoord=icoord+1
#.........这里部分代码省略.........

def layercake_volume_fromacis_mpiregularmap(filename=None):
    import sys
    import start as start
    #
    mpiflag,iproc,numproc,mpi   = start.start_mpi()
    #
    cfg                         = start.start_cfg(filename=filename)                       
    #
    from utilities import geo2utm, savegeometry
    #
    from math import sqrt
    #
    try:
        mpi.barrier()
    except:
        pass
    #
    #
    command = "comment '"+"PROC: "+str(iproc)+"/"+str(numproc)+" '"
    cubit.cmd(command)
    #
    #get the limit of the volume considering the cpu
    def xwebcut(x):
        command='create planar surface with plane xplane offset '+str(x)
        cubit.cmd(command)
        last_surface=cubit.get_last_id("surface")
        command="webcut volume all tool volume in surf "+str(last_surface)
        cubit.cmd(command)
        command="del surf "+str(last_surface)
        cubit.cmd(command)
        
    def ywebcut(x):
        command='create planar surface with plane yplane offset '+str(x)
        cubit.cmd(command)
        last_surface=cubit.get_last_id("surface")
        command="webcut volume all tool volume in surf "+str(last_surface)
        cubit.cmd(command)
        command="del surf "+str(last_surface)
        cubit.cmd(command)
        
    def translate2zero():
        ss=cubit.parse_cubit_list('surface','all')
        box = cubit.get_total_bounding_box("surface", ss)
        xmin=box[0]
        ymin=box[3]
        cubit.cmd('move surface all x '+str(-1*xmin)+' y '+str(-1*ymin))
        return xmin,ymin
        
    def translate2original(xmin,ymin):
        cubit.cmd('move surface all x '+str(xmin)+' y '+str(ymin))
        
    if mpiflag:
        icpux = iproc % cfg.nproc_xi
        icpuy = int(iproc / cfg.nproc_xi)
    else:
        icpuy=int(cfg.id_proc/cfg.nproc_xi)
        icpux=cfg.id_proc%cfg.nproc_xi
    #
    ner=cubit.get_error_count()
    #
    icurve=0
    isurf=0
    ivertex=0
    #
    xlength=(cfg.xmax-cfg.xmin)/float(cfg.cpux) #length of x slide for chunk
    ylength=(cfg.ymax-cfg.ymin)/float(cfg.cpuy) #length of y slide for chunk
    xmin_cpu=cfg.xmin+(xlength*(icpux))
    ymin_cpu=cfg.ymin+(ylength*(icpuy))
    xmax_cpu=xmin_cpu+xlength
    ymax_cpu=ymin_cpu+ylength
    #
    #importing the surfaces
    for inz in range(cfg.nz-2,-2,-1):
        if cfg.bottomflat and inz==-1:
            command = "create planar surface with plane zplane offset "+str(cfg.depth_bottom)
            cubit.cmd(command)
        else:
            command = "import cubit '"+cfg.filename[inz]+"'"
            cubit.cmd(command)
            
            
    #translate
    xmin,ymin=translate2zero()
    print 'translate ...', -xmin,-ymin
    xmin_cpu=xmin_cpu-xmin
    ymin_cpu=ymin_cpu-ymin
    xmax_cpu=xmax_cpu-xmin
    ymax_cpu=ymax_cpu-ymin
    
    ss=cubit.parse_cubit_list('surface','all')
    box = cubit.get_total_bounding_box("surface", ss)
    print 'dimension... ', box
    #cutting the surfaces
    xwebcut(xmin_cpu)
    xwebcut(xmax_cpu)
    ywebcut(ymin_cpu)
    ywebcut(ymax_cpu)
    #
    list_surface_all=cubit.parse_cubit_list("surface","all")
    #condisidering only the surfaces inside the boundaries
#.........这里部分代码省略.........

def surface_regular_grid(isurface=0,cfgname=None):
    """
    create an acis surface from a regular lon/lat/z grid
    """
    import sys,os
    from math import sqrt
    from utilities import geo2utm
    import start as start
    #
    #
    cfg                     = start.start_cfg(cfgname)
    numpy                   = start.start_numpy()
    #
    def create_line_u(ind,n,step,data,unit):
        last_curve_store=cubit.get_last_id("curve")
        command='create curve spline '
        for i in range(0,n):
            if i%step == 0:
                lon,lat,z=data[i+ind][0],data[i+ind][1],data[i+ind][2]
                x,y=geo2utm(lon,lat,unit)
                txt=' Position ' +   str(x)  +' '+  str(y) +' '+  str(z)
                command=command+txt
                #print command
        cubit.silent_cmd(command)
        last_curve=cubit.get_last_id("curve")
        if last_curve != last_curve_store:
            return last_curve
        else:
            return 0
    
    def create_line_v(ind,n,n2,step,data,unit):
        last_curve_store=cubit.get_last_id("curve")
        command='create curve spline '
        for i in range(0,n):
            if i%step == 0:
                lon,lat,z=data[n2*i+ind][0],data[n2*i+ind][1],data[n2*i+ind][2]
                x,y=geo2utm(lon,lat,unit)
                txt=' Position ' +   str(x)  +' '+  str(y) +' '+  str(z)
                command=command+txt
                #print command
        cubit.silent_cmd(command)
        last_curve=cubit.get_last_id("curve")
        if last_curve != last_curve_store:
            return last_curve
        else:
            return 0
    #
    #
    cubit.cmd("reset")
    #
    position=True
    #
    #
    nu= cfg.num_x[isurface]
    nv= cfg.num_y[isurface]
    ustep= cfg.xstep[isurface]
    vstep= cfg.ystep[isurface]
    exag=1.
    unit=cfg.unit2[isurface]
    #
    #
    data=numpy.loadtxt(cfg.surface_name[isurface])
    if len(data) > 100:
        command = "set echo off"
        cubit.cmd(command)
        command = "set journal off"
        cubit.cmd(command)
    #
    u_curve=[]
    v_curve=[]
    #
    for iv in range(0,nv):
        if iv%vstep == 0.:
            u=create_line_u(iv*(nu),nu,ustep,data,unit)
            u_curve.append(u)
    for iu in range(0,nu):
        if iu%ustep == 0.:
            v=create_line_v(iu,nv,nu,ustep,data,unit)
            v_curve.append(v)
    #
    umax=max(u_curve)
    umin=min(u_curve)
    vmax=max(v_curve)
    vmin=min(v_curve)
    cubitcommand= 'create surface net u curve '+ str( umin )+' to '+str( umax )+ ' v curve '+ str( vmin )+ ' to '+str( vmax )+' heal'
    cubit.cmd(cubitcommand)
    command = "del curve all"
    cubit.cmd(command)
    suff=cfg.surface_name[isurface].split('/')
    command = "save as '"+cfg.working_dir+"/surf_"+suff[-1]+".cub' overwrite"
    cubit.cmd(command)
    #
    #
    #        
    cubit.cmd("set info "+cfg.cubit_info)
    cubit.cmd("set echo "+cfg.echo_info)
    cubit.cmd("set journal "+cfg.jou_info)

def surface_skin(isurface=0,cfgname=None):
    """
    create an acis surface interpolating no-intersecting lines
    """
    import sys,os
    from math import sqrt
    from utilities import geo2utm
    import start as start
    #
    #
    cubit                   = start.start_cubit()
    cfg                     = start.start_cfg(cfgname)
    #
    def define_next_line(directionx,directiony,n,data):
        ndata=len(data)
        command=''
        ind=n
        try:
            record=data[ind]
        except:
            return False,False
        try:
            x,y,z=map(float,record.split())
        except:
            return False,False
        txt=' Position ' +   record
        command=command+txt
        x_store,y_store,z_store = x,y,z
        icount=1
        while True:
                    ind+=1
                    if ind >= ndata: return ind,command
                    record=data[ind]
                    try:
                        x,y,z=map(float,record.split())
                    except:
                        return ind,command
                    dx,dy = x-x_store,y-y_store        
                    if  directionx == 0 and dy/abs(dy) * directiony >= 0:
                        txt=' Position ' +   record
                        command=command+txt
                        icount+=1
                        x_store,y_store,z_store = x,y,z
                    elif  directiony == 0 and dx/abs(dx) == directionx :
                        txt=' Position ' +   record
                        command=command+txt
                        icount+=1
                        x_store,y_store,z_store = x,y,z
                    else:
                        if icount==1:
                           x,y,z=x_store+1e-4*directionx,y_store+1e-4*directiony,z_store
                           txt=' Position ' +str(x)+ ' '+str(y)+ ' '+str(z)
                           command=command+txt
                        return ind,command    
    def create_line(position):
        if position:
            last_curve_store=cubit.get_last_id("curve")
            command='create curve spline '+position
            cubit.silent_cmd(command)
            last_curve=cubit.get_last_id("curve")
            if last_curve != last_curve_store:
                return last_curve
            else:
                return False
        else:
            return False
                        
    command = "reset"
    cubit.cmd(command)
    #
    position=True
    #
    try:
         grdfile = open(cfg.surface_name[isurface], 'r')
    except:
         raise NameError, 'No such file or directory: '+  str( cfg.surface_name[isurface] )
    #
    directionx=cfg.directionx[isurface]
    directiony=cfg.directiony[isurface]
    step=cfg.step[isurface]
    position=True
    curveskin=[]
    count_line=0
    data=grdfile.read().split('\n')
    ndata=len(data)
    n=0
    #
    #
    command = "set echo off"
    cubit.cmd(command)
    command = "set journal off"
    cubit.cmd(command)
    command = "set info off"
    cubit.cmd(command)         
    #
    while position:
          index,position=define_next_line(directionx,directiony,n,data)
          if n%step == 0:
              curve=create_line(position)
              if curve: curveskin.append(curve)
#.........这里部分代码省略.........

def mesh_layercake_regularmap(filename=None):
    import sys,os
    import start as start
    mpiflag,iproc,numproc,mpi   = start.start_mpi()
    from utilities import  importgeometry,savemesh,get_v_h_list,cubit_command_check
    #
    numpy                       = start.start_numpy()
    cfg                         = start.start_cfg(filename=filename)
    from math import sqrt
    from sets import Set

    #
    class cubitvolume:
          def __init__(self,ID,intervalv,centerpoint,dimension):
              self.ID=ID
              self.intervalv=intervalv
              self.centerpoint=centerpoint
              self.dim=dimension
          
          def __repr__(self):
              msg="(vol:%3i, vertical interval: %4i, centerpoint: %8.2f)" % (self.ID, self.intervalv,self.centerpoint)
              return msg       
    #
    def by_z(x,y):
        return cmp(x.centerpoint,y.centerpoint)
    #
    #
    #
    list_vol=cubit.parse_cubit_list("volume","all")
    if len(list_vol) != 0:
        pass
    else:
        geometryfile='geometry_vol_'+str(iproc)+'.cub'
        importgeometry(geometryfile,iproc=iproc)
    #
    command = 'composite create curve all'
    cubit.cmd(command)
    print '###"No valid composites can be created from the specified curves."  is NOT a critical ERROR.'
    #
    command = "compress all"
    cubit.cmd(command)
    list_vol=cubit.parse_cubit_list("volume","all")
    nvol=len(list_vol)                                 
    vol=[]
    for id_vol in list_vol:
        p=cubit.get_center_point("volume",id_vol)
        vol.append(cubitvolume(id_vol,1,p[2],0))
    vol.sort(by_z)
    #
    for id_vol in range(0,nvol):
        vol[id_vol].intervalv=cfg.iv_interval[id_vol]
    #
    #
    surf_vertical=[]
    surf_or=[]
    top_surface=0
    top_surface_add=''
    bottom_surface=0
    #
    zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6]
    xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
    xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
    ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
    ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
    #
    #
    #interval assignement
    surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top = get_v_h_list(list_vol,chktop=cfg.chktop)
    print 'vertical surfaces: ',surf_vertical    
    
    for k in surf_vertical:
        command = "surface "+str(k)+" scheme submap"
        cubit.cmd(command)
    for k in surf_or:
        command = "surface "+str(k)+" scheme "+cfg.or_mesh_scheme
        cubit.cmd(command)
    #
    ucurve,vcurve=get_uv_curve(list_curve_or)
    schemepave=False
    #
    ucurve_interval={}
    for k in ucurve:
        length=cubit.get_curve_length(k)
        interval=int(2*round(.5*length/cfg.size,0))
        ucurve_interval[k]=interval
        command = "curve "+str(k)+" interval "+str(interval)
        cubit.cmd(command)
        #cubit_error_stop(iproc,command,ner)
        command = "curve "+str(k)+" scheme equal"
        cubit.cmd(command)
        #cubit_error_stop(iproc,command,ner)
    if max(ucurve_interval.values()) != min(ucurve_interval.values()):
        schemepave=True
        print 'mesh scheme is set to pave'
        for sk in surf_or:
            command = "surface "+str(sk)+" scheme pave"
            cubit.cmd(command)
    #
    vcurve_interval={}
    for k in vcurve:
#.........这里部分代码省略.........

def refinement(nvol,vol,filename=None):
    import start as start
    cfg                         = start.start_cfg(filename=filename)
    from utilities import get_v_h_list
    #
    #vertical refinement
    #for nvol = 3 
    #
    #___________________________ interface 4
    #                 
    #vol 2              
    #___________________________ interface 3
    #
    #vol 1
    #___________________________ interface 2
    #
    #vol 0
    #___________________________ interface 1
    #
    #
    if cfg.ntripl != 0:
        if len(cfg.refinement_depth) != 0:
            #get the topo surface....
            #surf=cubit.get_relatives('volume',vol[nvol-1].ID,'surface')
            #zstore=[-1,-999999999]
            #for s in surf:
            #     c=cubit.get_center_point('surface',s)
            #     z=c[2]
            #     print s,z
            #     if z > zstore[1]:
            #         zstore=[s,z]
            #tsurf=zstore[0]
            _,_,_,_,_,tsurf = get_v_h_list([vol[nvol-1].ID])
            tsurf=' '.join(str(x) for x in tsurf)
            for idepth in cfg.refinement_depth:
                 cubitcommand= 'refine node in surf  '+str(tsurf)+' numsplit 1 bias 1.0 depth '+str(idepth)
                 cubit.cmd(cubitcommand)
        else:
            for ir in cfg.tripl:
                if ir == 1:
                   command = "comment '"+"interface = 1 means that the refinement interface is at the bottom of the volume"+"'"
                   cubit.cmd(command)
                   txt=' all '
                   idepth = 1
                   cubitcommand= 'refine hex in vol  '+txt
                elif ir != nvol+1:
                   txt=''
                   for id_vol_ref in range(ir-1,nvol):
                       txt=txt+str(vol[id_vol_ref].ID)+' '
                   #txt=txt+'except hex in vol '+str(vol[ir-2].ID)
                   #idepth = 1
                   #try:
                   #     if  cfg.refine_basin:
                   #         idepth=2
                   #except:
                   #     pass
                   cubitcommand= 'refine hex in vol  '+txt
                else:
                   #refinement on the top surface
                   _,_,_,_,_,tsurf = get_v_h_list([vol[ir-2].ID])
                   tsurf=' '.join(str(x) for x in tsurf)
                   idepth=1
                   cubitcommand= 'refine node in surf '+str(tsurf)+' numsplit 1 bias 1.0 depth '+str(idepth)
                cubit.cmd(cubitcommand)

        if not nvol and cfg.volume_type == 'verticalsandwich_volume_ascii_regulargrid_mpiregularmap':
            # AAA
            # Volume 2 is always in between the 2nd and 3rd vertical surfaces from the left
            cubitcommand = "refine node in volume 2 numsplit 1 depth 0"
            cubit.cmd(cubitcommand)

def savemesh(mpiflag, iproc=0, filename=None):
    import start as start
    cfg = start.start_cfg(filename=filename)
    mpiflag, iproc, numproc, mpi = start.start_mpi()

    def runsave(meshfile, iproc, filename=None):
        import start as start
        cubit = start.start_cubit()
        cfg = start.start_cfg(filename=filename)
        flag = 0
        ner = cubit.get_error_count()
        cubitcommand = 'save as "' + cfg.output_dir + \
            '/' + meshfile + '.cub' + '" overwrite'
        cubit.cmd(cubitcommand)
        ner2 = cubit.get_error_count()
        if ner == ner2:
            cubitcommand = 'export mesh "' + cfg.output_dir + '/' + \
                meshfile + '.e' + '" dimension 3 block all overwrite'
            cubit.cmd(cubitcommand)
            ner2 = cubit.get_error_count()
        if ner == ner2:
            flag = 1
        return flag

    meshfile = 'mesh_vol_' + str(iproc)

    flagsaved = 0
    infosave = (iproc, flagsaved)

    mpi.barrier()
    total_saved = mpi.allgather(flagsaved)
    if isinstance(total_saved, int):
        total_saved = [total_saved]

    ind = 0
    saving = True
    while saving:
        if len(total_saved) != sum(total_saved):
            #
            if not flagsaved:
                flagsaved = runsave(meshfile, iproc, filename=filename)
                if flagsaved:
                    infosave = (iproc, flagsaved)
                    if numproc > 1:
                        f = open('mesh_saved' + str(iproc), 'w')
                        f.close()
            mpi.barrier()
            total_saved = mpi.allgather(flagsaved)
            if isinstance(total_saved, int):
                total_saved = [total_saved]
            ind = ind + 1
        else:
            saving = False
        if ind > len(total_saved) + 10:
            saving = False
        print sum(total_saved), '/', len(total_saved), ' saved'

    info_total_saved = mpi.allgather(infosave)
    if isinstance(info_total_saved, int):
        info_total_saved = [info_total_saved]

    if iproc == 0:
        f = open('mesh_saving.log', 'w')
        f.write('\n'.join(str(x) for x in info_total_saved))
        f.close()

    f = open(cfg.output_dir + '/' + 'blocks_' + str(iproc).zfill(5), 'w')
    blocks = cubit.get_block_id_list()

    for block in blocks:
        name = cubit.get_exodus_entity_name('block', block)
        element_count = cubit.get_exodus_element_count(block, "block")
        nattrib = cubit.get_block_attribute_count(block)
        attr = [cubit.get_block_attribute_value(
            block, x) for x in range(0, nattrib)]
        ty = cubit.get_block_element_type(block)
        f.write(str(block) + ' ; ' + name + ' ; nattr ' + str(nattrib) +
                ' ; ' + ' '.join(str(x) for x in attr) + ' ; ' + ty + ' ' +
                str(element_count) + '\n')
    f.close()

    import quality_log
    f = open(cfg.output_dir + '/' + 'quality_' + str(iproc).zfill(5), 'w')
    max_skewness, min_length = quality_log.quality_log(f)
    f.close()

    count_hex = [cubit.get_hex_count()]
    count_node = [cubit.get_node_count()]
    max_skew = [(iproc, max_skewness)]
    min_l = [(iproc, min_length)]

    mpi.barrier()
    total_min_l = mpi.gather(min_l)
    total_hex = mpi.gather(count_hex)
    total_node = mpi.gather(count_node)
    total_max_skew = mpi.gather(max_skew)

    mpi.barrier()
    if iproc == 0:
        min_total_min_l = min([ms[1] for ms in total_min_l])
#.........这里部分代码省略.........