Commit 1c51da8e authored by Antoine MARTIN's avatar Antoine MARTIN

modifs fractions intermediaires

parent 373457c0
......@@ -65,7 +65,7 @@ def verif_im(c,i):
def verif_eps(c,i,j,niter):
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
EpsHH=fi['EpsHH'][j,niter-1]
frac=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8]
frac=[0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.8]
f=frac[j]
plt.figure()
pp=plt.imshow(EpsHH[:,:,0,0],origin='lower')
......
......@@ -38,7 +38,9 @@ for c in [100,0.01]:
print("échantillon "+str(i)+'\n')
fname = "c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5"
fi = h5py.File(fname, "w")
im,f0=image_init(N,D,0.1)
im,centers_diam,f0=image_init(N,D,0.8)
l=len(centers_diam)
fi['centers']=np.zeros((8,l,3))
fi['images']=np.zeros((8,N,N))
fi['fractions']=np.zeros((8,))
fi['ChomEF']=np.zeros((8,3,3))
......@@ -46,19 +48,20 @@ for c in [100,0.01]:
fi['ChomHH']=np.zeros((8,10,2,3,3))
fi['EpsHH']=np.zeros((8,10,N,N,3,3))
for j,f in enumerate([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8]):
for j,f in enumerate([0.8,0.6,0.5,0.4,0.3,0.2,0.1,0.05]):
os.system("mkdir c="+str(c)+"/ech="+str(i)+"/frac="+str(f))
print("fraction volumique "+str(f)+'\n')
if f!=0.1:
f1=image_augmente(im,D,f)
if f!=0.8:
f1=image_diminue(im,centers_diam,f)
else:
f1=f0
# plt.figure()
# plt.imshow(im[:,:],origin='lower')
# plt.show()
# plt.close()
fi['images'][j]=im
fi['fractions'][j]=f1
fi['images'][7-j]=im
fi['centers'][7-j]=np.array(centers_diam)
fi['fractions'][7-j]=f1
tC=np.array([[C[0] for ii in range(N)] for jj in range(N)])
for ii in range(N):
......@@ -67,22 +70,22 @@ for c in [100,0.01]:
tC[ii,jj]=C[1]
ChomEF,CVoigt = reference(tC,256,c,i,f)
fi['ChomEF'][j]=ChomEF
fi['CVoigt'][j]=CVoigt
fi['ChomEF'][7-j]=ChomEF
fi['CVoigt'][7-j]=CVoigt
ChomHH_list,tA_list=homog_hierar(tC,CVoigt)
l=len(ChomHH_list)
if l<10:
for k in range(l):
fi['ChomHH'][j,k]=ChomHH_list[k]
fi['EpsHH'][j,k]=tA_list[k]
fi['ChomHH'][7-j,k]=ChomHH_list[k]
fi['EpsHH'][7-j,k]=tA_list[k]
for k in range(l,10):
fi['ChomHH'][j,k]=ChomHH_list[l-1]
fi['EpsHH'][j,k]=tA_list[l-1]
fi['ChomHH'][7-j,k]=ChomHH_list[l-1]
fi['EpsHH'][7-j,k]=tA_list[l-1]
else:
for k in range(10):
fi['ChomHH'][j,k]=ChomHH_list[k]
fi['EpsHH'][j,k]=tA_list[k]
fi['ChomHH'][7-j,k]=ChomHH_list[k]
fi['EpsHH'][7-j,k]=tA_list[k]
......
......@@ -15,17 +15,16 @@ def compteur(im):
return cpt/nn/mm
def ajoute(im,D,nb):
def ajoute(im,centers_diam,D,nb):
nn,mm=im.shape
centers=[[0.,0.] for i in range(nb)]
ll=len(centers_diam)
for i in range(nb):
centers[i]=[uniform(0,1),uniform(0,1)]
centers_diam.append([uniform(0,1),uniform(0,1),D])
for i in range(nb):
r=D/2
x=centers[i][0]
y=centers[i][1]
x=centers_diam[ll+i][0]
y=centers_diam[ll+i][1]
for j in range(int((x-r)*nn)-1,int((x+r)*nn)+1):
for k in range(int((y-r)*mm)-1,int((y+r)*mm)+1):
x1=(j+0.5)/nn
......@@ -41,6 +40,35 @@ def ajoute(im,D,nb):
m=0
im[j-l,k-m]=1
def enleve(im,centers_diam,nb):
nn,mm=im.shape
im2=np.zeros((nn,mm))
for i in range(nb):
centers_diam.pop()
ll=len(centers_diam)
for i in range(ll):
r=centers_diam[i][2]/2
x=centers_diam[i][0]
y=centers_diam[i][1]
for j in range(int((x-r)*nn)-1,int((x+r)*nn)+1):
for k in range(int((y-r)*mm)-1,int((y+r)*mm)+1):
x1=(j+0.5)/nn
y1=(k+0.5)/mm
if (x1-x)**2+(y1-y)**2<=r**2:
if j>=nn:
l=nn
else:
l=0
if k>=mm:
m=mm
else:
m=0
im2[j-l,k-m]=1
return im2
def image_init(N,D,f):
......@@ -49,16 +77,17 @@ def image_init(N,D,f):
n=int(f/np.pi/(D/2)**2)
centers_diam=[]
c=0
while c<f:
ajoute(im,D,n)
ajoute(im,centers_diam,D,n)
c=compteur(im)
n=int((f-c)/np.pi/(D/2)**2)+1
return im,c
return im,centers_diam,c
def image_augmente(im,D,f):
def image_augmente(im,centers_diam,D,f):
f0=compteur(im)
if f<=f0:
print('fraction volumique actuelle :'+str(f0))
......@@ -66,11 +95,34 @@ def image_augmente(im,D,f):
c=f0
n=int((f-c)/np.pi/(D/2)**2)+1
while c<f:
ajoute(im,D,n)
ajoute(im,centers_diam,D,n)
c=compteur(im)
n=int((f-c)/np.pi/(D/2)**2)+1
return c
def image_diminue(im,centers_diam,f):
f0=compteur(im)
if f>f0:
print('fraction volumique actuelle :'+str(f0))
else:
c=f0
im2=im
while c>f:
l=len(centers_diam)
cc=0.
i=0
while cc < c-f:
r=centers_diam[l-i-1][2]/2
cc+=np.pi*r**2
i+=1
im2=enleve(im2,centers_diam,i)
c=compteur(im2)
nn,mm=im.shape
for i in range(nn):
for j in range(mm):
im[i,j]=im2[i,j]
return c
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment