diff --git a/Reference.py b/Reference.py index 4b7b2d1f1f26475cd6b43b776e84f6e93c25f05a..f5b5d3b1fe64e7d7acb496e4c17962bb0654a91c 100644 --- a/Reference.py +++ b/Reference.py @@ -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' def reference(tabC,res1,c,ec,fr): #res=256 k,ll,mm,nn=tabC.shape - mesh = UnitSquareMesh(res1, res1, "crossed") + mesh = UnitSquareMesh(res1, res1, "left/right") class LocalC(UserExpression): def __init__(self, Cij, mesh): @@ -32,7 +43,6 @@ def reference(tabC,res1,c,ec,fr): 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): @@ -101,7 +111,7 @@ def reference(tabC,res1,c,ec,fr): 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 = VectorFunctionSpace(mesh,'DG',degree=1,dim=3) v_,lamb_ = TestFunctions(W) dv, dlamb = TrialFunctions(W) @@ -114,42 +124,42 @@ def reference(tabC,res1,c,ec,fr): Chom = np.zeros((3, 3)) 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[j] = 1 Eps.assign(Constant(ee)) solve(a == L, w, solver_parameters={"linear_solver": "mumps"})#, 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() + vf.assign(local_project(eps(v), V)) + ffile.write_checkpoint(vf, case, 0.) + + # plt.figure() + # pp=plot(vf[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(vf[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(vf[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() + ffile.close() return Chom Chom = reference_problem(Cvar, mesh,k) diff --git a/courbes.py b/courbes.py index 09c07d3b39747ce773ebc544f3fccb522f4a224c..3026ff1ec75cf546cace1d1d045a392b4fbc3260 100644 --- a/courbes.py +++ b/courbes.py @@ -2,7 +2,7 @@ import numpy as np import h5py import matplotlib.pyplot as plt 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): 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') @@ -38,11 +38,11 @@ def verif_CV(c,i,j): 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.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'] @@ -54,12 +54,12 @@ def verif_im(c,i): 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): @@ -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.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): + 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'] @@ -125,47 +125,46 @@ for c in [100,0.01]: 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.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.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.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.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 + + + + diff --git a/execute.py b/execute.py index ed475b9b38120e3c7d1a08a0aa1deb9f2f61621f..d2aa11fb9535f18b943dd68ff765adcb2152c849 100644 --- a/execute.py +++ b/execute.py @@ -36,7 +36,8 @@ for c in [100,0.01]: 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") + fname = "c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5" + fi = h5py.File(fname, "w") im,f0=image_init(N,D,0.1) fi['images']=np.zeros((8,N,N)) fi['fractions']=np.zeros((8,))