Commit 2897e76a authored by Jeremy BLEYER's avatar Jeremy BLEYER

Add Antoine files

parents
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from Reference import reference
from math import *
from random import uniform
def homog_hierar(tabC,CVoigt):
# 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
p=int(log(pp)/log(2))
gamma0=np.array([[[0. for i in range(12)] for j in range(12)] for n in range(p)])
gamma12=np.array([[[0. for i in range(12)] for j in range(12)] for n in range(p)])
t0=np.load('influence_tensors-nu=0.0-patch_size=02-max_depth=10.npy')
t12=np.load('influence_tensors-nu=0.5-patch_size=02-max_depth=10.npy')
for n in range(p):
for j in range(4):
for k in range(3):
for i in range(4):
for l in range(3):
gamma0[n,3*i+l,3*j+k]=t0[9-n,i,j,l,k]
gamma12[n,3*i+l,3*j+k]=t12[9-n,i,j,l,k]
#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):
M=np.array([[0. for i in range(12)] for j in range(12)])
for i in range(4):
S=np.linalg.inv(C_list[i]-K)
for k in range(3):
for l in range(3):
M[i*3+l,3*i+k]=S[l,k]
gamma_=gamma[n-1]
T=np.linalg.inv(M-gamma_)
Chom=np.array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
MT_=np.array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
T_=np.array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
E=np.array([np.array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]]) for i in range(4)])
MT=np.dot(M,T)
for k in range(3):
for l in range(3):
for j in range(4):
for i in range(4):
E[i,l,k]+=MT[3*i+l,3*j+k]
T_[l,k]+=0.25*T[3*i+l,3*j+k]
MT_[l,k]+=0.25*MT[3*i+l,3*j+k]
for i in range(4):
E[i]=np.dot(E[i],np.linalg.inv(MT_))
Chom=K+np.dot(T_,np.linalg.inv(MT_))
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):
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 j in range(2**(n-1)):
C_list=[tC_[2*i+1,2*j],tC_[2*i+1,2*j+1],tC_[2*i,2*j],tC_[2*i,2*j+1]]
t1[i,j],Aloc=auxiliary_problem(C_list, n,gamma,K)
for ii in range(2**(p-n)):
for jj in range(2**(p-n)):
tA[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj]=np.dot(tA[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj],Aloc[0])
for ii in range(2**(p-n)):
for jj in range(2**(p-n)):
tA[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj+2**(p-n)]=np.dot(tA[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj+2**(p-n)],Aloc[1])
for ii in range(2**(p-n)):
for jj in range(2**(p-n)):
tA[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj]=np.dot(tA[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj],Aloc[2])
for ii in range(2**(p-n)):
for jj in range(2**(p-n)):
tA[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj+2**(p-n)]=np.dot(tA[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj+2**(p-n)],Aloc[3])
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_):
# K[2,2]=K[0,0]-K[0,1]
# mu=np.matrix.trace(np.dot(K,K0))/4
# km=np.matrix.trace(np.dot(K,J0))/3
# nu=(3*km-2*mu)/2/(3*km)
C_11=(K_[0,0]+K_[1,1])/2
C_12=(K_[0,1]+K_[1,0])/2
nu=1/(C_11/C_12+1)
mu=C_11/2*(1-2*nu)/(1-nu)
C11=2*mu*(1-nu)/(1-2*nu)
C22=C11
C12=nu/(1-nu)*C11
C33=2*mu
K=np.array([[C11,C12,0.],[C12,C22,0.],[0.,0.,C33]])
gamma=-((1-2*nu)*gamma0+nu*gamma12)/(1-nu)/mu
t=tC0
tA=np.array([[np.array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]]) for i in range(2**p)] for j in range(2**p)])
for i in range(n):
t2=tC(t,tA,n-i,gamma,K)
t=t2
return(K,t[0,0],tA)
C_hom1=CVoigt
err=1.0
it=1
C_hom_list=[]
tA_list=[]
while err > 1e-2:
print("itération "+str(it)+'\n')
C_hom=C_hom1
C0,C_hom1,tA=Chom(tabC,p,C_hom)
err=abs(C_hom1[0,0]-C_hom[0,0])/C_hom[0,0]
it+=1
C_hom_list.append(np.array([C0,C_hom1]))
tA_list.append(tA)
print(C0,C_hom1)
return C_hom_list,tA_list
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
#cf 'Contenu des fichiers' à propos de 'reference'
def reference(tabC,res1,c,ec,fr):
#res=256
k,ll,mm,nn=tabC.shape
mesh = UnitSquareMesh(res1, res1, "crossed")
class LocalC(UserExpression):
def __init__(self, Cij, mesh):
self.Cij = Cij
self.Nx, self.Ny, a, b = Cij.shape
assert a==3 and b==3, "Wrong shape"
self.mesh = mesh
UserExpression.__init__(self)
def eval_cell(self, value, x, ufc_cell):
p = Cell(mesh, ufc_cell.index).midpoint()
i, j = int(p[0]*(self.Nx)), int(p[1]*(self.Ny))
value[:] = self.Cij[j, i,:,:].flatten()
def value_shape(self):
return (3, 3)
#Ici on crée, à partir du tableau des tenseurs d'élasticité correspondant au domaine dont on veut calculer le Chom, un tableau Cv représentant le même Cn mais adapté à une résolution plus fine 'res'
# Cv = np.array([[np.zeros((3,3)) for i in range(res)] for j in range(res)])
# for i0 in range(k):
# for j0 in range(k):
# for i in range(int(i0*res/k),int((i0+1)*res/k)):
# for j in range (int(j0*res/k),int((j0+1)*res/k)):
# Cv[i,j]=tabC[i0,j0]
# plt.imshow(tabC[:,:,0,0])
# plt.show()
# plt.close()
# plt.imshow(Cv[:,:,0,0])
# plt.show()
# plt.close()
#Cvar est la transformation du tableau Cv en une fonction définie sur 'mesh'
Cvar= LocalC(tabC,mesh)
# class used to define the periodic boundary map
class PeriodicBoundary(SubDomain):
def __init__(self, vertices, tolerance=DOLFIN_EPS):
""" vertices stores the coordinates of the 4 unit cell corners"""
SubDomain.__init__(self, tolerance)
self.tol = tolerance
self.vv = vertices
self.a1 = self.vv[1,:]-self.vv[0,:] # first vector generating periodicity
self.a2 = self.vv[3,:]-self.vv[0,:] # second vector generating periodicity
# check if UC vertices form indeed a parallelogram
assert np.linalg.norm(self.vv[2, :]-self.vv[3, :] - self.a1) <= self.tol
assert np.linalg.norm(self.vv[2, :]-self.vv[1, :] - self.a2) <= self.tol
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the
# bottom-right or top-left vertices
return bool((near(x[0], self.vv[0,0] + x[1]*self.a2[0]/self.vv[3,1], self.tol) or
near(x[1], self.vv[0,1] + x[0]*self.a1[1]/self.vv[1,0], self.tol)) and
(not ((near(x[0], self.vv[1,0], self.tol) and near(x[1], self.vv[1,1], self.tol)) or
(near(x[0], self.vv[3,0], self.tol) and near(x[1], self.vv[3,1], self.tol)))) and on_boundary)
def map(self, x, y):
if near(x[0], self.vv[2,0], self.tol) and near(x[1], self.vv[2,1], self.tol): # if on top-right corner
y[0] = x[0] - (self.a1[0]+self.a2[0])
y[1] = x[1] - (self.a1[1]+self.a2[1])
elif near(x[0], self.vv[1,0] + x[1]*self.a2[0]/self.vv[2,1], self.tol): # if on right boundary
y[0] = x[0] - self.a1[0]
y[1] = x[1] - self.a1[1]
else: # should be on top boundary
y[0] = x[0] - self.a2[0]
y[1] = x[1] - self.a2[1]
#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 eps(v):
return as_vector([v[0].dx(0),
v[1].dx(1),
(v[0].dx(1)+v[1].dx(0))/sqrt(2)])
def sigma(v, Eps):
return dot(C, eps(v)+Eps)
Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0)
vertices = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))
V = VectorFunctionSpace(mesh,'CG',degree=2,dim=3)
v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
w = Function(W)
Eps = Constant((0,)*3)
F = inner(sigma(dv, Eps), eps(v_))*dx
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
Chom = np.zeros((3, 3))
for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
ee = np.zeros((3,))
ee[j] = 1
Eps.assign(Constant(ee))
solve(a == L, w, [])#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"})
(v, lamb) = split(w)
ffile = XDMFFile('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Strain EF '+case+'.xdmf')
ffile.parameters["functions_share_mesh"] = True
vf = Function(V, name='Strain-'+case)
vf.assign(project(eps(v),V))
ffile.write(vf,0.)
plt.figure()
pp=plot(eps(v)[0])
plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 1'+str(j+1)+'.pdf')
#plt.show()
plt.close()
plt.figure()
pp=plot(eps(v)[1])
plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 2'+str(j+1)+'.pdf')
#plt.show()
plt.close()
plt.figure()
pp=plot(eps(v)[2])
plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf')
#plt.show()
plt.close()
Sigma = np.zeros((3,))
for kk in range(3):
Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx)
Chom[j, :] = Sigma
ffile.close()
return Chom
Chom = reference_problem(Cvar, mesh,k)
CVoigt = np.array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
for i in range(k):
for j in range(k):
CVoigt=CVoigt+tabC[i,j]
CVoigt=CVoigt/k/k
print(Chom)
return Chom,CVoigt
import numpy as np
import h5py
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
def verif_CV(c,i,j):
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
f=fi['fractions'][j]
ChomEF=fi['ChomEF'][j]
ChomHH=fi['ChomHH'][j]
niter=len(ChomHH)
iterations=np.zeros((niter,))
C_11=np.zeros((niter,))
C_22=np.zeros((niter,))
C_12=np.zeros((niter,))
C_33=np.zeros((niter,))
C_11_EF=ChomEF[0,0]
C_22_EF=ChomEF[1,1]
C_12_EF=ChomEF[0,1]
C_33_EF=ChomEF[2,2]/2
for k in range(niter):
iterations[k]=k+1
C_11[k]=ChomHH[k,1,0,0]
C_22[k]=ChomHH[k,1,1,1]
C_12[k]=ChomHH[k,1,0,1]
C_33[k]=ChomHH[k,1,2,2]/2
plt.figure()
g1,=plt.plot(iterations,C_11,color='red')
plt.axhline(y=C_11_EF,color='red')
g2,=plt.plot(iterations,C_22,color='blue')
plt.axhline(y=C_22_EF,color='blue')
g3,=plt.plot(iterations,C_12,color='yellow')
plt.axhline(y=C_12_EF,color='yellow')
g4,=plt.plot(iterations,C_33,color='green')
plt.axhline(y=C_33_EF,color='green')
plt.legend([g1,g2,g3,g4],['C_11','C_22','C_12','C_33'])
plt.title('Contraste '+str(c)+' échantillon '+str(i)+' fraction '+str(f))
plt.xlabel('itérations')
plt.show()
plt.close()
def verif_im(c,i):
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
images=fi['images']
l=len(images)
for j in range(l):
f=fi['fractions'][j]
plt.figure()
plt.imshow(images[j,:,:],origin='lower')
plt.savefig('c='+str(c)+'ech='+str(i)+'frac='+str(f)+'.pdf')
plt.show()
plt.close()
# for i in range(2):
# for j in range(3):
# verif_CV(100,i,j)
#verif_im(100,0)
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]
f=frac[j]
plt.figure()
pp=plt.imshow(EpsHH[:,:,0,0],origin='lower')
plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/Aloc HH 11 iter='+str(niter)+'.pdf')
#plt.show()
plt.close()
plt.figure()
pp=plt.imshow(EpsHH[:,:,0,1],origin='lower')
plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/Aloc HH 12 iter='+str(niter)+'.pdf')
#plt.show()
plt.close()
plt.figure()
pp=plt.imshow(EpsHH[:,:,1,1],origin='lower')
plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/Aloc HH 22 iter='+str(niter)+'.pdf')
#plt.show()
plt.close()
plt.figure()
pp=plt.imshow(EpsHH[:,:,1,1],origin='lower')
plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/Aloc HH 33 iter='+str(niter)+'.pdf')
#plt.show()
plt.close()
#verif_eps(100,0,1,10)
#verif_eps(100,0,3,10)
for c in [100,0.01]:
print("contraste "+str(c)+'\n')
for i in range(10):
print("échantillon "+str(i)+'\n')
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r")
fractions=fi['fractions']
ChomEF=fi['ChomEF']
ChomHH=fi['ChomHH']
C_11=np.zeros((8,))
C_11_EF=np.zeros((8,))
C_22=np.zeros((8,))
C_22_EF=np.zeros((8,))
C_12=np.zeros((8,))
C_12_EF=np.zeros((8,))
C_33=np.zeros((8,))
C_33_EF=np.zeros((8,))
for j,f in enumerate([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8]):
C_11[j]=ChomHH[j,-1,1,0,0]/1.11111*np.sqrt(c)/10
C_11_EF[j]=ChomEF[j,0,0]/1.11111*np.sqrt(c)/10
C_22[j]=ChomHH[j,-1,1,1,1]/1.11111*np.sqrt(c)/10
C_22_EF[j]=ChomEF[j,1,1]/1.11111*np.sqrt(c)/10
C_12[j]=ChomHH[j,-1,1,0,1]/0.277777*np.sqrt(c)/10
C_12_EF[j]=ChomEF[j,0,1]/0.277777*np.sqrt(c)/10
C_33[j]=ChomHH[j,-1,1,2,2]/2/0.416666*np.sqrt(c)/10
C_33_EF[j]=ChomEF[j,2,2]/2/0.416666*np.sqrt(c)/10
plt.figure()
g1,=plt.plot(fractions,C_11)
g11,=plt.plot(fractions,C_11_EF)
plt.legend([g1,g11],['HH','EF'])
plt.title('C_11 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1])
plt.show()
plt.close()
plt.figure()
g1,=plt.plot(fractions,C_22)
g11,=plt.plot(fractions,C_22_EF)
plt.legend([g1,g11],['HH','EF'])
plt.title('C_22 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1])
plt.show()
plt.close()
plt.figure()
g1,=plt.plot(fractions,C_12)
g11,=plt.plot(fractions,C_12_EF)
plt.legend([g1,g11],['HH','EF'])
plt.title('C_12 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1])
plt.show()
plt.close()
plt.figure()
g1,=plt.plot(fractions,C_33)
g11,=plt.plot(fractions,C_33_EF)
plt.legend([g1,g11],['HH','EF'])
plt.title('C_33 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1])
plt.show()
plt.close()
\ No newline at end of file
import numpy as np
import h5py
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from Reference import reference
from Homog_hierar import homog_hierar
from inclusions import image_init,image_augmente
import os
N=256
D=0.03
C=[0. for i in range(2)]
E=[0.,0.]
nu=[0.,0.]
E[0]=1
nu[0]=0.2
nu[1]=0.2
for c in [100,0.01]:
print("contraste "+str(c)+'\n')
os.system("mkdir c="+str(c))
E[1]=1/c
for i in range(2):
C[i]=np.array([[(1-nu[i])*E[i]/(1-2*nu[i])/(1+nu[i]),
nu[i]*E[i]/(1-2*nu[i])/(1+nu[i]),0.],
[nu[i]*E[i]/(1-2*nu[i])/(1+nu[i]),(1-nu[i])*E[i]/(1-2*nu[i])/(1+nu[i]),0.],
[0.,0.,E[i]/(1+nu[i])]])
for i in range(10):
os.system("mkdir c="+str(c)+"/ech="+str(i))
print("échantillon "+str(i)+'\n')
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","a")
im,f0=image_init(N,D,0.1)
fi['images']=np.zeros((8,N,N))
fi['fractions']=np.zeros((8,))
fi['ChomEF']=np.zeros((8,3,3))
fi['CVoigt']=np.zeros((8,3,3))
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]):
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)
else:
f1=f0
# plt.figure()
# plt.imshow(im[:,:],origin='lower')
# plt.show()
# plt.close()
fi['images'][j]=im
fi['fractions'][j]=f1
tC=np.array([[C[0] for ii in range(N)] for jj in range(N)])
for ii in range(N):
for jj in range(N):
if im[ii,jj]==1:
tC[ii,jj]=C[1]
ChomEF,CVoigt = reference(tC,128,c,i,f)
fi['ChomEF'][j]=ChomEF
fi['CVoigt'][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]
for k in range(l,10):
fi['ChomHH'][j,k]=ChomHH_list[l-1]
fi['EpsHH'][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]