Commit 373457c0 authored by Antoine MARTIN's avatar Antoine MARTIN
parents a007729c beb06cad
...@@ -10,11 +10,22 @@ from matplotlib.ticker import LinearLocator, FormatStrFormatter ...@@ -10,11 +10,22 @@ from matplotlib.ticker import LinearLocator, FormatStrFormatter
def local_project(v,V):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv,v_)*dx
b_proj = inner(v,v_)*dx
solver = LocalSolver(a_proj,b_proj)
solver.factorize()
u = Function(V)
solver.solve_local_rhs(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):
#res=256 #res=256
k,ll,mm,nn=tabC.shape k,ll,mm,nn=tabC.shape
mesh = UnitSquareMesh(res1, res1, "crossed") mesh = UnitSquareMesh(res1, res1, "left/right")
class LocalC(UserExpression): class LocalC(UserExpression):
def __init__(self, Cij, mesh): def __init__(self, Cij, mesh):
...@@ -32,7 +43,6 @@ def reference(tabC,res1,c,ec,fr): ...@@ -32,7 +43,6 @@ def reference(tabC,res1,c,ec,fr):
return (3, 3) 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' #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)]) # Cv = np.array([[np.zeros((3,3)) for i in range(res)] for j in range(res)])
# for i0 in range(k): # for i0 in range(k):
...@@ -101,7 +111,7 @@ def reference(tabC,res1,c,ec,fr): ...@@ -101,7 +111,7 @@ def reference(tabC,res1,c,ec,fr):
Re = VectorElement("R", mesh.ufl_cell(), 0) Re = VectorElement("R", mesh.ufl_cell(), 0)
vertices = np.array([[0, 0], [1, 0], [1, 1], [0, 1]]) vertices = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10)) W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))
V = VectorFunctionSpace(mesh,'CG',degree=2,dim=3) V = VectorFunctionSpace(mesh,'DG',degree=1,dim=3)
v_,lamb_ = TestFunctions(W) v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W) dv, dlamb = TrialFunctions(W)
...@@ -114,42 +124,42 @@ def reference(tabC,res1,c,ec,fr): ...@@ -114,42 +124,42 @@ def reference(tabC,res1,c,ec,fr):
Chom = np.zeros((3, 3)) Chom = np.zeros((3, 3))
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))
ee = np.zeros((3,)) ee = np.zeros((3,))
ee[j] = 1 ee[j] = 1
Eps.assign(Constant(ee)) Eps.assign(Constant(ee))
solve(a == L, w, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"}) solve(a == L, w, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"})
(v, lamb) = split(w) (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 = Function(V, name='Strain-'+case)
vf.assign(project(eps(v),V)) vf.assign(local_project(eps(v), V))
ffile.write(vf,0.) ffile.write_checkpoint(vf, case, 0.)
plt.figure() # plt.figure()
pp=plot(eps(v)[0]) # pp=plot(vf[0])
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 1'+str(j+1)+'.pdf') # plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 1'+str(j+1)+'.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plot(eps(v)[1]) # pp=plot(vf[1])
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 2'+str(j+1)+'.pdf') # plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 2'+str(j+1)+'.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
plt.figure() # plt.figure()
pp=plot(eps(v)[2]) # pp=plot(vf[2])
plt.colorbar(pp) # plt.colorbar(pp)
plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf') # plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf')
#plt.show() # #plt.show()
plt.close() # plt.close()
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 return Chom
Chom = reference_problem(Cvar, mesh,k) Chom = reference_problem(Cvar, mesh,k)
......
...@@ -2,7 +2,7 @@ import numpy as np ...@@ -2,7 +2,7 @@ import numpy as np
import h5py import h5py
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib import cm from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter from matplotlib.ticker import LinearLocator, FormatStrFormatter
...@@ -27,7 +27,7 @@ def verif_CV(c,i,j): ...@@ -27,7 +27,7 @@ def verif_CV(c,i,j):
C_22[k]=ChomHH[k,1,1,1] C_22[k]=ChomHH[k,1,1,1]
C_12[k]=ChomHH[k,1,0,1] C_12[k]=ChomHH[k,1,0,1]
C_33[k]=ChomHH[k,1,2,2]/2 C_33[k]=ChomHH[k,1,2,2]/2
plt.figure() plt.figure()
g1,=plt.plot(iterations,C_11,color='red') g1,=plt.plot(iterations,C_11,color='red')
plt.axhline(y=C_11_EF,color='red') plt.axhline(y=C_11_EF,color='red')
...@@ -38,11 +38,11 @@ def verif_CV(c,i,j): ...@@ -38,11 +38,11 @@ def verif_CV(c,i,j):
g4,=plt.plot(iterations,C_33,color='green') g4,=plt.plot(iterations,C_33,color='green')
plt.axhline(y=C_33_EF,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.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.title('Contraste '+str(c)+' échantillon '+str(i)+' fraction '+str(f))
plt.xlabel('itérations') plt.xlabel('itérations')
plt.show() plt.show()
plt.close() plt.close()
def verif_im(c,i): def verif_im(c,i):
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")
images=fi['images'] images=fi['images']
...@@ -54,12 +54,12 @@ def verif_im(c,i): ...@@ -54,12 +54,12 @@ def verif_im(c,i):
plt.savefig('c='+str(c)+'ech='+str(i)+'frac='+str(f)+'.pdf') plt.savefig('c='+str(c)+'ech='+str(i)+'frac='+str(f)+'.pdf')
plt.show() plt.show()
plt.close() plt.close()
# for i in range(2): # for i in range(2):
# for j in range(3): # for j in range(3):
# verif_CV(100,i,j) # verif_CV(100,i,j)
#verif_im(100,0) #verif_im(100,0)
def verif_eps(c,i,j,niter): def verif_eps(c,i,j,niter):
...@@ -91,14 +91,14 @@ def verif_eps(c,i,j,niter): ...@@ -91,14 +91,14 @@ def verif_eps(c,i,j,niter):
plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/Aloc HH 33 iter='+str(niter)+'.pdf') plt.savefig('c='+str(c)+'/ech='+str(i)+'/frac='+str(f)+'/Aloc HH 33 iter='+str(niter)+'.pdf')
#plt.show() #plt.show()
plt.close() plt.close()
#verif_eps(100,0,1,10) #verif_eps(100,0,1,10)
#verif_eps(100,0,3,10) #verif_eps(100,0,3,10)
for c in [100,0.01]: for c in [100,0.01]:
print("contraste "+str(c)+'\n') print("contraste "+str(c)+'\n')
for i in range(10): for i in range(10):
print("échantillon "+str(i)+'\n') print("échantillon "+str(i)+'\n')
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")
fractions=fi['fractions'] fractions=fi['fractions']
...@@ -125,47 +125,46 @@ for c in [100,0.01]: ...@@ -125,47 +125,46 @@ for c in [100,0.01]:
g1,=plt.plot(fractions,C_11) g1,=plt.plot(fractions,C_11)
g11,=plt.plot(fractions,C_11_EF) g11,=plt.plot(fractions,C_11_EF)
plt.legend([g1,g11],['HH','EF']) plt.legend([g1,g11],['HH','EF'])
plt.title('C_11 contraste '+str(c)+' échantillon '+str(i)) plt.title('C_11 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f') plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide') plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1]) plt.axis([0,1,0,1])
plt.show() plt.show()
plt.close() plt.close()
plt.figure() plt.figure()
g1,=plt.plot(fractions,C_22) g1,=plt.plot(fractions,C_22)
g11,=plt.plot(fractions,C_22_EF) g11,=plt.plot(fractions,C_22_EF)
plt.legend([g1,g11],['HH','EF']) plt.legend([g1,g11],['HH','EF'])
plt.title('C_22 contraste '+str(c)+' échantillon '+str(i)) plt.title('C_22 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f') plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide') plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1]) plt.axis([0,1,0,1])
plt.show() plt.show()
plt.close() plt.close()
plt.figure() plt.figure()
g1,=plt.plot(fractions,C_12) g1,=plt.plot(fractions,C_12)
g11,=plt.plot(fractions,C_12_EF) g11,=plt.plot(fractions,C_12_EF)
plt.legend([g1,g11],['HH','EF']) plt.legend([g1,g11],['HH','EF'])
plt.title('C_12 contraste '+str(c)+' échantillon '+str(i)) plt.title('C_12 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f') plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide') plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1]) plt.axis([0,1,0,1])
plt.show() plt.show()
plt.close() plt.close()
plt.figure() plt.figure()
g1,=plt.plot(fractions,C_33) g1,=plt.plot(fractions,C_33)
g11,=plt.plot(fractions,C_33_EF) g11,=plt.plot(fractions,C_33_EF)
plt.legend([g1,g11],['HH','EF']) plt.legend([g1,g11],['HH','EF'])
plt.title('C_33 contraste '+str(c)+' échantillon '+str(i)) plt.title('C_33 contraste '+str(c)+' échantillon '+str(i))
plt.xlabel('f') plt.xlabel('f')
plt.ylabel('Cij/Cij_rigide') plt.ylabel('Cij/Cij_rigide')
plt.axis([0,1,0,1]) plt.axis([0,1,0,1])
plt.show() plt.show()
plt.close() plt.close()
\ No newline at end of file
...@@ -36,7 +36,8 @@ for c in [100,0.01]: ...@@ -36,7 +36,8 @@ for c in [100,0.01]:
for i in range(10): for i in range(10):
os.system("mkdir c="+str(c)+"/ech="+str(i)) os.system("mkdir c="+str(c)+"/ech="+str(i))
print("échantillon "+str(i)+'\n') print("échantillon "+str(i)+'\n')
fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","a") 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,f0=image_init(N,D,0.1)
fi['images']=np.zeros((8,N,N)) fi['images']=np.zeros((8,N,N))
fi['fractions']=np.zeros((8,)) fi['fractions']=np.zeros((8,))
......
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