Commit e8146c2a authored by Antoine MARTIN's avatar Antoine MARTIN

add explo2 explohh

parent 34ee12d2
...@@ -8,16 +8,18 @@ from math import * ...@@ -8,16 +8,18 @@ from math import *
from random import uniform from random import uniform
def homog_hierar(tabC,CVoigt): def homog_hierar(tabC):
# K0=np.array([[0.5,-0.5, 0. ],
# [-0.5, 0.5, 0. ],
# [0. , 0. , 1.]])
# J0=np.array([[0.5, 0.5, 0. ],
# [0.5, 0.5, 0. ],
# [0. , 0. , 0.]])
pp,qq,rr,ss=tabC.shape pp,qq,rr,ss=tabC.shape
CVoigt = np.array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
for i in range(pp):
for j in range(pp):
CVoigt=CVoigt+tabC[i,j]
CVoigt=CVoigt/pp/pp
p=int(log(pp)/log(2)) p=int(log(pp)/log(2))
...@@ -37,8 +39,7 @@ def homog_hierar(tabC,CVoigt): ...@@ -37,8 +39,7 @@ def homog_hierar(tabC,CVoigt):
#la fonction auxiliary_problem calcule les A locaux : l'argument principal est C_list, la liste des tenseurs d'élasticité, et la fonction renvoie Aloc, un tableau à trois dimensions. Aloc[i] correspond à la moyenne du tenseur de localisation sur la case i, calculé par homogénéisation périodique def auxiliary_problem(C_list, n,gamma,K):
def auxiliary_problem(C_list, n,gamma,K):
M=np.array([[0. for i in range(12)] for j in range(12)]) M=np.array([[0. for i in range(12)] for j in range(12)])
...@@ -79,7 +80,6 @@ def homog_hierar(tabC,CVoigt): ...@@ -79,7 +80,6 @@ def homog_hierar(tabC,CVoigt):
return Chom,E return Chom,E
#la fonction tC prend en argument un tableau de type tCA et renvoie le tableau de tenseurs d'élasticité qui sera celui de l'étape d'après
def tC(tC_,tA,n,gamma,K): def tC(tC_,tA,n,gamma,K):
t1=np.array([[np.zeros((3,3)) for i in range(2**(n-1))] for j in range(2**(n-1))]) t1=np.array([[np.zeros((3,3)) for i in range(2**(n-1))] for j in range(2**(n-1))])
for i in range(2**(n-1)): for i in range(2**(n-1)):
...@@ -101,10 +101,7 @@ def homog_hierar(tabC,CVoigt): ...@@ -101,10 +101,7 @@ def homog_hierar(tabC,CVoigt):
return t1 return t1
#la fonction Chom prend en argument le tableau initial tC0 des tenseurs d'élasticité associés à la grille image de taille n, et appelle les fonctions décrites plus haut pour calculer les nouveaux tableaux de tenseurs d'élasticité associés à chaque étape. La fonction enregistre ces tableaux à chaque étape dans un fichier. Le nombre d'étapes 'e' doit être précisé pour des méthodes 10 ou 11.
def Chom(tC0,n,K_): def Chom(tC0,n,K_):
# K[2,2]=K[0,0]-K[0,1] # K[2,2]=K[0,0]-K[0,1]
......
...@@ -22,31 +22,11 @@ def local_project(v,V): ...@@ -22,31 +22,11 @@ def local_project(v,V):
return u return u
#cf 'Contenu des fichiers' à propos de 'reference' #cf 'Contenu des fichiers' à propos de 'reference'
def reference(tabC,res1,c,ec,fr): def reference(tabC,res1,c,ec,fr,tA):
#res=256 #res=256
k,ll,mm,nn=tabC.shape k,ll,mm,nn=tabC.shape
mesh = UnitSquareMesh(res1, res1, "left/right") mesh = UnitSquareMesh(res1, res1, "left/right")
def pixel_0(x,on_boundary):
return x[0]<=0.5 and x[1]<=0.5
pixels=np.array([[pixel_0 for i in range(res1)] for j in range(res1)])
for i in range(res1):
for j in range(res1):
def pixel(x,on_boundary):
return x[0]<=(j+1)/res1 and x[0]>=j/res1 and x[1]<=(i+1)/res1 and x[1]>=i/res1
pixels[i,j]=pixel()
domains = MeshFunction("size_t", mesh,2)
domains.set_all(0)
for i in range(res1):
for j in range(res1):
AutoSubDomain(pixels[i,j]).mark(domains, i*res1+j)
dx = Measure("dx", domain=mesh, subdomain_data=domains)
class LocalC(UserExpression): class LocalC(UserExpression):
def __init__(self, Cij, mesh): def __init__(self, Cij, mesh):
self.Cij = Cij self.Cij = Cij
...@@ -80,6 +60,7 @@ def reference(tabC,res1,c,ec,fr): ...@@ -80,6 +60,7 @@ def reference(tabC,res1,c,ec,fr):
#Cvar est la transformation du tableau Cv en une fonction définie sur 'mesh' #Cvar est la transformation du tableau Cv en une fonction définie sur 'mesh'
Cvar= LocalC(tabC,mesh) Cvar= LocalC(tabC,mesh)
Avar= LocalC(tA,mesh)
# class used to define the periodic boundary map # class used to define the periodic boundary map
class PeriodicBoundary(SubDomain): class PeriodicBoundary(SubDomain):
...@@ -116,7 +97,7 @@ def reference(tabC,res1,c,ec,fr): ...@@ -116,7 +97,7 @@ def reference(tabC,res1,c,ec,fr):
#la fonction reference_problem prend en argument C, donnant le tenseur d'élasticité en un point de 'mesh', et renvoie le Chom du domaine. k étant le nombre de phases différentes sur le domaine, qui correspond à la taille du tableau Cn pris en argument par la fonction 'reference' #la fonction reference_problem prend en argument C, donnant le tenseur d'élasticité en un point de 'mesh', et renvoie le Chom du domaine. k étant le nombre de phases différentes sur le domaine, qui correspond à la taille du tableau Cn pris en argument par la fonction 'reference'
def reference_problem(C, mesh,k): def reference_problem(C, mesh,k,A):
def eps(v): def eps(v):
return as_vector([v[0].dx(0), return as_vector([v[0].dx(0),
...@@ -126,6 +107,9 @@ def reference(tabC,res1,c,ec,fr): ...@@ -126,6 +107,9 @@ def reference(tabC,res1,c,ec,fr):
def sigma(v, Eps): def sigma(v, Eps):
return dot(C, eps(v)+Eps) return dot(C, eps(v)+Eps)
def sigma_HH(Eps):
return dot(C,dot(A,Eps))
Ve = VectorElement("CG", mesh.ufl_cell(), 2) Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0) Re = VectorElement("R", mesh.ufl_cell(), 0)
...@@ -143,9 +127,7 @@ def reference(tabC,res1,c,ec,fr): ...@@ -143,9 +127,7 @@ def reference(tabC,res1,c,ec,fr):
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
Chom = np.zeros((3, 3)) Chom = np.zeros((3, 3))
tA=np.array([[np.array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]]) for i in range(res1)] for j in range(res1)])
for (j, case) in enumerate(["Exx", "Eyy", "Exy"]): for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
ffile = XDMFFile('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Strain_EF_{}.xdmf'.format(case)) ffile = XDMFFile('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Strain_EF_{}.xdmf'.format(case))
...@@ -158,39 +140,36 @@ def reference(tabC,res1,c,ec,fr): ...@@ -158,39 +140,36 @@ def reference(tabC,res1,c,ec,fr):
vf.assign(local_project(eps(v), V)) vf.assign(local_project(eps(v), V))
ffile.write_checkpoint(vf, case, 0.) ffile.write_checkpoint(vf, case, 0.)
# plt.figure() diff_sig=sigma(v,Eps)-sigma_HH(Eps)
# pp=plot(vf[0]) plt.figure()
# plt.colorbar(pp) pp=plot(diff_sig[0])
# plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 1'+str(j+1)+'.pdf') plt.colorbar(pp)
# #plt.show() plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig 11'+case+'.pdf')
# plt.close() #plt.show()
# plt.figure() plt.close()
# pp=plot(vf[1]) plt.figure()
# plt.colorbar(pp) pp=plot(diff_sig[1])
# plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 2'+str(j+1)+'.pdf') plt.colorbar(pp)
# #plt.show() plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig 22'+case+'.pdf')
# plt.close() #plt.show()
# plt.figure() plt.close()
# pp=plot(vf[2]) plt.figure()
# plt.colorbar(pp) pp=plot(diff_sig[2])
# plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf') plt.colorbar(pp)
# #plt.show() plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig 12'+case+'.pdf')
# plt.close() #plt.show()
plt.close()
for ii in range(res1):
for jj in range(res1):
for kk in range(3):
tA[res1-ii,jj][kk,j]=assemble(eps(v)[kk]*dx(ii*res1+jj))+Eps[kk]
Sigma = np.zeros((3,)) Sigma = np.zeros((3,))
for kk in range(3): for kk in range(3):
Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx) Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx)
Chom[j, :] = Sigma Chom[j, :] = Sigma
ffile.close() ffile.close()
return Chom,tA return Chom
Chom,tA = reference_problem(Cvar, mesh,k) Chom = reference_problem(Cvar, mesh,k,Avar)
CVoigt = np.array([[0., 0., 0.], CVoigt = np.array([[0., 0., 0.],
[0., 0., 0.], [0., 0., 0.],
[0., 0., 0.]]) [0., 0., 0.]])
...@@ -201,7 +180,7 @@ def reference(tabC,res1,c,ec,fr): ...@@ -201,7 +180,7 @@ def reference(tabC,res1,c,ec,fr):
CVoigt=CVoigt/k/k CVoigt=CVoigt/k/k
print(Chom) print(Chom)
return Chom,CVoigt,tA return Chom,CVoigt
......
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import numpy as np import numpy as np
import h5py import h5py
...@@ -40,7 +40,7 @@ for c in [100,0.01]: ...@@ -40,7 +40,7 @@ for c in [100,0.01]:
fi = h5py.File(fname, "r") fi = h5py.File(fname, "r")
#im,centers_diam,f0=image_init(N,D,0.9) #im,centers_diam,f0=image_init(N,D,0.9)
#l=len(centers_diam) #l=len(centers_diam)
centers_diam=fi['centers'][7] centers_diam=(fi['centers'][7]).tolist()
l0=len(centers_diam) l0=len(centers_diam)
im=fi['images'][7] im=fi['images'][7]
f0=image_augmente(im,centers_diam,D,0.9) f0=image_augmente(im,centers_diam,D,0.9)
...@@ -52,7 +52,7 @@ for c in [100,0.01]: ...@@ -52,7 +52,7 @@ for c in [100,0.01]:
fractions_2=np.zeros((10,)) fractions_2=np.zeros((10,))
ChomEF_2=np.zeros((10,3,3)) ChomEF_2=np.zeros((10,3,3))
CVoigt_2=np.zeros((10,3,3)) CVoigt_2=np.zeros((10,3,3))
ChomHH_2=np.zeros((10,10,N,N,3,3) ChomHH_2=np.zeros((10,10,2,3,3))
EpsHH_2=np.zeros((10,10,N,N,3,3)) EpsHH_2=np.zeros((10,10,N,N,3,3))
for j1,j2 in enumerate([0,1,2,3,4,5,6,-1,7,-1]): for j1,j2 in enumerate([0,1,2,3,4,5,6,-1,7,-1]):
...@@ -66,7 +66,8 @@ for c in [100,0.01]: ...@@ -66,7 +66,8 @@ for c in [100,0.01]:
EpsHH_2[j1]=fi['EpsHH'][j2] EpsHH_2[j1]=fi['EpsHH'][j2]
for j,f in enumerate([0.9,0.,0.67,0.,0.,0.,0.,0.,0.,0.]): for j,f in enumerate([0.9,0.,0.67,0.,0.,0.,0.,0.,0.,0.]):
os.system("mkdir c="+str(c)+"/ech="+str(i)+"/frac="+str(f)) if f!=0:
os.system("mkdir c="+str(c)+"/ech="+str(i)+"/frac="+str(f))
print("fraction volumique "+str(f)+'\n') print("fraction volumique "+str(f)+'\n')
if f!=0.9 and f!=0: if f!=0.9 and f!=0:
......
...@@ -66,7 +66,7 @@ def verif_im(c,i): ...@@ -66,7 +66,7 @@ def verif_im(c,i):
def verif_eps(c,i,j,niter): def verif_eps(c,i,j,niter):
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
EpsHH=fi['EpsHH'][j,niter-1] EpsHH=fi['EpsHH'][j,niter-1]
frac=[0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9] frac=[0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,0.95]
f=frac[j] f=frac[j]
plt.figure() plt.figure()
pp=plt.imshow(EpsHH[:,:,0,0],origin='lower') pp=plt.imshow(EpsHH[:,:,0,0],origin='lower')
...@@ -103,27 +103,27 @@ def courbes(c,i): ...@@ -103,27 +103,27 @@ def courbes(c,i):
fr=fi['fractions'] fr=fi['fractions']
ChomEF=fi['ChomEF'] ChomEF=fi['ChomEF']
ChomHH=fi['ChomHH'] ChomHH=fi['ChomHH']
C_11=np.ones((12,)) C_11=np.ones((13,))
C_11_EF=np.ones((12,)) C_11_EF=np.ones((13,))
C_22=np.ones((12,)) C_22=np.ones((13,))
C_22_EF=np.ones((12,)) C_22_EF=np.ones((13,))
C_12=np.ones((12,)) C_12=np.ones((13,))
C_12_EF=np.ones((12,)) C_12_EF=np.ones((13,))
C_33=np.ones((12,)) C_33=np.ones((13,))
C_33_EF=np.ones((12,)) C_33_EF=np.ones((13,))
C_11_0=np.ones((12,)) C_11_0=np.ones((13,))
C_22_0=np.ones((12,)) C_22_0=np.ones((13,))
C_12_0=np.ones((12,)) C_12_0=np.ones((13,))
C_33_0=np.ones((12,)) C_33_0=np.ones((13,))
fractions=np.ones((12,)) fractions=np.ones((13,))
for j,f in enumerate([0.,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,1.]): for j,f in enumerate([0.,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,0.95,1.]):
if c==0.01: if c==0.01:
jj=11-j jj=12-j
else: else:
jj=j jj=j
if f==0. or f==1.: if f==0. or f==1.:
fractions[jj]=f fractions[jj]=f
if jj==11: if jj==12:
C_11[jj]=0.01 C_11[jj]=0.01
C_11_0[jj]=0.01 C_11_0[jj]=0.01
C_11_EF[jj]=0.01 C_11_EF[jj]=0.01
...@@ -209,29 +209,29 @@ def courbe_moyenne(c): ...@@ -209,29 +209,29 @@ def courbe_moyenne(c):
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
fichiers.append(fi) fichiers.append(fi)
C_11=np.ones((12,)) C_11=np.ones((13,))
C_11_EF=np.ones((12,)) C_11_EF=np.ones((13,))
C_22=np.ones((12,)) C_22=np.ones((13,))
C_22_EF=np.ones((12,)) C_22_EF=np.ones((13,))
C_12=np.ones((12,)) C_12=np.ones((13,))
C_12_EF=np.ones((12,)) C_12_EF=np.ones((13,))
C_33=np.ones((12,)) C_33=np.ones((13,))
C_33_EF=np.ones((12,)) C_33_EF=np.ones((13,))
C_11_0=np.ones((12,)) C_11_0=np.ones((13,))
C_22_0=np.ones((12,)) C_22_0=np.ones((13,))
C_12_0=np.ones((12,)) C_12_0=np.ones((13,))
C_33_0=np.ones((12,)) C_33_0=np.ones((13,))
fractions=np.ones((12,)) fractions=np.ones((13,))
for i in range(10): for i in range(10):
ChomHH=fichiers[i]['ChomHH'] ChomHH=fichiers[i]['ChomHH']
ChomEF=fichiers[i]['ChomEF'] ChomEF=fichiers[i]['ChomEF']
for j,f in enumerate([0.,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,1.]): for j,f in enumerate([0.,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,0.95,1.]):
if c==0.01: if c==0.01:
jj=11-j jj=12-j
else: else:
jj=j jj=j
fractions[jj]=f fractions[jj]=f
if jj==11: if jj==12:
C_11[jj]=0.01 C_11[jj]=0.01
C_11_0[jj]=0.01 C_11_0[jj]=0.01
C_11_EF[jj]=0.01 C_11_EF[jj]=0.01
...@@ -311,71 +311,71 @@ def courbe_moyenne(c): ...@@ -311,71 +311,71 @@ def courbe_moyenne(c):
plt.close() plt.close()
def erreur_abs_def(c,i,j): # def erreur_abs_def(c,i,j):
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") # fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
tAHH=fi['EpsHH'][j,9] # tAHH=fi['EpsHH'][j,9]
tAEF=fi['EpsEF'][j] # tAEF=fi['EpsEF'][j]
frac=[0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9] # frac=[0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9]
f=frac[j] # f=frac[j]
plt.figure() # plt.figure()
pp=plt.imshow(tAHH[:,:,0,0]-tAEF[:,:,0,0],origin='lower') # pp=plt.imshow(tAHH[:,:,0,0]-tAEF[:,:,0,0],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 11.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 11.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plt.imshow(tAHH[:,:,0,1]-tAEF[:,:,0,1],origin='lower') # pp=plt.imshow(tAHH[:,:,0,1]-tAEF[:,:,0,1],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 12.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 12.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plt.imshow(tAHH[:,:,1,1]-tAEF[:,:,1,1],origin='lower') # pp=plt.imshow(tAHH[:,:,1,1]-tAEF[:,:,1,1],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 22.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 22.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plt.imshow(tAHH[:,:,2,2]-tAEF[:,:,2,2],origin='lower') # pp=plt.imshow(tAHH[:,:,2,2]-tAEF[:,:,2,2],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 33.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur abs eps 33.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
def erreur_rel_def(c,i,j): # def erreur_rel_def(c,i,j):
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") # fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
tAHH=fi['EpsHH'][j,9] # tAHH=fi['EpsHH'][j,9]
tAEF=fi['EpsEF'][j] # tAEF=fi['EpsEF'][j]
ii,jj,kk,ll=tAHH.shape # ii,jj,kk,ll=tAHH.shape
err=np.zeros((ii,jj,kk,ll)) # err=np.zeros((ii,jj,kk,ll))
for k in range(ii): # for k in range(ii):
for l in range(jj): # for l in range(jj):
for m in range(kk): # for m in range(kk):
for n in range(ll): # for n in range(ll):
err[k,l,m,n]=tAHH[k,l,m,n]/tAEF[k,l,m,n]-1 # err[k,l,m,n]=tAHH[k,l,m,n]/tAEF[k,l,m,n]-1
frac=[0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9] # frac=[0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9]
f=frac[j] # f=frac[j]
plt.figure() # plt.figure()
pp=plt.imshow(err[:,:,0,0],origin='lower') # pp=plt.imshow(err[:,:,0,0],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 11.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 11.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plt.imshow(err[:,:,0,1],origin='lower') # pp=plt.imshow(err[:,:,0,1],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 12.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 12.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plt.imshow(err[:,:,1,1],origin='lower') # pp=plt.imshow(err[:,:,1,1],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 22.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 22.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plt.imshow(err[:,:,2,2],origin='lower') # pp=plt.imshow(err[:,:,2,2],origin='lower')
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 33.pdf') # plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/erreur rel eps 33.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
\ No newline at end of file \ No newline at end of file
...@@ -40,31 +40,31 @@ for c in [100,0.01]: ...@@ -40,31 +40,31 @@ for c in [100,0.01]:
print("échantillon "+str(i)+'\n') print("échantillon "+str(i)+'\n')
fname = "c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5" fname = "c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5"
fi = h5py.File(fname, "w") fi = h5py.File(fname, "w")
im,centers_diam,f1=image_init(N,D,0.9) im,centers_diam,f1=image_init(N,D,0.95)
l=len(centers_diam) l=len(centers_diam)
fi['centers']=np.zeros((10,l,3)) fi['centers']=np.zeros((11,l,3))
fi['images']=np.zeros((10,N,N)) fi['images']=np.zeros((11,N,N))
fi['fractions']=np.zeros((10,)) fi['fractions']=np.zeros((11,))
fi['ChomEF']=np.zeros((10,3,3)) fi['ChomEF']=np.zeros((11,3,3))
fi['CVoigt']=np.zeros((10,3,3)) fi['CVoigt']=np.zeros((11,3,3))
fi['ChomHH']=np.zeros((10,10,2,3,3)) fi['ChomHH']=np.zeros((11,10,2,3,3))
fi['EpsHH']=np.zeros((10,10,N,N,3,3)) fi['EpsHH']=np.zeros((11,10,N,N,3,3))
fi['EpsEF']=np.zeros((10,N,N,3,3))
for j,f in enumerate([0.9,0.8,0.67,0.6,0.5,0.4,0.3,0.2,0.1,0.05]): for j,f in enumerate([0.95,0.9,0.8,0.67,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)) os.system("mkdir c="+str(c)+"/ech="+str(i)+"/frac="+str(f))
print("fraction volumique "+str(f)+'\n') print("fraction volumique "+str(f)+'\n')
if f!=0.9: if f!=0.95:
f1=image_diminue(im,centers_diam,f) f1=image_diminue(im,centers_diam,f)
# plt.figure() # plt.figure()
# plt.imshow(im[:,:],origin='lower') # plt.imshow(im[:,:],origin='lower')
# plt.show() # plt.show()
# plt.close() # plt.close()
fi['images'][9-j]=im fi['images'][10-j]=im
l=len(centers_diam) l=len(centers_diam)
for ii in range(l): for ii in range(l):