You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
131 lines
2.9 KiB
Python
131 lines
2.9 KiB
Python
import numpy as np
|
|
import os
|
|
import time
|
|
from JoinCmaps import *
|
|
#k[x,y,z]
|
|
|
|
|
|
def div_veccon(kc,kh,nbl,rundir):
|
|
|
|
t0=time.time()
|
|
kc=np.where(kc==kh,1,0).astype(int)
|
|
|
|
tcmaps=time.time()
|
|
kc=get_smallCmap(kc,nbl,rundir)
|
|
tcmaps=time.time()-tcmaps
|
|
#if s_scale<kc.shape[0]:
|
|
kc=join(kc,nbl)
|
|
|
|
|
|
y = np.bincount(kc.reshape(-1))
|
|
ii = np.nonzero(y)[0]
|
|
cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia
|
|
if cf[0,0]==0:
|
|
cf=cf[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero
|
|
nclus=cf.shape[0] #cantidad de clusters
|
|
nper=np.sum(cf[:,1]) #num de celdas permeables
|
|
|
|
print(nbl,nclus,float(nper)/(kc.size), time.time()-t0)
|
|
|
|
return np.array([nbl,nclus,float(nper)/(kc.size),time.time()-t0, tcmaps,tcmaps/(time.time()-t0)])
|
|
|
|
|
|
|
|
def get_smallCmap(vec,nbl,rundir):
|
|
|
|
|
|
Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2]
|
|
sx,sy,sz = Nx//nbl,Ny//nbl,Nz//nbl
|
|
params, execCon = ConConfig(sx,sy,sz,Nz,rundir)
|
|
if Nz==1:
|
|
nblz=1
|
|
sz=1
|
|
else:
|
|
nblz=nbl
|
|
for i in range(nbl):
|
|
for j in range(nbl):
|
|
for k in range(nblz):
|
|
vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz]=connec(vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz],execCon,params,rundir)
|
|
return vec
|
|
|
|
|
|
def connec(vec,execCon,params,rundir):
|
|
np.savetxt(rundir+params[2],vec.reshape(-1))
|
|
os.system(rundir+execCon +'>/dev/null') #'cd ' +exeDir++'>/dev/null'
|
|
vec=np.loadtxt(params[-2]).reshape(vec.shape[0],vec.shape[1],vec.shape[2]).astype(int)
|
|
return vec
|
|
|
|
def ConConfig(sx,sy,sz,Nz,rundir):
|
|
|
|
params=[]
|
|
if Nz==1:
|
|
params=['1','4','vecconec.txt',str(sx)+' '+str(sy),'1.0 1.0','pardol.STA','pardol.CCO','pardol.COF']
|
|
execCon='conec2d'
|
|
|
|
else:
|
|
params=['1','6','vecconec.txt',str(sx)+' '+str(sy)+' ' +str(sz),'1.0 1.0 1.0','30','pardol.STA','pardol.CCO','pardol.COF']
|
|
execCon='conec3d'
|
|
|
|
|
|
with open(rundir+'coninput.txt', 'w') as f:
|
|
for item in params:
|
|
f.write("%s\n" % item)
|
|
|
|
return params, execCon
|
|
|
|
def join(vec,nbl):
|
|
|
|
|
|
Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2]
|
|
sx,sy,sz = Nx//nbl,Ny//nbl,Nz//nbl
|
|
ex,ey,ez=np.log2(Nx),np.log2(Ny),np.log2(Nz)
|
|
|
|
if Nz==1:
|
|
sz=1
|
|
nbz=1
|
|
ez=1
|
|
esz=1
|
|
else:
|
|
esz=np.log2(sz)
|
|
|
|
|
|
esx,esy=np.log2(sx),np.log2(sy)
|
|
|
|
|
|
|
|
for bs in range(0,int(ex-esx)):
|
|
|
|
nbx,nby = int(2**(ex-esx-bs-1)),int(2**(ey-esy-bs-1))
|
|
if Nz==1:
|
|
sz=1
|
|
nbz=1
|
|
else:
|
|
nbz=int(2**(ez-esz-bs-1))
|
|
sz=Nz//nbz
|
|
sx,sy=Nx//nbx,Ny//nby
|
|
|
|
for i in range(nbx):
|
|
for j in range(nby):
|
|
for k in range(nbz):
|
|
a=2
|
|
vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz]=joinBox(vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz],True,False)
|
|
|
|
return vec
|
|
|
|
'''
|
|
job=0
|
|
k=np.load('../../data/'+str(job)+'/k.npy')
|
|
div_veccon(k,100,1,'./')
|
|
div_veccon(k,100,2,'./')
|
|
div_veccon(k,100,4,'./')
|
|
'''
|
|
|
|
for job in range(6):
|
|
k=np.load('../../data/'+str(job)+'/k.npy')
|
|
print(job)
|
|
res=div_veccon(k,100,4,'./')
|
|
np.savetxt('../../data/'+str(job)+'/Cmap_res.txt',res)
|
|
res=div_veccon(k,100,1,'./')
|
|
#div_veccon(k,100,64,'./')
|
|
#div_veccon(k,100,128,'./')
|