From 19f9155deb67c5b6155cd27146bc6d6ef68cdad1 Mon Sep 17 00:00:00 2001 From: Antoine Martin Date: Mon, 10 Aug 2020 19:10:58 +0200 Subject: [PATCH] add explo2 explohh --- explo2/.vi.swp | Bin 0 -> 12288 bytes explo2/Champs_bis.py | 155 ++ explo2/Homog_hierar.py | 205 +++ explo2/Reference.py | 209 +++ explo2/courbes.py | 407 +++++ explo2/execute.py | 127 ++ explo2/inclusions.py | 132 ++ explo2/inclusions_elliptiques.py | 144 ++ ...sors-nu=0.0-patch_size=02-max_depth=10.npy | Bin 0 -> 18128 bytes ...sors-nu=0.5-patch_size=02-max_depth=10.npy | Bin 0 -> 18128 bytes explo2/output.log | 359 +++++ explohh/.vi.swp | Bin 0 -> 12288 bytes explohh/Aloc.pdf | Bin 0 -> 10233 bytes explohh/Champs_admissibles.py | 186 +++ explohh/Homog_hierar.py | 204 +++ explohh/Reference.py | 209 +++ explohh/courbes.py | 407 +++++ explohh/execute.py | 132 ++ explohh/inclusions.py | 132 ++ explohh/inclusions_elliptiques.py | 144 ++ ...sors-nu=0.0-patch_size=02-max_depth=10.npy | Bin 0 -> 18128 bytes ...sors-nu=0.5-patch_size=02-max_depth=10.npy | Bin 0 -> 18128 bytes explohh/output.log | 1331 +++++++++++++++++ 23 files changed, 4483 insertions(+) create mode 100644 explo2/.vi.swp create mode 100644 explo2/Champs_bis.py create mode 100644 explo2/Homog_hierar.py create mode 100644 explo2/Reference.py create mode 100644 explo2/courbes.py create mode 100644 explo2/execute.py create mode 100644 explo2/inclusions.py create mode 100644 explo2/inclusions_elliptiques.py create mode 100644 explo2/influence_tensors-nu=0.0-patch_size=02-max_depth=10.npy create mode 100644 explo2/influence_tensors-nu=0.5-patch_size=02-max_depth=10.npy create mode 100644 explo2/output.log create mode 100644 explohh/.vi.swp create mode 100644 explohh/Aloc.pdf create mode 100644 explohh/Champs_admissibles.py create mode 100644 explohh/Homog_hierar.py create mode 100644 explohh/Reference.py create mode 100644 explohh/courbes.py create mode 100644 explohh/execute.py create mode 100644 explohh/inclusions.py create mode 100644 explohh/inclusions_elliptiques.py create mode 100644 explohh/influence_tensors-nu=0.0-patch_size=02-max_depth=10.npy create mode 100644 explohh/influence_tensors-nu=0.5-patch_size=02-max_depth=10.npy create mode 100644 explohh/output.log diff --git a/explo2/.vi.swp b/explo2/.vi.swp new file mode 100644 index 0000000000000000000000000000000000000000..5cc85bfb757c71a9396852008b6ba1ca295b63e6 GIT binary patch literal 12288 zcmeI%y$*sf5Ww-n-Dva$PFflNYboe(ZZMo)p z`E5-DcfZ})o24hZjmyu!y+aX|sB;mg^{>`3%&x(uHQEFoRYK=*&rewtIk&h8GD84? zo&|>f+FQBhA&dFUo=*Ip&m%no2q1s}0tg_000Ib@Bv9m*Jas{Xswk^vG^u84K>z^+ z5I_I{1Q0*~0R#|00D*1-O26~JMj{{G%lH5K{{LHdOT}&k5I_I{1Q0*~0R#|0009IL JFj3%`$`=K}7?S`1 literal 0 HcmV?d00001 diff --git a/explo2/Champs_bis.py b/explo2/Champs_bis.py new file mode 100644 index 0000000..8229a03 --- /dev/null +++ b/explo2/Champs_bis.py @@ -0,0 +1,155 @@ + + + +from dolfin import * +import numpy as np +import matplotlib +matplotlib.use('Agg') +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 champs_admissibles(tabC,tab_tau,C0,C_HH,tA,res1): + k,ll,mm,nn=tab_tau.shape + tabChh=np.zeros((k,k,3,3)) + tabC0=np.zeros((k,k,3,3)) + for i in range(k): + for j in range(k): + tabC0[i,j]=C0 + tabChh[i,j]=C_HH + mesh = UnitSquareMesh(res1, res1, "left/right") + + 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) + + Cvar=LocalC(tabC,mesh) + C0var=LocalC(tabC0,mesh) + CHH=LocalC(tabChh,mesh) + tau_var= LocalC(tab_tau,mesh) + AHH= LocalC(tA,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] + + + def reference_problem(C,tau,C0,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 sigma1(v1): + Eps1=as_vector([1.,0.,0.]) + return dot(C0,eps(v1))+dot(tau,Eps1) + + def sigma2(v2): + Eps2=as_vector([0.,1.,0.]) + return dot(C0,eps(v2))+dot(tau,Eps2) + + def sigma3(v3): + Eps3=as_vector([0.,0.,1.]) + return dot(C0,eps(v3))+dot(tau,Eps3) + + def A_CA(e1,e2,e3): + return as_matrix([[e1[0],e2[0],e3[0]],[e1[1],e2[1],e3[1]],[e1[2],e2[2],e3[2]]]) + + def B_SA(A): + return dot(dot(C0,A-AHH)+dot(C,AHH),inv(CHH)) + + def t(A): + return as_matrix([[A[0,0],A[1,0],A[2,0]],[A[0,1],A[1,1],A[2,1]],[A[0,2],A[1,2],A[2,2]]]) + + def E_CA(A): + return dot(t(A),dot(C,A)) + + def E_SA(B): + return dot(t(B),dot(inv(C),B)) + + 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,Ve,Ve, Re,Re,Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10)) + V = VectorFunctionSpace(mesh,'DG',degree=1,dim=3) + + v_1,v_2,v_3,lamb_1,lamb_2,lamb_3 = TestFunctions(W) + dv1,dv2,dv3, dlamb1,dlamb2,dlamb3 = TrialFunctions(W) + w = Function(W) + + + F = inner(sigma1(dv1),eps(v_1))*dx+inner(sigma2(dv2),eps(v_2))*dx+inner(sigma3(dv3),eps(v_3))*dx + a, L = lhs(F), rhs(F) + a += dot(lamb_1,dv1)*dx + dot(dlamb1,v_1)*dx+dot(lamb_2,dv2)*dx+dot(dlamb2,v_2)*dx+dot(lamb_3,dv3)*dx+dot(dlamb3,v_3)*dx + solve(a == L, w)#, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"}) + (v1,v2,v3,lamb1,lamb2,lamb3) = split(w) + e1=eps(v1)+as_vector([1.,0.,0.]) + e2=eps(v2)+as_vector([0.,1.,0.]) + e3=eps(v3)+as_vector([0.,0.,1.]) + A=A_CA(e1,e2,e3) + B=B_SA(A) + E1=E_CA(A) + E2=E_SA(B) + C_CA=np.zeros((3,3)) + S_SA=np.zeros((3,3)) + for i in range(3): + for j in range(3): + C_CA[i,j]=assemble(E1[i,j]*dx) + S_SA[i,j]=assemble(E2[i,j]*dx) + return C_CA,np.linalg.inv(S_SA) + + C_CA,C_SA = reference_problem(Cvar,tau_var,C0var,k) + + + print(C_CA,C_SA) + return C_CA,C_SA + + + + + + diff --git a/explo2/Homog_hierar.py b/explo2/Homog_hierar.py new file mode 100644 index 0000000..e59f5c8 --- /dev/null +++ b/explo2/Homog_hierar.py @@ -0,0 +1,205 @@ + + +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,scheme,schemebis): + + 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)) + + + 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+1,j+1,l,k] + gamma12[n,3*i+l,3*j+k]=t12[9-n,i+1,j+1,l,k] + + + + + + def auxiliary_problem(C_list, n,gamma,K,scheme): + + 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.]]) + Chom_DEF=np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]]) + Chom_CONT=np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]]) + 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)]) + B=np.zeros((4,3,3)) + 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_)) + for i in range(4): + B[i]=np.dot(np.dot(C_list[i],E[i]),np.linalg.inv(Chom_)) + + if scheme=='CONT' or scheme=='DEF': + for i in range(4): + Chom_DEF+=0.25*np.dot(np.transpose(E[i]),np.dot(C_list[i],E[i])) + if scheme=='CONT': + Chom_CONT=np.dot(Chom_,np.dot(np.linalg.inv(Chom_DEF),np.transpose(Chom_))) + + if scheme=='HH': + return Chom_,E,B + if scheme=='DEF': + return Chom_DEF,E,B + if scheme=='CONT': + return Chom_CONT,E,B + + + def tC(tC_,tA,tB,n,gamma,K,scheme): + 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,Bloc=auxiliary_problem(C_list, n,gamma,K,scheme) + 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]) + tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj]=np.dot(tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj],Bloc[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]) + tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj+2**(p-n)]=np.dot(tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj+2**(p-n)],Bloc[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]) + tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj]=np.dot(tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj],Bloc[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]) + tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj+2**(p-n)]=np.dot(tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj+2**(p-n)],Bloc[3]) + return t1 + + + + + def Chom(tC0,n,K_,scheme): + # 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)]) + tB=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,tB,n-i,gamma,K,scheme) + t=t2 + return(K,t[0,0],tA,tB) + + + 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,tB=Chom(tabC,p,C_hom,scheme) + if schemebis=='HHbis': + C_hom2=np.zeros((3,3)) + for i in range(pp): + for j in range(qq): + C_hom2+=np.dot(tabC[i,j],tA[i,j])/pp/qq + if schemebis=='DEFbis': + C_hom2=np.zeros((3,3)) + for i in range(pp): + for j in range(qq): + C_hom2+=np.dot(np.transpose(tA[i,j]),np.dot(tabC[i,j],tA[i,j]))/pp/qq + if schemebis=='CONTbis': + Scont=np.zeros((3,3)) + for i in range(pp): + for j in range(qq): + Scont+=np.dot(np.transpose(tB[i,j]),np.dot(np.linalg.inv(tabC[i,j]),tB[i,j]))/pp/qq + C_hom2=np.linalg.inv(Scont) + C_hom1=C_hom2 + 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(C_hom1) + CA=np.zeros((3,3)) + CDEF=np.zeros((3,3)) + SCONT=np.zeros((3,3)) + + for i in range(pp): + for j in range(qq): + CA+=np.dot(tabC[i,j],tA[i,j])/pp/qq + CDEF+=np.dot(np.transpose(tA[i,j]),np.dot(tabC[i,j],tA[i,j]))/pp/qq + SCONT+=np.dot(np.transpose(tB[i,j]),np.dot(np.linalg.inv(tabC[i,j]),tB[i,j]))/pp/qq + CCONT=np.linalg.inv(SCONT) + + + return CA,CDEF,CCONT,tA + + + + + + + + + + diff --git a/explo2/Reference.py b/explo2/Reference.py new file mode 100644 index 0000000..02676b1 --- /dev/null +++ b/explo2/Reference.py @@ -0,0 +1,209 @@ + + + +from dolfin import * +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from matplotlib import cm +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, "left/right") + + 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) + #Avar= LocalC(tA,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) + + #def sigma_HH(Eps): + # return dot(C,dot(A,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,'DG',degree=1,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"]): + 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) + + vf = Function(V, name='Strain-'+case) + vf.assign(local_project(eps(v), V)) + ffile.write_checkpoint(vf, case, 0.) + + + #sEF=sigma(v,Eps) + #sHH=sigma_HH(Eps) + #diff_sig=sEF-sHH + #plt.figure(figsize=(80.,50.)) + #pp=plot(diff_sig[0],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig_11'+case+'.jpg') + #plt.show() + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(sEF[0],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/sigEF_11'+case+'.jpg') + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(sEF[0]-diff_sig[0],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/sigHH_11'+case+'.jpg') + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(diff_sig[1],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig_22'+case+'.jpg') + #plt.show() + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(diff_sig[2],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig_12'+case+'.jpg') + #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 + + + + + + diff --git a/explo2/courbes.py b/explo2/courbes.py new file mode 100644 index 0000000..8ff286c --- /dev/null +++ b/explo2/courbes.py @@ -0,0 +1,407 @@ +import numpy as np +import h5py +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import cm +from matplotlib.ticker import LinearLocator, FormatStrFormatter +import os + + + +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.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,0.95] + 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[:,:,2,2],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,0,10) +#verif_eps(100,0,3,10) + + +def courbes(c,i): + os.system('mkdir courbes') + fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") + fr=fi['fractions'] + ChomEF=fi['ChomEF'] + ChomHH=fi['ChomHH'] + C_11=np.ones((13,)) + C_11_EF=np.ones((13,)) + C_22=np.ones((13,)) + C_22_EF=np.ones((13,)) + C_12=np.ones((13,)) + C_12_EF=np.ones((13,)) + C_33=np.ones((13,)) + C_33_EF=np.ones((13,)) + C_11_0=np.ones((13,)) + C_22_0=np.ones((13,)) + C_12_0=np.ones((13,)) + C_33_0=np.ones((13,)) + 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,0.95,1.]): + if c==0.01: + jj=12-j + else: + jj=j + if f==0. or f==1.: + if c==0.01: + fractions[jj]=1-f + else: + fractions[jj]=f + if jj==12: + C_11[jj]=0.01 + C_11_0[jj]=0.01 + C_11_EF[jj]=0.01 + C_22[jj]=0.01 + C_22_0[jj]=0.01 + C_22_EF[jj]=0.01 + C_12[jj]=0.01 + C_12_0[jj]=0.01 + C_12_EF[jj]=0.01 + C_33[jj]=0.01 + C_33_0[jj]=0.01 + C_33_EF[jj]=0.01 + if f!=0. and f!=1.: + if c==0.01: + fractions[jj]=1.-fr[j-1] + else: + fractions[jj]=fr[j-1] + C_11[jj]=ChomHH[j-1,-1,1,0,0]/1.11111*np.sqrt(c)/10 + C_11_0[jj]=ChomHH[j-1,-1,0,0,0]/1.11111*np.sqrt(c)/10 + C_11_EF[jj]=ChomEF[j-1,0,0]/1.11111*np.sqrt(c)/10 + C_22[jj]=ChomHH[j-1,-1,1,1,1]/1.11111*np.sqrt(c)/10 + C_22_0[jj]=ChomHH[j-1,-1,0,1,1]/1.11111*np.sqrt(c)/10 + C_22_EF[jj]=ChomEF[j-1,1,1]/1.11111*np.sqrt(c)/10 + C_12[jj]=ChomHH[j-1,-1,1,0,1]/0.277777*np.sqrt(c)/10 + C_12_0[jj]=ChomHH[j-1,-1,0,0,1]/0.277777*np.sqrt(c)/10 + C_12_EF[jj]=ChomEF[j-1,0,1]/0.277777*np.sqrt(c)/10 + C_33[jj]=ChomHH[j-1,-1,1,2,2]/2/0.416666*np.sqrt(c)/10 + C_33_0[jj]=ChomHH[j-1,-1,0,2,2]/2/0.416666*np.sqrt(c)/10 + C_33_EF[jj]=ChomEF[j-1,2,2]/2/0.416666*np.sqrt(c)/10 + plt.figure() + g1,=plt.plot(fractions,C_11,'+') + g2,=plt.plot(fractions,C_11_0,'+') + g11,=plt.plot(fractions,C_11_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_11 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_11 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_22,'+') + g2,=plt.plot(fractions,C_22_0,'+') + g11,=plt.plot(fractions,C_22_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_22 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_22 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_12,'+') + g2,=plt.plot(fractions,C_12_0,'+') + g11,=plt.plot(fractions,C_12_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_12 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_12 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_33,'+') + g2,=plt.plot(fractions,C_33_0,'+') + g11,=plt.plot(fractions,C_33_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_33 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_33 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + +def courbe_moyenne(c,nb): + os.system('mkdir courbes') + fichiers=[] + for i in range(nb): + fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") + fichiers.append(fi) + + C_11=np.ones((13,)) + C_11_EF=np.ones((13,)) + C_22=np.ones((13,)) + C_22_EF=np.ones((13,)) + C_12=np.ones((13,)) + C_12_EF=np.ones((13,)) + C_33=np.ones((13,)) + C_33_EF=np.ones((13,)) + C_11_DEF=np.ones((13,)) + C_22_DEF=np.ones((13,)) + C_12_DEF=np.ones((13,)) + C_33_DEF=np.ones((13,)) + C_11_CONT=np.ones((13,)) + C_22_CONT=np.ones((13,)) + C_12_CONT=np.ones((13,)) + C_33_CONT=np.ones((13,)) + fractions=np.ones((13,)) + for i in range(nb): + ChomHH=fichiers[i]['ChomHH'] + ChomCA=fichiers[i]['C_CA'] + ChomSA=fichiers[i]['C_SA'] + 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,0.95,1.]): + if c==0.01: + jj=12-j + fractions[jj]=1.-f + else: + jj=j + fractions[jj]=f + if jj==12: + C_11[jj]=0.01 + C_11_DEF[jj]=0.01 + C_11_CONT[jj]=0.01 + C_11_EF[jj]=0.01 + C_22[jj]=0.01 + C_22_DEF[jj]=0.01 + C_22_CONT[jj]=0.01 + C_22_EF[jj]=0.01 + C_12[jj]=0.01 + C_12_DEF[jj]=0.01 + C_12_CONT[jj]=0.01 + C_12_EF[jj]=0.01 + C_33[jj]=0.01 + C_33_DEF[jj]=0.01 + C_33_CONT[jj]=0.01 + C_33_EF[jj]=0.01 + if f!=0. and f!=1.: + C_11[jj]+=ChomHH[j-1,0,0,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_11_DEF[jj]+=ChomCA[j-1,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_11_CONT[jj]+=ChomSA[j-1,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_11_EF[jj]+=ChomEF[j-1,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22[jj]+=ChomHH[j-1,0,0,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22_DEF[jj]+=ChomCA[j-1,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22_CONT[jj]+=ChomSA[j-1,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22_EF[jj]+=ChomEF[j-1,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_12[jj]+=ChomHH[j-1,0,0,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_12_DEF[jj]+=ChomCA[j-1,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_12_CONT[jj]+=ChomSA[j-1,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_12_EF[jj]+=ChomEF[j-1,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_33[jj]+=ChomHH[j-1,0,0,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + C_33_DEF[jj]+=ChomCA[j-1,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + C_33_CONT[jj]+=ChomSA[j-1,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + C_33_EF[jj]+=ChomEF[j-1,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + + plt.figure() + g1,=plt.plot(fractions,C_11,'+-') + g2,=plt.plot(fractions,C_11_DEF,'+-') + g3,=plt.plot(fractions,C_11_CONT,'+-') + g11,=plt.plot(fractions,C_11_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_11 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_11 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_22,'+-') + g2,=plt.plot(fractions,C_22_DEF,'+-') + g3,=plt.plot(fractions,C_22_CONT,'+-') + g11,=plt.plot(fractions,C_22_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_22 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_22 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_12,'+-') + g2,=plt.plot(fractions,C_12_DEF,'+-') + g3,=plt.plot(fractions,C_12_CONT,'+-') + g11,=plt.plot(fractions,C_12_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_12 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_12 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_33,'+-') + g2,=plt.plot(fractions,C_33_DEF,'+-') + g3,=plt.plot(fractions,C_33_CONT,'+-') + g11,=plt.plot(fractions,C_33_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_33 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_33 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + +courbe_moyenne(100,1) +#courbe_moyenne(0.01,1) +#for i in range(10): +# courbes(100,i) +# courbes(0.01,i) + +Id=np.array([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]]) + +def err(c,i): + os.system('mkdir courbes') + fi=h5py.File('c='+str(c)+'/c='+str(c)+'ech'+str(i)+'.hdf5','r') + Chom=fi['ChomHH'] + err_def_EF=fi['err_def_EF'] + err_cont_EF=fi['err_cont_EF'] + err_def=fi['err_def'] + err_cont=fi['err_cont'] + fr=fi['fractions'] + err_def_EF_11=np.zeros((11,)) + err_cont_EF_11=np.zeros((11,)) + err_def_11=np.zeros((11,)) + err_cont_11=np.zeros((11,)) + fractions=np.zeros((11,)) + for j,f in enumerate([0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,0.95]): + if c==0.01: + jj=10-j + fractions[10-j]=1-fr[j] + else: + jj=j + fractions[j]=fr[j] + err_def_EF_11[jj]=err_def_EF[j]*np.sqrt(c)*10/1.1111 + err_cont_EF_11[jj]=err_cont_EF[j]/np.sqrt(c)*10/0.96 + err_def_11[jj]=err_def[j][0,0]*np.sqrt(c)*10/1.1111 + err_cont_11[jj]=err_cont[j][0,0]/np.sqrt(c)*10/0.96 + + plt.figure() + g1,=plt.plot(fractions,err_def_11,'+-') + g2,=plt.plot(fractions,err_cont_11,'+-') + plt.legend([g1,g2],['E_HH_souple_11','E_HH_rigide_11']) + plt.title('Energie_HH_11_contraste '+str(c)+'ech '+str(i)) + plt.xlabel('f_souple') + plt.savefig('courbes/Energie_11_contraste '+str(c)+'ech '+str(i)+'.pdf') + plt.close() + + plt.figure() + g1,=plt.plot(fractions,err_def_EF_11,'+-') + g2,=plt.plot(fractions,err_cont_EF_11,'+-') + plt.legend([g1,g2],['E_EF_souple_11','E_EF_rigide_11']) + plt.title('Energie_EF_11_contraste '+str(c)+'ech '+str(i)) + plt.xlabel('f_souple') + plt.savefig('courbes/Energie_EF_11_contraste '+str(c)+'ech '+str(i)+'.pdf') + plt.close() + +#err(100,0) +#err(0.01,0) + + + + + + + + + + + + + diff --git a/explo2/execute.py b/explo2/execute.py new file mode 100644 index 0000000..fd42635 --- /dev/null +++ b/explo2/execute.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +import numpy as np +import h5py +import matplotlib +matplotlib.use('Agg') +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_elliptiques import image_init,image_augmente,image_diminue +from Champs_bis import champs_admissibles +import os + + +N=32 +D=0.1 +nb=1 + +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(nb): + os.system("mkdir c="+str(c)+"/ech="+str(i)) + print("échantillon "+str(i)+'\n') + fname = "c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5" + fi = h5py.File(fname, "w") + im,centers_diam,f1=image_init(N,0.03,0.03,0,True,0.95) + l=len(centers_diam) + fi['centers']=np.zeros((11,l,5)) + fi['images']=np.zeros((11,N,N)) + fi['fractions']=np.zeros((11,)) + fi['ChomEF']=np.zeros((11,3,3)) + fi['CVoigt']=np.zeros((11,3,3)) + fi['C_CA']=np.zeros((11,3,3)) + fi['C_SA']=np.zeros((11,3,3)) + fi['ChomHH']=np.zeros((11,3,3,3,3)) + #fi['EpsHH']=np.zeros((11,10,N,N,3,3)) + fi['Chom_DEF']=np.zeros((11,3,3,3,3)) + fi['Chom_CONT']=np.zeros((11,3,3,3,3)) + + + 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)) + print("fraction volumique "+str(f)+'\n') + if f!=0.95: + f1=image_diminue(im,centers_diam,f) + fi['images'][10-j]=im + l=len(centers_diam) + for ii in range(l): + fi['centers'][10-j][ii]=centers_diam[ii] + fi['fractions'][10-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] + + CA,CDEF,CCONT,tA=homog_hierar(tC,'HH','HHbis') + fi['ChomHH'][10-j,0]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'HH','DEFbis') + # fi['ChomHH'][10-j,1]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'HH','CONTbis') + # fi['ChomHH'][10-j,2]=[CA,CDEF,CCONT] + + # CA,CDEF,CCONT=homog_hierar(tC,'DEF','HHbis') + # fi['Chom_DEF'][10-j,0]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'DEF','DEFbis') + # fi['Chom_DEF'][10-j,1]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'DEF','CONTbis') + # fi['Chom_DEF'][10-j,2]=[CA,CDEF,CCONT] + + # CA,CDEF,CCONT=homog_hierar(tC,'CONT','HHbis') + # fi['Chom_CONT'][10-j,0]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'CONT','DEFbis') + # fi['Chom_CONT'][10-j,1]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'CONT','CONTbis') + # fi['Chom_CONT'][10-j,2]=[CA,CDEF,CCONT] + + ttau=np.zeros((N,N,3,3)) + for ii in range(N): + for jj in range(N): + ttau[ii,jj]=np.dot(tC[ii,jj]-CA,tA[ii,jj]) + + C_CA,C_SA=champs_admissibles(tC,ttau,CA,CA,tA,128) + fi['C_CA'][10-j]=C_CA + fi['C_SA'][10-j]=C_SA + ChomEF,CVoigt = reference(tC,128,c,i,f) + fi['ChomEF'][10-j]=ChomEF + fi['CVoigt'][10-j]=CVoigt + + + + + + + + + + + + + + + + + + diff --git a/explo2/inclusions.py b/explo2/inclusions.py new file mode 100644 index 0000000..d1aceef --- /dev/null +++ b/explo2/inclusions.py @@ -0,0 +1,132 @@ +import numpy as np +import matplotlib.pyplot as plt + + +from random import uniform + + +def compteur(im): + nn,mm=im.shape + cpt=0 + for i in range(nn): + for j in range(mm): + if im[i,j]==1: + cpt=cpt+1 + return cpt/nn/mm + + +def ajoute(im,centers_diam,D,nb): + nn,mm=im.shape + ll=len(centers_diam) + for i in range(nb): + centers_diam.append([uniform(0,1),uniform(0,1),D]) + + for i in range(nb): + r=D/2 + 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 + 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 + 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): + + im=np.zeros((N,N)) + + n=int(f/np.pi/(D/2)**2) + + centers_diam=[] + c=0 + while cf0: + 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 + + + + + + + diff --git a/explo2/inclusions_elliptiques.py b/explo2/inclusions_elliptiques.py new file mode 100644 index 0000000..470bc8d --- /dev/null +++ b/explo2/inclusions_elliptiques.py @@ -0,0 +1,144 @@ +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + + +from random import uniform + + +def compteur(im): + nn,mm=im.shape + cpt=0 + for i in range(nn): + for j in range(mm): + if im[i,j]==1: + cpt=cpt+1 + return cpt/nn/mm + + +def ajoute(im,centers,a,b,th,iso,nb): + nn,mm=im.shape + ll=len(centers) + for i in range(nb): + if iso==True: + centers.append([uniform(0,1),uniform(0,1),uniform(0,np.pi),a,b]) + else: + centers.append([uniform(0,1),uniform(0,1),th,a,b]) + + for i in range(nb): + c=max(centers[ll+i][3],centers[ll+i][4]) + x=centers[ll+i][0] + y=centers[ll+i][1] + the=centers[ll+i][2] + for j in range(int((x-c)*nn)-1,int((x+c)*nn)+1): + for k in range(int((y-c)*mm)-1,int((y+c)*mm)+1): + x1=(j+0.5)/nn + y1=(k+0.5)/mm + x2=(x1-x)*np.cos(the)+(y1-y)*np.sin(the) + y2=(y1-y)*np.cos(the)-(x1-x)*np.sin(the) + if x2**2/a**2+y2**2/b**2<=1: + if j>=nn: + l=nn + else: + l=0 + if k>=mm: + m=mm + else: + m=0 + im[j-l,k-m]=1 + +def enleve(im,centers,nb): + nn,mm=im.shape + im2=np.zeros((nn,mm)) + + for i in range(nb): + centers.pop() + + ll=len(centers) + for i in range(ll): + a=centers[i][3] + b=centers[i][4] + c=max(a,b) + x=centers[i][0] + y=centers[i][1] + the=centers[i][2] + for j in range(int((x-c)*nn)-1,int((x+c)*nn)+1): + for k in range(int((y-c)*mm)-1,int((y+c)*mm)+1): + x1=(j+0.5)/nn + y1=(k+0.5)/mm + x2=(x1-x)*np.cos(the)+(y1-y)*np.sin(the) + y2=(y1-y)*np.cos(the)-(x1-x)*np.sin(the) + if x2**2/a**2+y2**2/b**2<=1: + 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,a,b,th,iso,f): + + im=np.zeros((N,N)) + + n=int(f/np.pi/a/b) + centers=[] + c=0 + while cf0: + print('fraction volumique actuelle :'+str(f0)) + else: + c=f0 + im2=im + while c>f: + l=len(centers) + cc=0. + i=0 + while cc < c-f: + a=centers[l-i-1][3] + b=centers[l-i-1][4] + cc+=np.pi*a*b + i+=1 + im2=enleve(im2,centers,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 + + + + + diff --git a/explo2/influence_tensors-nu=0.0-patch_size=02-max_depth=10.npy b/explo2/influence_tensors-nu=0.0-patch_size=02-max_depth=10.npy new file mode 100644 index 0000000000000000000000000000000000000000..eacec1c37c20863cdd55228dcd346c8b3e90ce49 GIT binary patch literal 18128 zcmeI233yFsyM|Z#H>Gn*l~$X^n1moC2wC)r4ByzjTx zHq7q${4??0{FJfE>%sl{kLo)zxM@(ZIv^@IEGT$D(#X*xdk^cGG_qfRb3L|q;;8<{ z`lvy@hxa!)q<(lZ?^NJbTiH{+xg|O`^MUO+Q^?IN&O(! zsgt`sZpNACd>wH%AkODwb!I>aGuR9c1c?@uT?@`3~yj?+;Ikd^h(i5w&+s#!H=Y z^^Z%k>efh{C`>&sbx$;TZ`!?`<3_Kx))S7Mj;3y2P0ljz^G}>DxR3XB-Ff_*bElR% z$o!m;<14PGbl&KoPG+3XTU0AIHFlaXughYm>8M{PQ->ceoownmW9RGik-jnxnL6*I zP96fEnfF5YGyX<_+PuPQ2YDa-vDW?9MpGxVF7|`{sj_HLrF!%IH8S;G6&!KdPi}Wm zC$nFypMG#Y$kb7XYre>)Zu;ZhuVMGZ)=0coa!bD^4G#I(e&D!+I+=CR$F7)X`bVbD zb?+<9akzpYib3Q765oUbY0 zReE!0{euqbWcG{m!*gJsn~-%ETy?Q-W7(ejWfhe=nRONX&F2Q^`8NLi<@3e7Sy$re zy)(-`f7MkN>wazX=rSGaeea-7W?kIJI9Vd1Q@nqXD-J*Br+o$GpUe6qn!1be{JEuH z8w=+rUb_31^ObNY@NjI~xzYL1$?RAB#&W$<<`>#>_rI@nKa}kZRLuY0a*#7LsRt4C&Pg%jY;bNu zoF5RUvfS?^-v|A?jPImQ-iSCCxz8)kd2j{jd%WV17b8v{_Fch!bCEx}6m_{0L|DPqOhH}%o$KD$t#KDd@84C#=RDt^x_N#SC0xUE z-lfZR{u-|F|8OqE_%%gc`mWB=)efv++CMthNU0;nq1dbr~mI zZExP!f4bm4-p9`&zlzV1>oJJKT`_70#93E!e%jYr?+WwUZQ;B6HtCN}rp|TmtEDfT zt9$#r6>;b@b>2_i#qgPVKjYz#qL9yd_|w4UudJ+r#=7u6_>cRA_l4z)4}V;K<-kAG zBT$bZ)1QCUFU|$}`bz)E)G^1dev&!Q^pU>G&#kJG*#}Mxbu#nee3(~r zet3S!x$uYU-dFnLyI(P`^M&`4>%7Q3crIczo-gLjbHIGb?1NnA`Qomd=PL`(g<`Mg z2AMvy->i#FpVy;)`oXy;(;u#TU+ItUeg$d1`^D!@5T6@bJDhj=!~5tXnV%!qdA{g# z1J@j~AAJ7NK{9o&dtZ6JQcJ^9C(J0eon1F@#PuNs&OHC&em#`A(d5!%srT|v@at-y z?4VB0EceL|%70z9&Aw#rD|CizG~;C3D(qbo)}$cQM!s-G>cxBHI(2f;qh_2U=UBw4 zh&T~4j+6X8;&j0}b#ev7c@J^wBThK>lfOcoidd&ku8TP15vMNVl$HB!9~_VZew8sPmNr$4oEF5(cUw#Mg=?|y}=zRwrmqx3t*b-sorG;8+o)be8Q_FsN}we1;a zxxkZ;1fE>uAcuC7dVu0*{Iol$lQTL_o_z9bC2`1#Q|z3)!;Dj4yMI=dLkBgdjoiAc z8K>C!G}ft;8(MKLJM$1{Fyd5}`<>)4#Ho*U>g10RM?;)Bh!cqYgX7wgo?9S}!H zoZg5ND)-ySk%)5-)~ScG*YXE!2;O#`VdK%@)@>{eI8sjhTYn6gtR}R^;Ublj$?(%<_wK7~fL*XRv2J5i$=$m-$7goWs#x^TWD1 zKlJA+{Gm?f{M7b&zHZkqo-6h%z}2t6cE0%j43J)6h3zojUm$#5sw4sw1B**iT-CIGeFfot%$2-LdZy@;M^++sNmjuf;lb zax27HwqlF6`0YZy!_cQr)u@v#Qui$xaXB|rke5T}y7%?7#p9FL8SlZya_+@<^7(`J zpsJB^Ztz|ULBCa8L$q1)IWHRDziLL$=X3O01L~ZN9K!ctq)JYDdFqAB6TIJpk(ftC zsbk((J6R0JzM39!0#tHCe2z7D_O8eyK*lk+Hsbs?YUKISA18}q^Llr(K12GWlg~lt zy7$%ES3O;)&u!uJdx%3`1fRL?eWlNhkaviCUG`{=JjC6f9XMZxZq;dUDC+XvuaMjH z3+G&6-vYFm=ojm^u7*F;F+b$r;ZF+Isgo^#4S(3LVE4L==MVWd{gQRL^@~1pe(^c3 z`Cj2U&pLI?wN9Uz59h`5??i2l?2%7Qyv)PUt$hYpgU=UHzjZbIL0<)#{&C&=YV9-p zX@v9jt;hKa(a4_X>(+jyd7UrLih-N3Uq@3XcinOo%s#Q7jnbpUB3JEnC81*Wu-65W4M;p*snRLe z?)}<7Le^t%BJ(l$0P5mAzXPAw!cX#G_?&}v>g4M1xxm-_TyE*gr-gey4adrS4BhhC z;Qh#l^-qC6i{Jx!KK#kXI(4$OuUWdj+SmEAN8F}gc)lCD?0a49OPsIP$h*LOe(bF^ zvgi3~EbljT8K<$f8+8r!n&;)n|2Fe1sAHZ*Rrov`^-qF7><4)U^iRC{MSk1T({;{| z^c~NQ$73|IN52}4RIPdp-MSjge(`+CbAmX%!d3E4_wzL}O2#pCD~`ci;E(U~73k{M zrykG4PDxg+wr;5L)8Qara7w-GJ3+?J3l8e! z8Y${a!NIQz`$lu$v-(p%nsI827ea=e-5U0^AZKrodK0WuCwHG}#_6UfSou7w$NnJW z#LE3H&i_Kj>CSar$$^M70r>?qF?{aiteHFB5ksx(Bz=MO zK_~Zvz6|Ts$x=_(=fTHA@cBdRC%+GUCDy5v_rvFsYSN#3XSw;4RY4`MaQEkhW~irv zyc+zVb>Ad`_p)cddI{BMzZBeW=+fs3D*N>Y>VFOXAYVc513gpLC8(39S$evzA#dqB z|6W*C$(7vu71L0YdEotSF#9#lYo49(lk>a?`c|w{C+EXw=JSX9{EQD)$(7ynvEQ)% z-LhinR$c}(pCr^j4E|)n2QuepIo7F@rGKvPQJ(!eA$*^&Ud^RHhHm9!F#F~6eC37s zK3{qGo@(e;9D^_6T+D`#hdt+6kkg^BLS2|=K`wyLufv}S@Mi`3LhcV8^P*EHCtLc> z`BB`?MbNU|D)|Zbe(k7c`C{p2VFq7AU6*{FuWwptKF`;vs?rxjxBM~qp?`M1f*Sbl zS6(A|zoAQ?8@c+mtHbfAMU`iZe+;>}?eT91>5o`HjylMN9i=|b{4w`Il{)!o|CXs= z?D#+|h%@(%(*H~{!qFhEjmV!KZ}=oL#$ILM+`LMysH1>8M*Hr#Azq@>*Ul< zQtyLx>f~1tCmH!%N1XQ9PhNq1RIF1c*Fu~}kk3@aX(9Iu@>Jy073peVEkY9qn73C zeO4Zl{Z^{u`jCd_kFdS@#J=ZBbM?z}*606nWrBQfo`_V*X%)7gse5jfTp!*=wc;7N zbv5|m(4>^Z#S=sd_GP%osnlI155wnJ_g(KR^DuPsuOTu{dgV8YhUp{4{s(P|>wCVf zeg0U*H&bQUb@43y7lGs_^Z|IiK2A9YDbcT$2v_bug;LlR@ zg**;=8rG?k$HAZJ$a}h5U3sUUQ_1fAQeGTv)nn>bT?U^;zb3$+6t_6`!);Wu`}qp0 zA^$$u)Xg|1*Z*hdE2y3Ce%YUv{u;WKkHPHMq0W=lqgU668G+C5ShQk-e6F4-tCGib zk^1xIk23h)V(8>O-Aj@?+^31g*1nl~lPAnLJ;kc7Z4&)<^cLjG(5GOXI(eBDCqw_% z%4eqjJK`*n`(6A6^bfI4o!tFJGoK8d`NSj6O6(`sL_R54r%qmke15>^6?oM}UITrW zjN_`C+!^`&cy;@WZ7%Q9Lym|eO*hXKXZ(g#y(@dZAb$dVH`b|>tNKSR*-+-NUPsQo z_((o~@E+XjAUB3?^*Mc*ls>q4Tbk5=z&dquM&l7brk`J;)5j$Ed>#AA2jKHhSf@_D za<}xk8}hz<%FUngeu1V-Ew|(awX6C z(Wzs8^rP^H=e?X;zw!cJ*2wR<_p6Kio?z;hj|R6#zy7zKFWKk+|MmZh1Ff!E?^E%M zJQwmkm=L6rS3x)Hu{XuPdpC4)H~Bslh@;~^`MygC(#Sm?Ju&C|z{S?O5l7Ckp(Bnq zro5@==q1?4`%GPNi-UZuJpYY2%;%Q6kPqWn^%$Sys_WsB_i|Iu=ga%7B>vtncgmkT x)i3z`xBah$*@mNXVgCY~>!SJ(e1BDNFsCI6nwhr{?!ZXww6@O09 z0f%1H%ET!QKPc-~N*oo8Q&m-_A9}5D6^c{msvJK!^kj+COV$lQ9A8{0UftqAH$OR! zql8NwXR*Y&M%E3MIKf}zj-7YlDnInK!ZYiaPU>8y zxLH3e^MgYp&eD#T*LbPx)3#X8+($*~yU^ftT_sxJ(1yR+^<~r-mg}tcIx*`UQe0=v zDUQ~tGxA@+b>cYNb~w=BfK4E0Hn`qY;CBvjT1{Xl<)B@UUD`d-)%4Zif1t3T>CXQ>|?dh%J$SHI_zT+|(X zF#Alu`eNFL;n2tnb;kTsUcvO(&yBov&P#Jzr#0#vyFMwDw?20sIQgmS>>o+tcv>4C zYkuO+i(k0A<(|{}H7+<}{UfP+{ov51-?5$>=DwjmzxZR-CzSIm&5hx-PHWT${Xlhlp=*t0z{!1Mb=Z)}gS@#X%s3iT>*?|rV zucgdazGJB2&;upTdRezu;#61WtB(V{Rd^j`zHab?L${MSSa%q49;~jvI{2aE)^Xl> zZOzw6$56weM@gK0>)JO7Y;rxZ{M^WDU_pb^9Og9TIn*3_%2QwOj66|V`Fua=Ia6|p z3w_dg3DsqP%_)x7axUziCx58~)US#i%bV@`_z&g?uu&PFkbRy`(;3 zE8{?epPUZjDAhyigMOetwBDBxyI&IPdfHKT+rX^4t*1 z^$9h1p}{FaavQF%=FphunDt5O7<+!^b!30dDUR00KYD)9AM6}DcE7MH!$MwI2S2QU z2FLYNE&3j#Uy};uCDZXgSd`I8vXjyirYq24E-{&dpX(g)aeagSeyR5AG(I{ z>MbT+Sie)D0}eemwdebBAN~}$VwLfC!|7}TMh2NSX#aOevt zPHDKYtg9_?va9xs>mR(^51k}DR5dRyM-@2W(1=q@*2#SXh4HCz-GieY=*xxgJ+C5e zLS%vi4t=id=leG~uV%8J3C>r0RvdlEfvzHPx_wkS`X4C+9dPI>5~uf1cdwZ`)b+&j zJkR;gfd)7Ij@T`D{xpZa;k+)-t@@zU+so&|7V8MD|yBJ||kJGxC@GB#u(h$ou5_ zuu@U=K|jzRUa`l42IuP~{Qj=?h%Ir@FXWGU@cN9{^-*}gIAO_O*9Uo_KO86X41Gtx zGNn!@`qfkFELQ3hw4lK`erxWR!uwBacYw?@>WI1`FKBSS-d6U3yxLkbWDcb+C(n{AtUSjOTxf8LqYlY)S$ZK?{kKZXqx91O4IkOB`r$x*qyk>b!)l*67#WQXjEgAC+oB zgLC{;*+;76m+BJN^@*8RD77*#aNbWt_SYObAbsB|brvgiwp?g%@+;5U)-`3G zaX+c2#8Cy;0S&|D4&Rh7A`=mKE z`i1j@^Ta&<4ChOJ@1x=+P6Lbb*XIlKiut8DzDq4=aJpWW!ueoqL)*Bl!2gSkh4h%3eWeXqa<=j+KFM{{VLi%xQW#7f@-E;P9K@jg^QY;tAa z(GT>8))&{o=dsG?CC8D!zt9{S{mMxHVqc5Q6g!Qrm*smXzDCwOKXhZ^4<7vdw&#vL z<5Rgy~o6{LtmXC;$1*mD>xR@qn z8x3i?J#%P**O;x$;SoPHIQgk7W&V!%p_>@L9!2&Q5=Sh@$!ldnrwJ#&aq@GV=FnG4 zoc&`)?>!h_WQ*mzp1j6|1}DEmz1d%LXyQqWyez4+Sl&-&7Z)0w{MO0)sOHeQQs;Ep zw^*ron7c3 zg_B>XJ^O18OT*AJm-VY;({r)DiQ-e6JrGoUTXk_osXLd})pO=`QmumhYS0!G-Q7oZ=+N z=S$6@bEM8PCpPBAaNqe|E9VE?To2umN&bdI?=|(+^+CUohv7>11Nx-AR5RWO`lU5- z-7l8mQG#XP z%hp@%>R0Bce(34KpMJ``7W=`Wtw}3#e>mF5UUA68N%2b9+J7~#sc-qAUl$%;tzR=E z4mfl@iqqU%CF|ajINQt2Yl|QHP2t_>^_2q-y+`6KmUVAPoJ;BTk^}vk@VQ0iwa5X7 zZYpt>$-3ndr>l95eaC^8{JaeFTDi#&4*jjfd1!Fw%j=gr-UPO0p7Wkv9dx0=_4Tq> zTQ#-7q1$X2c;}BhOTAlp?rpJr-m~jl(BRrnsRwFvzcq)J=e)gczx~0xYYM!zZ25cz z(_QHIg_B=~{QgmM==FwoD5{LJAu3Mr^&AUYEbTM2;b!tT+{7`uRN_3kZtLo$tDmr^ zv8DA9b-e{WRXE-MvfdQOaOlZ~AA8XYNu9-VUK{$j(BO3co5Wvp=uT4SwYZC$XmIjtlF$B{Llf8aiP^7` zi(~gIIf?x>r#M=pUmwYQ$T{P_th~~NmUG7YJdnozno}IDalSCGm|x=ayI9cRbp7s5 z?5{a=7pXJmM=aep;nhJr7ptG}z>VLQ*9K7>!=Y)ubbavL*cg2-3a_)E(NDU6gRUH;OC z4sv^3a<$Ncr<*?VFw&~fW*Om z>dAiS_ptB&%Maa9_Cx#Vb=U!iE|NI8a=&cZ59PgZzXN@p@Wl8IiSco3{ov38CC-M0 z+tsH(ukys|a}lcc%-DPj8eI3sOKCGy_lY_5jHWXt)T;WimoLv9u{1Nj`TIrLG9GkoqN#WUB=vBgqenXN}z(BQ@|)Mga<8xB2${9p7wl{#-?%Xxis zvkNVj;;$S%l=IQt#4&oF)VbO$uFv^wIZi?+7g{XEQHfV^9L-G}qpOo&k&S+!Kg19F zV)rZW8?LkFv`%aEYrNFwKDK<{UHLBb{ldxb6L}xh9J)a2gMNvnyn3M$>(4dN;Knal^H%aV9D2Ug`9Z1A1h!P4k#CH! zpv6-B!qSyzNQ3YpRb%&v=5~@ z#nHOPSP#X69@=ex!+$C{c_`6aOlY&n*EINlXAb?? zP2pVwaz6Hs%JW?;&9e%==|F><{)BSn=M>GMXYskP56)Tf?ykc#yhChh{v#LIE;Kmx zd3u26&~WJSyO{7YHs2vt>;Q!=m+{^aBy2o3mV*9_f6w-Za6gh^@7yr`KbDYcJ;NO!Hu7qoJ)Bb z4n0HaGe`1zIy$cz$rd#FY5YPZ^6v#ThaM++#mPL2r8-}nHr9d$*M7kn!|tK|7!Hkj zE|C39iQdoQ4i>anuCv#L;wa6{ew5ajpZBFcYoqF;>fdHTi=};r9v#eaG&f<4_M|@O z2l`{Us-5UUgVXh^=D&=QC^ghX({9HK?nnOPzbv`2Vd@yRBBU47Y z(9!ceuJahm$8ZzJ=sRSdF&|>Np0{UO&@xA~&$6BJbDZWTj?s@xeQ++ONB66}Ps)(ZwrwT{JgorQN2Ei8I1$$JRFQZ?)d6 z=|YR;_rZ7fIpFj@sP#mO^Q8BK#5o{w%zb=^n!3>C!p%Oe{;)CC#c=3J5~qi(%Z`c@ zdcBql-BY-^|IZi7=P}Kp$4i{|W!<`{IKkO4tVy zJ+Y>q*7UTIE_5T|`g&Q?;yX=V#GwzJKl18P7k}(MEAOLXX`VxW>g+<#uY8}?b-S|- z)yHrhhx7tI7k2A^wHaK#Zl0IPmiikhYH2~Y5pMj{b^QMg8V=oW!-)5`JoS`)iL6VC ziW3@flM8)mRGh@txA1;6H*t*KBXM5+=Q}e8yfxW=g)QZE=t!;$Etd8fTGor=7;fSi zU1!ACYnyI=#Ty}Y6l>N8vpNm2pux@c!D@Y_3FXtAEY9z3>#etshkoh$Pkk@|R1^~qS#*@Xtz zern4(gQ+ftLuW~S&@br=%}-wAzAkihzd}22pg4w`I7aW1e&PJ!JQ@DqJYRJG|NrCv zOa*qyc@S&nU)^z40~Z?HT)!Z-A$0JBw|_+SFuwok3;RYrOuYY~K4;7; zI)1FYPU+Vfx^-JZ+(6%FK|=Cvz@R1 E1xhBregFUf literal 0 HcmV?d00001 diff --git a/explo2/output.log b/explo2/output.log new file mode 100644 index 0000000..9fe6bc4 --- /dev/null +++ b/explo2/output.log @@ -0,0 +1,359 @@ +-------------------------------------------------------------------------- +[[10849,1],0]: A high-performance Open MPI point-to-point messaging module +was unable to find any relevant network interfaces: + +Module: OpenFabrics (openib) + Host: doubs + +Another transport will be used instead, although this may result in +lower performance. + +NOTE: You can disable this warning by setting the MCA parameter +btl_base_warn_component_unused to 0. +-------------------------------------------------------------------------- +contraste 100 + +échantillon 0 + +fraction volumique 0.95 + +itération 1 + +itération 2 + +itération 3 + +[[ 1.20992624e-02 3.02791793e-03 -1.02529145e-06] + [ 3.02788658e-03 1.20077618e-02 -1.16644965e-06] + [ -1.18473687e-06 -1.24265480e-06 8.93749155e-03]] +Solving linear variational problem. +[[ 1.70424052e-02 4.89825441e-03 -7.69472093e-05] + [ 4.89825441e-03 1.61316669e-02 6.91756718e-05] + [ -7.69472093e-05 6.91756718e-05 1.16113030e-02]] [[ 1.20775752e-02 3.04137065e-03 7.37730654e-06] + [ 3.04137065e-03 1.20558624e-02 4.92922601e-06] + [ 7.37730654e-06 4.92922601e-06 8.96177485e-03]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 1.22108038e-02 3.05125606e-03 1.32488415e-05] + [ 3.05125606e-03 1.21846486e-02 1.73639432e-05] + [ 1.32488415e-05 1.73639432e-05 9.07862742e-03]] +fraction volumique 0.9 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +[[ 1.32740169e-02 3.33357254e-03 1.43881667e-05] + [ 3.33334384e-03 1.31571064e-02 1.64203022e-05] + [ 1.49196341e-05 1.80195471e-05 9.65563827e-03]] +Solving linear variational problem. +[[ 0.02615363 0.00834348 -0.00019011] + [ 0.00834348 0.02625687 -0.00022318] + [-0.00019011 -0.00022318 0.01685809]] [[ 1.32744971e-02 3.35686441e-03 -1.50519521e-05] + [ 3.35686441e-03 1.32579785e-02 -6.05919346e-05] + [ -1.50519521e-05 -6.05919346e-05 9.67913656e-03]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 1.36476883e-02 3.40746336e-03 -6.49054630e-05] + [ 3.40746336e-03 1.36438581e-02 -1.36349536e-04] + [ -6.49054630e-05 -1.36349536e-04 1.00430480e-02]] +fraction volumique 0.8 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +[[ 1.62232788e-02 4.13371351e-03 4.14889838e-05] + [ 4.13370762e-03 1.64868325e-02 4.66268077e-05] + [ 4.37297149e-05 4.76379900e-05 1.14602014e-02]] +Solving linear variational problem. +[[ 0.05332394 0.01859865 -0.00012353] + [ 0.01859865 0.06727432 -0.00025963] + [-0.00012353 -0.00025963 0.0331858 ]] [[ 1.60117130e-02 4.14049795e-03 -7.66289361e-05] + [ 4.14049795e-03 1.63609283e-02 -1.07075702e-04] + [ -7.66289361e-05 -1.07075702e-04 1.14297353e-02]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.01729981 0.00444632 -0.00021 ] + [ 0.00444632 0.01799472 -0.00016343] + [-0.00021 -0.00016343 0.01300395]] +fraction volumique 0.67 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +[[ 2.26335969e-02 5.95560161e-03 1.00373253e-04] + [ 5.96073989e-03 2.41277834e-02 9.72608542e-05] + [ 9.35996553e-05 8.39156181e-05 1.54202215e-02]] +Solving linear variational problem. +[[ 0.11083471 0.04055113 -0.00243946] + [ 0.04055113 0.15663066 -0.00354049] + [-0.00243946 -0.00354049 0.07158718]] [[ 0.02077976 0.00540522 -0.00033962] + [ 0.00540522 0.02240423 -0.00038031] + [-0.00033962 -0.00038031 0.01435685]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.02519844 0.00673205 -0.00151563] + [ 0.00673205 0.02923942 -0.00148107] + [-0.00151563 -0.00148107 0.01978838]] +fraction volumique 0.6 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +[[ 0.0294117 0.00767509 0.00014536] + [ 0.0076841 0.03093682 0.00013157] + [ 0.00013884 0.00011608 0.0189558 ]] +Solving linear variational problem. +[[ 0.18423775 0.0666173 -0.00374027] + [ 0.0666173 0.2350969 -0.0017721 ] + [-0.00374027 -0.0017721 0.10732952]] [[ 0.02463847 0.0063728 -0.00059286] + [ 0.0063728 0.02587365 -0.00046276] + [-0.00059286 -0.00046276 0.01622806]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.03480306 0.01012984 -0.0014995 ] + [ 0.01012984 0.03969563 -0.00133254] + [-0.0014995 -0.00133254 0.02854458]] +fraction volumique 0.5 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +itération 8 + +[[ 0.04771381 0.01283749 0.00015452] + [ 0.01286646 0.05090793 0.00020884] + [ 0.0001266 0.00020566 0.02918238]] +Solving linear variational problem. +[[ 0.31695038 0.11489078 -0.01177378] + [ 0.11489078 0.3677751 -0.00220521] + [-0.01177378 -0.00220521 0.18297196]] [[ 0.03176117 0.00778496 -0.00099927] + [ 0.00778496 0.03110973 -0.00050345] + [-0.00099927 -0.00050345 0.01857851]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.06250104 0.01698269 -0.00673719] + [ 0.01698269 0.06457141 -0.00473635] + [-0.00673719 -0.00473635 0.05033523]] +fraction volumique 0.4 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +itération 8 + +itération 9 + +[[ 9.70392655e-02 2.57731214e-02 1.03421106e-04] + [ 2.58248034e-02 1.04844200e-01 1.34385696e-04] + [ -1.09949766e-04 -1.20666308e-04 5.64342259e-02]] +Solving linear variational problem. +[[ 0.49922347 0.1728598 -0.00267058] + [ 0.1728598 0.54624358 0.01748384] + [-0.00267058 0.01748384 0.28303892]] [[ 0.03890438 0.00880798 -0.00062893] + [ 0.00880798 0.03865878 -0.00013273] + [-0.00062893 -0.00013273 0.01982406]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.12201314 0.02845497 -0.00316679] + [ 0.02845497 0.12434439 0.0040708 ] + [-0.00316679 0.0040708 0.08763849]] +fraction volumique 0.3 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +[[ 2.58191999e-01 6.55998524e-02 5.46356287e-04] + [ 6.56098766e-02 2.64548486e-01 -2.21261072e-04] + [ -9.29187706e-05 -2.39399410e-04 1.46156530e-01]] +Solving linear variational problem. +[[ 0.5781815 0.17155351 0.00843673] + [ 0.17155351 0.61044493 0.00685478] + [ 0.00843673 0.00685478 0.34142827]] [[ 5.21567469e-02 1.22840092e-02 -8.10121913e-05] + [ 1.22840092e-02 6.14301381e-02 -1.00467824e-03] + [ -8.10121913e-05 -1.00467824e-03 2.46848320e-02]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.23234596 0.05480467 0.00374799] + [ 0.05480467 0.28968937 0.00177932] + [ 0.00374799 0.00177932 0.1413434 ]] +fraction volumique 0.2 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +[[ 5.06687771e-01 1.31044630e-01 -9.57354472e-04] + [ 1.31078807e-01 5.12980697e-01 6.33871025e-04] + [ -1.48204770e-03 8.52852166e-05 3.18945588e-01]] +Solving linear variational problem. +[[ 0.61094853 0.17580768 -0.00738516] + [ 0.17580768 0.70086005 -0.0070871 ] + [-0.00738516 -0.0070871 0.42429732]] [[ 0.0948014 0.02283239 0.00012848] + [ 0.02283239 0.10807511 0.00071642] + [ 0.00012848 0.00071642 0.04846542]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.45237788 0.11777589 -0.01071644] + [ 0.11777589 0.53993993 -0.00342461] + [-0.01071644 -0.00342461 0.31072088]] +fraction volumique 0.1 + +itération 1 + +itération 2 + +itération 3 + +[[ 0.83209966 0.20619228 -0.00191994] + [ 0.20606532 0.80117647 -0.00188859] + [-0.00223212 -0.00199422 0.55322644]] +Solving linear variational problem. +[[ 0.84887624 0.21506862 0.0079913 ] + [ 0.21506862 0.83172375 0.00272714] + [ 0.0079913 0.00272714 0.58063397]] [[ 0.28188858 0.04491097 -0.00269438] + [ 0.04491097 0.23804175 -0.00043351] + [-0.00269438 -0.00043351 0.11891348]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.80122863 0.19298084 0.01075216] + [ 0.19298084 0.77530804 0.00661389] + [ 0.01075216 0.00661389 0.54227169]] +fraction volumique 0.05 + +itération 1 + +itération 2 + +[[ 0.96463944 0.24212291 -0.00227569] + [ 0.24209127 0.96617804 -0.00196716] + [-0.00234799 -0.0018619 0.68913631]] +Solving linear variational problem. +[[ 9.60473491e-01 2.44031943e-01 1.51905255e-03] + [ 2.44031943e-01 9.81025000e-01 -9.49666819e-04] + [ 1.51905255e-03 -9.49666819e-04 6.92409569e-01]] [[ 0.46374746 0.08506615 -0.0175395 ] + [ 0.08506615 0.46521117 -0.00618941] + [-0.0175395 -0.00618941 0.2447475 ]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 9.36323397e-01 2.33307506e-01 2.30107497e-03] + [ 2.33307506e-01 9.57055165e-01 5.50400438e-05] + [ 2.30107497e-03 5.50400438e-05 6.76286135e-01]] +contraste 0.01 + +échantillon 0 + +fraction volumique 0.95 + +itération 1 + +itération 2 + +[[ 95.72300274 24.21735067 -0.19385434] + [ 24.2191621 95.82953228 -0.17668746] + [ -0.18689262 -0.17564719 68.93125031]] +Solving linear variational problem. +[[ 95.97885033 24.42943992 0.47491726] + [ 24.42943992 95.52703255 0.42717941] + [ 0.47491726 0.42717941 69.48659933]] [[ 47.31987511 9.78605674 -1.41510456] + [ 9.78605674 47.96950138 -1.08631823] + [ -1.41510456 -1.08631823 23.8513985 ]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 93.65780098 23.53595537 0.63523796] + [ 23.53595537 93.3306977 0.58826397] + [ 0.63523796 0.58826397 67.96223466]] +fraction volumique 0.9 + +itération 1 + +itération 2 + +itération 3 + +[[ 81.47094469 20.39878704 -0.24766306] + [ 20.39319326 79.77704287 -0.34479898] + [ -0.24795965 -0.33423204 54.87336945]] +Solving linear variational problem. diff --git a/explohh/.vi.swp b/explohh/.vi.swp new file mode 100644 index 0000000000000000000000000000000000000000..5cc85bfb757c71a9396852008b6ba1ca295b63e6 GIT binary patch literal 12288 zcmeI%y$*sf5Ww-n-Dva$PFflNYboe(ZZMo)p z`E5-DcfZ})o24hZjmyu!y+aX|sB;mg^{>`3%&x(uHQEFoRYK=*&rewtIk&h8GD84? zo&|>f+FQBhA&dFUo=*Ip&m%no2q1s}0tg_000Ib@Bv9m*Jas{Xswk^vG^u84K>z^+ z5I_I{1Q0*~0R#|00D*1-O26~JMj{{G%lH5K{{LHdOT}&k5I_I{1Q0*~0R#|0009IL JFj3%`$`=K}7?S`1 literal 0 HcmV?d00001 diff --git a/explohh/Aloc.pdf b/explohh/Aloc.pdf new file mode 100644 index 0000000000000000000000000000000000000000..0b00375a6155a2aebeb4cee76cb66f1606a77a8c GIT binary patch literal 10233 zcmb_i2|SeB`)_mA*pg5vQ$n)LZf0ajw(JZc6oaudjIkti)4fI6#=aDlbSbGuDU{_# zX3#>GMy0Z~kwPKK|Ga}z-QWHD{Xc&`&*!{z-gD0LJ4lJL&1@7s*5*FT^)|FAbU`P;Al`of+MW=2T`}_a0kFgS09W|38TXihJJwO z&bd-&t~7-2g(Ivy$!=c$9&pv!?sn!b-V|3l9M2v6I~)`eRLoC6+uxr`4+1Pee;ojY z;tzfiX9qB#QG)~E&@^?d;0Rqxh?gscWB^!$U)D6Te^3CKM)7yuKR@%?_J9$>l;Y+^ z-boDu<3h}^sz?ZbI?*pnJkm(CUz&jzM%9h3p3lI}x9zL6t=D9|@g?CEWd6q$n}v z&N{>RkRYE^IxXzLR*Asl*d3$e#!ru$l?_V?@Wx)PGdR01YXz^~`xO#`Be!ZpCpPsi zUtQdPfMj2s&)nyl#2H0mI!=olzvP$~Y}_5B=&+CS)SjWbt1Eg{A=`LmuMD&Kw}#Qa z2V${@!$xukU!SR}*@OJFF_w^caQwYMlul4ntXc9Z(X6u%v7u83t3##`Z$$_9$vSs= zZR+p1e{V^p3^s4&yRS9+e$Od!#RX;x;pb}%zl*xrmQ_C*nQ@YlL_ zjxtsVJT};LQO@tR>BTLUw~cHEJnZ~)MT70F zF6$czDDK*~?|i_KGm|$TA-(rrX68rBUTJx}z9VAa6AAihxA+8!^o+`WW3l@$DXdLt zU91FHdJFA(_z_%QM&QQBBuS3tZJLDX4Yw?w4T8f*T1+G~Owlfv2( z9ZO@+yMI*VimpH0En3)Jq8nKE9PO?BPG(zUud;PZPVlMKN%DqDyCVByJpDGsb{F^B zB}TR0-V$`=ex$kj{|jCRp{hQ{Zk>Z0wNQf?B@K38o7O-%+j?!4==q#X3vrhh)1E7C ziqR396&rb2+22Ql)zup%sy|E=-mzh4y`P>LwoDcw9~vOUyL=GOajBEs%MJfy6Xj5IY3m&ofo zX7VVC6I1N_qtkq%TOGIkgG?9NyMNoO1mq9m(~GMVs6npOJ#fx>MUcekY`Odf-`lR@igZY>Joq`Vygj-$sIL^DxJR~Gg z<8IUXv46m|SQ~=D5PkrLBN0oU?%;nL3EK<*J>>3*mpnyL@17k=b%NmF+6vXIMi1uF>>6_+78|_BxM*77jMD%J>*c=E`DOESY;vN^Pbeq zRrI0bo8O(OS!39I!ey`8XCY0xQem}fZ_}j&qtCk(eN4JS^9{Y4H>Z?j=uz?1NR1mi zMpYUQ?z;=uw(8jG74hKwN!P8}>nVSRN@eq78>Egli>M#PtVtb<+#2O%ELX9Wp(QAA z`k>Rs{J0apQI`g_>9B6UJe092F_XunE z)fdBA>uN)5Pd{r6mJ}nd7Twjl&f2ta@6C<3ctfPzmyMpqDIH$_cDa@-3?CjCW^#u! z6ucg{;y|a5xw{2gQ9O`c$72=zh?RBQ{9GL0kmG#;R=IEoW!>4Fv^#Ns$e3;-?K=HN zt_Y)@v?Bk4Wq!O=ebnnFnUmz*GT!`Cj>eBWgGBC0V7H51P$w1`*q&%&*YhdzPT#y+ zqh~E~L4Ea_fa=FO(W~BYGw@df!fdSzcJ+ zAodReTo{!6s{n!bU+5fh7EJ;=sFjr-FtNh5%MQ8a+T04B@$q}Goh&@F(~gz8a(@x~ zp2BNmO*IptGL`BpoB1vV%MLbX!*dmfPiHZ&v^?r@4v>)B_q!P0`IT3^1dje9`kYGI z&O`OG@3SI;54@@}UmxbZo))q~ScmPr%DY0$)b6T!UXt0;o|&lI@09yaCyzGfN={0i z7<#|G)c=LNcgW@znkc%h6x**elFfg&9cqd35U1|q^KG%6fH^#^eWJrq3AJ!upQ9EX7; z+`U35P>eGh3w|303X1f-XhC#{6~KdlKB6^yNy8{HE~Vm2lPz~@~h#D89n+Ci>R2uDN{p@27F7sboNlMc-b z0t~%TMq^aD5t>I39D|-ky@>h#3>Ad8Ae#1a)%N%BrNEIOvLuIb0RntjG>Am$6h9k? z@Anxj5E;Vm8?XiQfF=gb!vU^}ME-ve?Eg9r21S5l!TN#t7)}6AM-`(ABjS*7Jn)QY zJPwW_BDsVEIHC7YdyFb-w#49}_Wwq(`S%#m79wbjdBR~(a1;^}009dK3`8Wr@p$kX zgNEY(b38;iG#pLD!2rQkajI|}4h#bH6c085o&W$K!824s?=gTtJcS4Sp)rXBA}D9; zpprWo3?2rR5C~{e5UB!9gZkpY#An9h*0D2ow?tk1HP8!4 zCl(M_A0VBWC$6r{=>((?3kkx-{SC%|bYq^NnX4jkU=0CQNC*Cbhy<`EkX`^IzL4g1 zV!nig0@fYIEurTH#9iyz_fVT}1nC8|s!+N2n~{U;|Gbgo+H&q`0Gvu8V-F z5Uyl63hIc^M$b1x{pV&2SS)aY-}x?>?DxurryoNAb~d^3{XMs0UPL>+=s^|~nhw=3 z00;=!ec&8*sJ>L1RRGzQ>vZOl1~A|mpo2O%`}+mWLWi0G=fLd{><9cX25c>m^Z>a7 zSFXe3Mk>H<{=;d1&aQKN%$@wvb6J>gzKxqGfz8QqJ{f|T^&U`1Fr{xff}y8PiEiV2 z?tE}+CkCHCU-{QxF9$w6X{dX+0$VG%S-WxBl^u3sX3NAjMW0FvxV}T|GX0M}X9dZ+ zWvgOzV{Ypy?NX}UJb1?0b#uV)b+vUE6ZLzJ>f>4Oa1&uy8Ch5(4siZM+Jof8=rK9=Y{A~Ox#^@dK3HB z!Nxn65k4xp&lI=8rqgT_ozif64iWtiFG({@-QT#_m^-GMhI}j=j2!t?w3xz{$0!b`+8qJ?iHh0EPhEb~FtF-YF$j!Ygub5UMKOf9&|IeRRPg(|cc zJ)8PzLl;Y2*PT>CRbfj^Dh{hvsN`I2zxHH(=*P6e1CN#Ym+xUX1+2@G%zjnAsc2hH z#GwedmC?yE&ox<@y?)<${tKjHC5N73|34_p^()&mb}MXlpv;t`+Ve4#{fzuQ_n8Dz#Jf2yun_!um(E3zp6U=J&`rpPvtE%BOH zx;zPCS8{|IdsCS)b?5Pq^Zyqt&7}`y&?44lb-c>y;Wu5(JrteZzdAI^WQ*rbdsaDj zv|0&Ok8~o~*;Kz{+KGy{Pcqr=#!8nwai0h9m=D7Suf+MA9SkqB)K=~md91L>d$W1e z?wHteS3VymySHgZ#wYi{yN8p}HE^@f%?4H3h>4ev6C34h^?S0Uu;Ctp)PG*Ep8+RO zR*$EIpQ~nQ4wYm$oRtg6uE>qmf4G*Ae6g;qDCY2+p_;U-e7C#OSMGnOxzlWwZb-GS z`T2<2nWr2?PQ>WIQOeyv<=bK>YEni9$Lp;d?|Gl9B)l(Ml=47We5RjNPnRh%mF$wv z9a6|SZ7&zl|BmU)j0z8h)@Ne9bhyPW5Aq1mJV=*TxGjP2{RhirK_knHE z(nRPwdVQ4VVa@rYl{t}ntnfW9>ES)6=tf>e5-<6V)qLi=PuDx;uF3mF`M8CWkRPWb zlv1P^RU31$W<^a`eOiR}cxEo)b>X38W`VrfB5^>wEREsCV2%Cpg0!skx2C;UHs@Qt zlJog<-AfKZ^3jB_UA?8IQb|a8D@OY%?TnToqxyM7R1aM%epL6>Zq-9sr~FSJtWEPy zLJ-=XC7FI5FVdNe9zXao+LuPTfpNOrZp~b^NUi99V&>Nr^dPNkOA_ggcu(Bu*Z%EH zcA`{OCSsyB$ebm2Q9rNiYeDslebw`^=cfRJvw~LUW8^o1GU){-%Vuh@hcx69+c7bBIlz3AF+r;B(w zy7xU(z5NSC)+x2xJ{BB2zI80pL8WeIqQAqX@FWpg#+17-r`{5;SLIAtrA|$tS{Jpe zoUBnvU=b3rnRK_iLxBjVfHxh7ne4kBHX}Tp&+Jo8Z9cmUs*K6YDK2@wVbg!McMB2P z1H*h(XmkKj!yxebH2P+0zVDRfO_OzE2ySQcOcWwe7GWuVWpHXlcL*|mq-vo!p8`YcwXwQI?=}pNP;khuwRL&`m7k7(q3zX`cAV;EhqNrqyCNAC2Bd+>ph<@&NxHv#Ky*N`q`-Q>8w**k_}+ zdAG=|s5t&nmD1uD%sfDdRT^6@f4u+m@wJ&tVx-Qgw3SaWH`wiYZkKN<8csWN!P`bH zyt4klljLBsL8i8$~j+#GV&V5Z@%oT`|E6KZrJPgh}!b- zJYQOAHv_Tmt`W!0kGa*Siuv_WMD3^hZ}qlLrPrnI_>j5}9`x~k9X_0X-(k_a7HpSN znp1&C^l$MKh%b1FTY@faF?*Jm7t8R-1%?-N*Kg|dOjW++0o;w|}w^_cI=sbXoJ4=028OtNTN~wwMEod?cg%qFl-3gDzV8 zCj>ZseLV#cwYGxlayfB%;?G{O?@y=r{=GQ=xaY0arHD7(mBoo4UrbZ;h;k|g`JJ9O z$6o>I&~T}yeYpbVs#|zIrL6GTvcE=aL?NY!FcfvYL!v^foEr`S6(0DL zVku0hJUCW&EuwZg-)E!qW$SuZYW#Mn!O*`nHsZl&&rMDN?*^kdy3>JY$EV-(aQ+*+ zQ{MhcI}teC$>_wPX>(2=r}_(h>`ZsJ{MYws5w)S$h1hb}-_%FcKDECTez{q|)I+U9 z{=9xd%Gp$JoCx;Bs{#>5yke7gp9gYnjv~pZ{qW<-Xfp$*oaxUW0 zI8npreTQA~X-@7NVfk`w1M`t=0X8dEh#5Dkt8ezE=8}@A-z`=VDPM1>s!K=YfuLh7 z78a+?_EEO4i2c)1JdN~5%w>nd2w$U6vQOC>;6MidAr4$<7opevpLFv=o}I71^4T&! z5oudM*b_0(`7w=OC9he0@LlZQ)d~e9Wpa$v(G1>EtLMH|>Y37BYdEFCllsYyACs0f z+?ILSpxkHcf9`WG{~ed#A$x6%@g+|&#cj0BeR zX&kwBWAbfA2urR~h|0>&&0du(-%IAr;3u_dh@C&&uAbfs!e;NT7)9pOL|6D{cssO;p;KYFhHO z{=f)NsC2SMYwo@D;&Rl!Msbtn$3q8Xn5O+vh9ymuhv z-9DzHA!lSzzNqU=;uA*(smsA{fKT)C*g4;kFY=nn2v32r-1x-m_rZIDyaSL<0imra zruy-ACYHQr%%&l+&azOJ+)HzXWJA%CK99jNY2yGER4bBl`?Ux_F_*&sRf@!)mz?VDO zOl8ZNIu(lYt&BWYN;-7s12KIi{dp}O!;CtXCyTUA0Uo1Z0IOjnk`F`#1 zzU#+G>Ob5kcq+BkZJTDHpweAVNxd%nQ(>~Hm`&YC1}3lV@WyGyaDUGtMA@j4jVtG@ zgUJZb?iQ9Dt=mY@p&<2Z?*NFW@p)~dnx(n91>plDy_3-fA&#l0t-GE)4nLQNZz3VP z7a2NAM@|+U7)D?3ea|oz94$|Nuu~{F>>#1Vpt}FXwB<#4q(rc+aB0mg1T4_vRwsA4 zZjsv7`rW2wvNYD?d#+aY-#J%$S@3H1f3h1SGTz!O)-e?Hw_R3uj&i!<=*NK#rdr4kT&VG1vEQbm{Y{fxr^-+o zt~J;7E%O}ductoOHp>*ydym?afBseo6PGp#V{bH_} zj>r`KEdFKiWwnw|6^@XLnCQs@K3{0&Wah!i4_X?sImPYAC!^!9#fr7!f04%LJlXJHgVt91nj znmk%tw~X~Q!Ip#>EW&3)j&g>Vapr2;ulHlzycv|w!xFHV`HUo#OPkN}=}?0~Vh^Ml zjlJBU3o!2YL!|GQT(jRmeNPaA?Dm4r;3FKm7u8=E`s8V=t{NJNMj{DF0t$`9pwWs* zv^)|iuK+{nK<&US*H%9=J;0Yr_w{mtqg7BUC?&WjogScuKzx6pLZx|tF5t$_E!Y)& z2(++IfSWtqh3x7BGOM%w6&55x!)O$D7`R%&!a!nv?ho8+;6N7G9saEjx&ws%;QsS< zMCgL%r#dW36$t6ac32d+zWccjO$74zsT~Ho%=@_xp!lf{iB$cy4g>C0er|_Gg3R=f z;~|NJU+SE*{ufz+JBeTTqA=)R>#*RG>!IdzVwsTGuMax50;2?sQ>@~ literal 0 HcmV?d00001 diff --git a/explohh/Champs_admissibles.py b/explohh/Champs_admissibles.py new file mode 100644 index 0000000..dbe7a74 --- /dev/null +++ b/explohh/Champs_admissibles.py @@ -0,0 +1,186 @@ + + + +from dolfin import * +import numpy as np +import matplotlib +matplotlib.use('Agg') +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 champs_admissibles(tabC,tA,tab_tau,C0,CHH,res1): + k,ll,mm,nn=tab_tau.shape + tabC0=np.zeros((k,k,3,3)) + for i in range(k): + for j in range(k): + tabC0[i,j]=C0 + mesh = UnitSquareMesh(res1, res1, "left/right") + + 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) + + Cvar=LocalC(tabC,mesh) + C0var=LocalC(tabC0,mesh) + tau_var= LocalC(tab_tau,mesh) + #Avar= LocalC(tA,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] + + class pix0(SubDomain): + def inside(self,x,on_boundary): + return x[0]<1 + + pix=np.array([[pix0 for i in range(k)] for j in range(k)]) + + for i in range(k): + for j in range(k): + class pixij(SubDomain): + def inside(self,x,on_boundary): + return x[0]>=j/k and x[0]<=(j+1)/k and x[1]<=1-i/k and x[1]>=1-(i+1)/k + pix[i,j]=pixij() + + subdomains=MeshFunction("size_t",mesh,2) + for i in range(k): + for j in range(k): + pix[i,j].mark(subdomains,i*k+j) + + dx=Measure('dx',subdomain_data=subdomains) + def reference_problem(C,tau,C0,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 sigma1(v1): + Eps1=as_vector([1.,0.,0.]) + return dot(C0,eps(v1))+dot(tau,Eps1) + + def sigma2(v2): + Eps2=as_vector([0.,1.,0.]) + return dot(C0,eps(v2))+dot(tau,Eps2) + + def sigma3(v3): + Eps3=as_vector([0.,0.,1.]) + return dot(C0,eps(v3))+dot(tau,Eps3) + + def A_CA(e1,e2,e3): + return as_matrix([[e1[0],e2[0],e3[0]],[e1[1],e2[1],e3[1]],[e1[2],e2[2],e3[2]]]) + + def sig(Eps): + M=inv(C-C0) + return dot(C0,dot(M,dot(tau,Eps)))+dot(tau,Eps) + + def En(v,Eps): + return dot(sigma(v,Eps)+dot(C0,Eps),eps(v)+Eps) + + def t(A): + return as_matrix([[A[0,0],A[1,0],A[2,0]],[A[0,1],A[1,1],A[2,1]],[A[0,2],A[1,2],A[2,2]]]) + + def E_CA(A): + return dot(t(A),dot(C,A)) + + 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,Ve,Ve,Re,Re,Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10)) + V = VectorFunctionSpace(mesh,'DG',degree=1,dim=3) + + v_1,v_2,v_3,lamb_1,lamb_2,lamb_3 = TestFunctions(W) + dv1,dv2,dv3, dlamb1,dlamb2,dlamb3 = TrialFunctions(W) + w = Function(W) + + Eps = Constant((0,)*3) + + F = inner(sigma1(dv1), eps(v_1))*dx+inner(sigma2(dv2),eps(v_2))*dx+inner(sigma3(dv3),eps(v_3))*dx + a, L = lhs(F), rhs(F) + a += dot(lamb_1,dv1)*dx + dot(dlamb1,v_1)*dx + dot(lamb_2,dv2)*dx +dot(dlamb2,v_2)*dx + dot(lamb_3,dv3)*dx + dot(dlamb3,v_3)*dx + + + Aloc=np.zeros((k,k,3,3)) + Bloc=np.zeros((k,k,3,3)) + + solve(a == L, w)#, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"}) + (v1,v2,v3, lamb1,lamb2,lamb3) = split(w) + e1=eps(v1)+as_vector([1.,0.,0.]) + e2=eps(v2)+as_vector([0.,1.,0.]) + e3=eps(v3)+as_vector([0.,0.,1.]) + A=A_CA(e1,e2,e3) + for ii in range(k): + print(ii) + for jj in range(k): + for j in range(3): + for iii in range(3): + Aloc[ii,jj,iii,j]=assemble(A[iii,j]*dx(ii*k+jj))*k*k + + for ii in range(k): + for jj in range(k): + Bloc[ii,jj]=np.dot(np.dot(tabC0[ii,jj],Aloc[k-1-ii,jj]-tA[ii,jj])+np.dot(tabC[ii,jj],tA[ii,jj]),np.linalg.inv(CHH)) + plt.figure() + plt.imshow(Aloc[:,:,0,0]) + plt.savefig('Aloc.pdf') + plt.close() + C_CA=np.zeros((3,3)) + S_SA=np.zeros((3,3)) + for ii in range(k): + for jj in range(k): + C_CA+=np.dot(np.transpose(Aloc[k-1-ii,jj]),np.dot(tabC[ii,jj],Aloc[k-1-ii,jj]))/k/k + S_SA+=np.dot(np.transpose(Bloc[ii,jj]),np.dot(np.linalg.inv(tabC[ii,jj]),Bloc[ii,jj]))/k/k + return C_CA,np.linalg.inv(S_SA),Aloc,Bloc + + C_CA,C_SA,A_CA,B_SA = reference_problem(Cvar,tau_var,C0var,k) + + print(C_CA,C_SA) + return C_CA,C_SA,A_CA,B_SA + + + + + + diff --git a/explohh/Homog_hierar.py b/explohh/Homog_hierar.py new file mode 100644 index 0000000..55d6750 --- /dev/null +++ b/explohh/Homog_hierar.py @@ -0,0 +1,204 @@ + + +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,scheme,schemebis): + + 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)) + + + 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+1,j+1,l,k] + gamma12[n,3*i+l,3*j+k]=t12[9-n,i+1,j+1,l,k] + + + + + + def auxiliary_problem(C_list, n,gamma,K,scheme): + + 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.]]) + Chom_DEF=np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]]) + Chom_CONT=np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]]) + 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)]) + B=np.zeros((4,3,3)) + 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_)) + for i in range(4): + B[i]=np.dot(np.dot(C_list[i],E[i]),np.linalg.inv(Chom_)) + + if scheme=='CONT' or scheme=='DEF': + for i in range(4): + Chom_DEF+=0.25*np.dot(np.transpose(E[i]),np.dot(C_list[i],E[i])) + if scheme=='CONT': + Chom_CONT=np.dot(Chom_,np.dot(np.linalg.inv(Chom_DEF),np.transpose(Chom_))) + + if scheme=='HH': + return Chom_,E,B + if scheme=='DEF': + return Chom_DEF,E,B + if scheme=='CONT': + return Chom_CONT,E,B + + + def tC(tC_,tA,tB,n,gamma,K,scheme): + 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,Bloc=auxiliary_problem(C_list, n,gamma,K,scheme) + 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]) + tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj]=np.dot(tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj],Bloc[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]) + tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj+2**(p-n)]=np.dot(tB[i*2**(p-n+1)+ii+2**(p-n),j*2**(p-n+1)+jj+2**(p-n)],Bloc[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]) + tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj]=np.dot(tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj],Bloc[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]) + tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj+2**(p-n)]=np.dot(tB[i*2**(p-n+1)+ii,j*2**(p-n+1)+jj+2**(p-n)],Bloc[3]) + return t1 + + + + + def Chom(tC0,n,K_,scheme): + # 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)]) + tB=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,tB,n-i,gamma,K,scheme) + t=t2 + return(K,t[0,0],tA,tB) + + + 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,tB=Chom(tabC,p,C_hom,scheme) + if schemebis=='HHbis': + C_hom2=np.zeros((3,3)) + for i in range(pp): + for j in range(qq): + C_hom2+=np.dot(tabC[i,j],tA[i,j])/pp/qq + if schemebis=='DEFbis': + C_hom2=np.zeros((3,3)) + for i in range(pp): + for j in range(qq): + C_hom2+=np.dot(np.transpose(tA[i,j]),np.dot(tabC[i,j],tA[i,j]))/pp/qq + if schemebis=='CONTbis': + Scont=np.zeros((3,3)) + for i in range(pp): + for j in range(qq): + Scont+=np.dot(np.transpose(tB[i,j]),np.dot(np.linalg.inv(tabC[i,j]),tB[i,j]))/pp/qq + C_hom2=np.linalg.inv(Scont) + C_hom1=C_hom2 + 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(C_hom1) + CA=np.zeros((3,3)) + CDEF=np.zeros((3,3)) + SCONT=np.zeros((3,3)) + + for i in range(pp): + for j in range(qq): + CA+=np.dot(tabC[i,j],tA[i,j])/pp/qq + CDEF+=np.dot(np.transpose(tA[i,j]),np.dot(tabC[i,j],tA[i,j]))/pp/qq + SCONT+=np.dot(np.transpose(tB[i,j]),np.dot(np.linalg.inv(tabC[i,j]),tB[i,j]))/pp/qq + CCONT=np.linalg.inv(SCONT) + + return CA,CDEF,CCONT,tA + + + + + + + + + + diff --git a/explohh/Reference.py b/explohh/Reference.py new file mode 100644 index 0000000..85f469f --- /dev/null +++ b/explohh/Reference.py @@ -0,0 +1,209 @@ + + + +from dolfin import * +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from matplotlib import cm +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, "left/right") + + 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) + #Avar= LocalC(tA,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) + + #def sigma_HH(Eps): + # return dot(C,dot(A,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,'DG',degree=1,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"]): + 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) + + vf = Function(V, name='Strain-'+case) + vf.assign(local_project(eps(v), V)) + ffile.write_checkpoint(vf, case, 0.) + + + #sEF=sigma(v,Eps) + #sHH=sigma_HH(Eps) + #diff_sig=sEF-sHH + #plt.figure(figsize=(80.,50.)) + #pp=plot(diff_sig[0],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig_11'+case+'.jpg') + #plt.show() + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(sEF[0],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/sigEF_11'+case+'.jpg') + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(sEF[0]-diff_sig[0],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/sigHH_11'+case+'.jpg') + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(diff_sig[1],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig_22'+case+'.jpg') + #plt.show() + #plt.close() + #plt.figure(figsize=(80.,50.)) + #pp=plot(diff_sig[2],cmap='hot') + #co=plt.colorbar(pp,shrink=1.) + #co.ax.tick_params(labelsize=100.) + #plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/diff_sig_12'+case+'.jpg') + #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 + + + + + + diff --git a/explohh/courbes.py b/explohh/courbes.py new file mode 100644 index 0000000..fc5e5c7 --- /dev/null +++ b/explohh/courbes.py @@ -0,0 +1,407 @@ +import numpy as np +import h5py +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import cm +from matplotlib.ticker import LinearLocator, FormatStrFormatter +import os + + + +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.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,0.95] + 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[:,:,2,2],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,0,10) +#verif_eps(100,0,3,10) + + +def courbes(c,i): + os.system('mkdir courbes') + fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") + fr=fi['fractions'] + ChomEF=fi['ChomEF'] + ChomHH=fi['ChomHH'] + C_11=np.ones((13,)) + C_11_EF=np.ones((13,)) + C_22=np.ones((13,)) + C_22_EF=np.ones((13,)) + C_12=np.ones((13,)) + C_12_EF=np.ones((13,)) + C_33=np.ones((13,)) + C_33_EF=np.ones((13,)) + C_11_0=np.ones((13,)) + C_22_0=np.ones((13,)) + C_12_0=np.ones((13,)) + C_33_0=np.ones((13,)) + 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,0.95,1.]): + if c==0.01: + jj=12-j + else: + jj=j + if f==0. or f==1.: + if c==0.01: + fractions[jj]=1-f + else: + fractions[jj]=f + if jj==12: + C_11[jj]=0.01 + C_11_0[jj]=0.01 + C_11_EF[jj]=0.01 + C_22[jj]=0.01 + C_22_0[jj]=0.01 + C_22_EF[jj]=0.01 + C_12[jj]=0.01 + C_12_0[jj]=0.01 + C_12_EF[jj]=0.01 + C_33[jj]=0.01 + C_33_0[jj]=0.01 + C_33_EF[jj]=0.01 + if f!=0. and f!=1.: + if c==0.01: + fractions[jj]=1.-fr[j-1] + else: + fractions[jj]=fr[j-1] + C_11[jj]=ChomHH[j-1,-1,1,0,0]/1.11111*np.sqrt(c)/10 + C_11_0[jj]=ChomHH[j-1,-1,0,0,0]/1.11111*np.sqrt(c)/10 + C_11_EF[jj]=ChomEF[j-1,0,0]/1.11111*np.sqrt(c)/10 + C_22[jj]=ChomHH[j-1,-1,1,1,1]/1.11111*np.sqrt(c)/10 + C_22_0[jj]=ChomHH[j-1,-1,0,1,1]/1.11111*np.sqrt(c)/10 + C_22_EF[jj]=ChomEF[j-1,1,1]/1.11111*np.sqrt(c)/10 + C_12[jj]=ChomHH[j-1,-1,1,0,1]/0.277777*np.sqrt(c)/10 + C_12_0[jj]=ChomHH[j-1,-1,0,0,1]/0.277777*np.sqrt(c)/10 + C_12_EF[jj]=ChomEF[j-1,0,1]/0.277777*np.sqrt(c)/10 + C_33[jj]=ChomHH[j-1,-1,1,2,2]/2/0.416666*np.sqrt(c)/10 + C_33_0[jj]=ChomHH[j-1,-1,0,2,2]/2/0.416666*np.sqrt(c)/10 + C_33_EF[jj]=ChomEF[j-1,2,2]/2/0.416666*np.sqrt(c)/10 + plt.figure() + g1,=plt.plot(fractions,C_11,'+') + g2,=plt.plot(fractions,C_11_0,'+') + g11,=plt.plot(fractions,C_11_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_11 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_11 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_22,'+') + g2,=plt.plot(fractions,C_22_0,'+') + g11,=plt.plot(fractions,C_22_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_22 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_22 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_12,'+') + g2,=plt.plot(fractions,C_12_0,'+') + g11,=plt.plot(fractions,C_12_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_12 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_12 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_33,'+') + g2,=plt.plot(fractions,C_33_0,'+') + g11,=plt.plot(fractions,C_33_EF,'+') + plt.legend([g1,g2,g11],['HH','HH C0','EF']) + plt.title('C_33 contraste '+str(c)+' échantillon '+str(i)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_33 contraste '+str(c)+' échantillon '+str(i)+'.pdf') + #plt.show() + plt.close() + + +def courbe_moyenne(c,nb): + os.system('mkdir courbes') + fichiers=[] + for i in range(nb): + fi = h5py.File("c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5","r") + fichiers.append(fi) + + C_11=np.ones((13,)) + C_11_EF=np.ones((13,)) + C_22=np.ones((13,)) + C_22_EF=np.ones((13,)) + C_12=np.ones((13,)) + C_12_EF=np.ones((13,)) + C_33=np.ones((13,)) + C_33_EF=np.ones((13,)) + C_11_DEF=np.ones((13,)) + C_22_DEF=np.ones((13,)) + C_12_DEF=np.ones((13,)) + C_33_DEF=np.ones((13,)) + C_11_CONT=np.ones((13,)) + C_22_CONT=np.ones((13,)) + C_12_CONT=np.ones((13,)) + C_33_CONT=np.ones((13,)) + fractions=np.ones((13,)) + for i in range(nb): + ChomHH=fichiers[i]['ChomHH'] + ChomCA=fichiers[i]['C_CA'] + ChomSA=fichiers[i]['C_SA'] + 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,0.95,1.]): + if c==0.01: + jj=12-j + fractions[jj]=1.-f + else: + jj=j + fractions[jj]=f + if jj==12: + C_11[jj]=0.01 + C_11_DEF[jj]=0.01 + C_11_CONT[jj]=0.01 + C_11_EF[jj]=0.01 + C_22[jj]=0.01 + C_22_DEF[jj]=0.01 + C_22_CONT[jj]=0.01 + C_22_EF[jj]=0.01 + C_12[jj]=0.01 + C_12_DEF[jj]=0.01 + C_12_CONT[jj]=0.01 + C_12_EF[jj]=0.01 + C_33[jj]=0.01 + C_33_DEF[jj]=0.01 + C_33_CONT[jj]=0.01 + C_33_EF[jj]=0.01 + if f!=0. and f!=1.: + C_11[jj]+=ChomHH[j-1,0,0,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_11_DEF[jj]+=ChomCA[j-1,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_11_CONT[jj]+=ChomSA[j-1,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_11_EF[jj]+=ChomEF[j-1,0,0]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22[jj]+=ChomHH[j-1,0,0,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22_DEF[jj]+=ChomCA[j-1,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22_CONT[jj]+=ChomSA[j-1,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_22_EF[jj]+=ChomEF[j-1,1,1]/1.11111*np.sqrt(c)/10/nb-1/nb + C_12[jj]+=ChomHH[j-1,0,0,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_12_DEF[jj]+=ChomCA[j-1,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_12_CONT[jj]+=ChomSA[j-1,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_12_EF[jj]+=ChomEF[j-1,0,1]/0.277777*np.sqrt(c)/10/nb-1/nb + C_33[jj]+=ChomHH[j-1,0,0,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + C_33_DEF[jj]+=ChomCA[j-1,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + C_33_CONT[jj]+=ChomSA[j-1,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + C_33_EF[jj]+=ChomEF[j-1,2,2]/2/0.416666*np.sqrt(c)/10/nb-1/nb + + plt.figure() + g1,=plt.plot(fractions,C_11,'+-') + g2,=plt.plot(fractions,C_11_DEF,'+-') + g3,=plt.plot(fractions,C_11_CONT,'+-') + g11,=plt.plot(fractions,C_11_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_11 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_11 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_22,'+-') + g2,=plt.plot(fractions,C_22_DEF,'+-') + g3,=plt.plot(fractions,C_22_CONT,'+-') + g11,=plt.plot(fractions,C_22_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_22 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_22 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_12,'+-') + g2,=plt.plot(fractions,C_12_DEF,'+-') + g3,=plt.plot(fractions,C_12_CONT,'+-') + g11,=plt.plot(fractions,C_12_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_12 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_12 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + + plt.figure() + g1,=plt.plot(fractions,C_33,'+-') + g2,=plt.plot(fractions,C_33_DEF,'+-') + g3,=plt.plot(fractions,C_33_CONT,'+-') + g11,=plt.plot(fractions,C_33_EF,'+-') + plt.legend([g1,g2,g3,g11],['HH','C_CA','C_SA','EF']) + plt.title('C_33 contraste '+str(c)) + plt.xlabel('f_souple') + plt.ylabel('Cij/Cij_rigide') + plt.axis([0,1,0,1]) + plt.savefig('courbes/C_33 moyen (nb=1) contraste '+str(c)+'.pdf') + #plt.show() + plt.close() + +courbe_moyenne(100,1) +courbe_moyenne(0.01,1) +#for i in range(10): +# courbes(100,i) +# courbes(0.01,i) + +Id=np.array([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]]) + +def err(c,i): + os.system('mkdir courbes') + fi=h5py.File('c='+str(c)+'/c='+str(c)+'ech'+str(i)+'.hdf5','r') + Chom=fi['ChomHH'] + err_def_EF=fi['err_def_EF'] + err_cont_EF=fi['err_cont_EF'] + err_def=fi['err_def'] + err_cont=fi['err_cont'] + fr=fi['fractions'] + err_def_EF_11=np.zeros((11,)) + err_cont_EF_11=np.zeros((11,)) + err_def_11=np.zeros((11,)) + err_cont_11=np.zeros((11,)) + fractions=np.zeros((11,)) + for j,f in enumerate([0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.67,0.8,0.9,0.95]): + if c==0.01: + jj=10-j + fractions[10-j]=1-fr[j] + else: + jj=j + fractions[j]=fr[j] + err_def_EF_11[jj]=err_def_EF[j]*np.sqrt(c)*10/1.1111 + err_cont_EF_11[jj]=err_cont_EF[j]/np.sqrt(c)*10/0.96 + err_def_11[jj]=err_def[j][0,0]*np.sqrt(c)*10/1.1111 + err_cont_11[jj]=err_cont[j][0,0]/np.sqrt(c)*10/0.96 + + plt.figure() + g1,=plt.plot(fractions,err_def_11,'+-') + g2,=plt.plot(fractions,err_cont_11,'+-') + plt.legend([g1,g2],['E_HH_souple_11','E_HH_rigide_11']) + plt.title('Energie_HH_11_contraste '+str(c)+'ech '+str(i)) + plt.xlabel('f_souple') + plt.savefig('courbes/Energie_11_contraste '+str(c)+'ech '+str(i)+'.pdf') + plt.close() + + plt.figure() + g1,=plt.plot(fractions,err_def_EF_11,'+-') + g2,=plt.plot(fractions,err_cont_EF_11,'+-') + plt.legend([g1,g2],['E_EF_souple_11','E_EF_rigide_11']) + plt.title('Energie_EF_11_contraste '+str(c)+'ech '+str(i)) + plt.xlabel('f_souple') + plt.savefig('courbes/Energie_EF_11_contraste '+str(c)+'ech '+str(i)+'.pdf') + plt.close() + +#err(100,0) +#err(0.01,0) + + + + + + + + + + + + + diff --git a/explohh/execute.py b/explohh/execute.py new file mode 100644 index 0000000..d55b549 --- /dev/null +++ b/explohh/execute.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +import numpy as np +import h5py +import matplotlib +matplotlib.use('Agg') +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_elliptiques import image_init,image_augmente,image_diminue +from Champs_admissibles import champs_admissibles +import os + + +N=32 +D=0.03 +nb=1 + +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(nb): + os.system("mkdir c="+str(c)+"/ech="+str(i)) + print("échantillon "+str(i)+'\n') + fname = "c="+str(c)+"/c="+str(c)+"ech"+str(i)+".hdf5" + fi = h5py.File(fname, "w") + im,centers_diam,f1=image_init(N,0.03,0.03,0,True,0.95) + l=len(centers_diam) + fi['centers']=np.zeros((11,l,5)) + fi['images']=np.zeros((11,N,N)) + fi['fractions']=np.zeros((11,)) + fi['ChomEF']=np.zeros((11,3,3)) + fi['CVoigt']=np.zeros((11,3,3)) + fi['C_CA']=np.zeros((11,3,3)) + fi['C_SA']=np.zeros((11,3,3)) + fi['ChomHH']=np.zeros((11,3,3,3,3)) + fi['EpsHH']=np.zeros((11,N,N,3,3)) + fi['A_CA']=np.zeros((11,N,N,3,3)) + fi['B_SA']=np.zeros((11,N,N,3,3)) + fi['Chom_DEF']=np.zeros((11,3,3,3,3)) + fi['Chom_CONT']=np.zeros((11,3,3,3,3)) + + + 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)) + print("fraction volumique "+str(f)+'\n') + if f!=0.95: + f1=image_diminue(im,centers_diam,f) + fi['images'][10-j]=im + l=len(centers_diam) + for ii in range(l): + fi['centers'][10-j][ii]=centers_diam[ii] + fi['fractions'][10-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] + + CA,CDEF,CCONT,tA=homog_hierar(tC,'HH','HHbis') + fi['ChomHH'][10-j,0]=[CA,CDEF,CCONT] + fi['EpsHH'][10-j]=tA + # CA,CDEF,CCONT=homog_hierar(tC,'HH','DEFbis') + # fi['ChomHH'][10-j,1]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'HH','CONTbis') + # fi['ChomHH'][10-j,2]=[CA,CDEF,CCONT] + + # CA,CDEF,CCONT=homog_hierar(tC,'DEF','HHbis') + # fi['Chom_DEF'][10-j,0]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'DEF','DEFbis') + # fi['Chom_DEF'][10-j,1]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'DEF','CONTbis') + # fi['Chom_DEF'][10-j,2]=[CA,CDEF,CCONT] + + # CA,CDEF,CCONT=homog_hierar(tC,'CONT','HHbis') + # fi['Chom_CONT'][10-j,0]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'CONT','DEFbis') + # fi['Chom_CONT'][10-j,1]=[CA,CDEF,CCONT] + # CA,CDEF,CCONT=homog_hierar(tC,'CONT','CONTbis') + # fi['Chom_CONT'][10-j,2]=[CA,CDEF,CCONT] + + ttau=np.zeros((N,N,3,3)) + for ii in range(N): + for jj in range(N): + ttau[ii,jj]=np.dot(tC[ii,jj]-CA,tA[ii,jj]) + + C_CA,C_SA,A_CA,B_SA=champs_admissibles(tC,tA,ttau,CA,CA,32) + fi['C_CA'][10-j]=C_CA + fi['C_SA'][10-j]=C_SA + fi['A_CA'][10-j]=A_CA + fi['B_SA'][10-j]=B_SA + ChomEF,CVoigt = reference(tC,128,c,i,f) + fi['ChomEF'][10-j]=ChomEF + fi['CVoigt'][10-j]=CVoigt + + + + + + + + + + + + + + + + + + diff --git a/explohh/inclusions.py b/explohh/inclusions.py new file mode 100644 index 0000000..d1aceef --- /dev/null +++ b/explohh/inclusions.py @@ -0,0 +1,132 @@ +import numpy as np +import matplotlib.pyplot as plt + + +from random import uniform + + +def compteur(im): + nn,mm=im.shape + cpt=0 + for i in range(nn): + for j in range(mm): + if im[i,j]==1: + cpt=cpt+1 + return cpt/nn/mm + + +def ajoute(im,centers_diam,D,nb): + nn,mm=im.shape + ll=len(centers_diam) + for i in range(nb): + centers_diam.append([uniform(0,1),uniform(0,1),D]) + + for i in range(nb): + r=D/2 + 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 + 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 + 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): + + im=np.zeros((N,N)) + + n=int(f/np.pi/(D/2)**2) + + centers_diam=[] + c=0 + while cf0: + 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 + + + + + + + diff --git a/explohh/inclusions_elliptiques.py b/explohh/inclusions_elliptiques.py new file mode 100644 index 0000000..470bc8d --- /dev/null +++ b/explohh/inclusions_elliptiques.py @@ -0,0 +1,144 @@ +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + + +from random import uniform + + +def compteur(im): + nn,mm=im.shape + cpt=0 + for i in range(nn): + for j in range(mm): + if im[i,j]==1: + cpt=cpt+1 + return cpt/nn/mm + + +def ajoute(im,centers,a,b,th,iso,nb): + nn,mm=im.shape + ll=len(centers) + for i in range(nb): + if iso==True: + centers.append([uniform(0,1),uniform(0,1),uniform(0,np.pi),a,b]) + else: + centers.append([uniform(0,1),uniform(0,1),th,a,b]) + + for i in range(nb): + c=max(centers[ll+i][3],centers[ll+i][4]) + x=centers[ll+i][0] + y=centers[ll+i][1] + the=centers[ll+i][2] + for j in range(int((x-c)*nn)-1,int((x+c)*nn)+1): + for k in range(int((y-c)*mm)-1,int((y+c)*mm)+1): + x1=(j+0.5)/nn + y1=(k+0.5)/mm + x2=(x1-x)*np.cos(the)+(y1-y)*np.sin(the) + y2=(y1-y)*np.cos(the)-(x1-x)*np.sin(the) + if x2**2/a**2+y2**2/b**2<=1: + if j>=nn: + l=nn + else: + l=0 + if k>=mm: + m=mm + else: + m=0 + im[j-l,k-m]=1 + +def enleve(im,centers,nb): + nn,mm=im.shape + im2=np.zeros((nn,mm)) + + for i in range(nb): + centers.pop() + + ll=len(centers) + for i in range(ll): + a=centers[i][3] + b=centers[i][4] + c=max(a,b) + x=centers[i][0] + y=centers[i][1] + the=centers[i][2] + for j in range(int((x-c)*nn)-1,int((x+c)*nn)+1): + for k in range(int((y-c)*mm)-1,int((y+c)*mm)+1): + x1=(j+0.5)/nn + y1=(k+0.5)/mm + x2=(x1-x)*np.cos(the)+(y1-y)*np.sin(the) + y2=(y1-y)*np.cos(the)-(x1-x)*np.sin(the) + if x2**2/a**2+y2**2/b**2<=1: + 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,a,b,th,iso,f): + + im=np.zeros((N,N)) + + n=int(f/np.pi/a/b) + centers=[] + c=0 + while cf0: + print('fraction volumique actuelle :'+str(f0)) + else: + c=f0 + im2=im + while c>f: + l=len(centers) + cc=0. + i=0 + while cc < c-f: + a=centers[l-i-1][3] + b=centers[l-i-1][4] + cc+=np.pi*a*b + i+=1 + im2=enleve(im2,centers,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 + + + + + diff --git a/explohh/influence_tensors-nu=0.0-patch_size=02-max_depth=10.npy b/explohh/influence_tensors-nu=0.0-patch_size=02-max_depth=10.npy new file mode 100644 index 0000000000000000000000000000000000000000..eacec1c37c20863cdd55228dcd346c8b3e90ce49 GIT binary patch literal 18128 zcmeI233yFsyM|Z#H>Gn*l~$X^n1moC2wC)r4ByzjTx zHq7q${4??0{FJfE>%sl{kLo)zxM@(ZIv^@IEGT$D(#X*xdk^cGG_qfRb3L|q;;8<{ z`lvy@hxa!)q<(lZ?^NJbTiH{+xg|O`^MUO+Q^?IN&O(! zsgt`sZpNACd>wH%AkODwb!I>aGuR9c1c?@uT?@`3~yj?+;Ikd^h(i5w&+s#!H=Y z^^Z%k>efh{C`>&sbx$;TZ`!?`<3_Kx))S7Mj;3y2P0ljz^G}>DxR3XB-Ff_*bElR% z$o!m;<14PGbl&KoPG+3XTU0AIHFlaXughYm>8M{PQ->ceoownmW9RGik-jnxnL6*I zP96fEnfF5YGyX<_+PuPQ2YDa-vDW?9MpGxVF7|`{sj_HLrF!%IH8S;G6&!KdPi}Wm zC$nFypMG#Y$kb7XYre>)Zu;ZhuVMGZ)=0coa!bD^4G#I(e&D!+I+=CR$F7)X`bVbD zb?+<9akzpYib3Q765oUbY0 zReE!0{euqbWcG{m!*gJsn~-%ETy?Q-W7(ejWfhe=nRONX&F2Q^`8NLi<@3e7Sy$re zy)(-`f7MkN>wazX=rSGaeea-7W?kIJI9Vd1Q@nqXD-J*Br+o$GpUe6qn!1be{JEuH z8w=+rUb_31^ObNY@NjI~xzYL1$?RAB#&W$<<`>#>_rI@nKa}kZRLuY0a*#7LsRt4C&Pg%jY;bNu zoF5RUvfS?^-v|A?jPImQ-iSCCxz8)kd2j{jd%WV17b8v{_Fch!bCEx}6m_{0L|DPqOhH}%o$KD$t#KDd@84C#=RDt^x_N#SC0xUE z-lfZR{u-|F|8OqE_%%gc`mWB=)efv++CMthNU0;nq1dbr~mI zZExP!f4bm4-p9`&zlzV1>oJJKT`_70#93E!e%jYr?+WwUZQ;B6HtCN}rp|TmtEDfT zt9$#r6>;b@b>2_i#qgPVKjYz#qL9yd_|w4UudJ+r#=7u6_>cRA_l4z)4}V;K<-kAG zBT$bZ)1QCUFU|$}`bz)E)G^1dev&!Q^pU>G&#kJG*#}Mxbu#nee3(~r zet3S!x$uYU-dFnLyI(P`^M&`4>%7Q3crIczo-gLjbHIGb?1NnA`Qomd=PL`(g<`Mg z2AMvy->i#FpVy;)`oXy;(;u#TU+ItUeg$d1`^D!@5T6@bJDhj=!~5tXnV%!qdA{g# z1J@j~AAJ7NK{9o&dtZ6JQcJ^9C(J0eon1F@#PuNs&OHC&em#`A(d5!%srT|v@at-y z?4VB0EceL|%70z9&Aw#rD|CizG~;C3D(qbo)}$cQM!s-G>cxBHI(2f;qh_2U=UBw4 zh&T~4j+6X8;&j0}b#ev7c@J^wBThK>lfOcoidd&ku8TP15vMNVl$HB!9~_VZew8sPmNr$4oEF5(cUw#Mg=?|y}=zRwrmqx3t*b-sorG;8+o)be8Q_FsN}we1;a zxxkZ;1fE>uAcuC7dVu0*{Iol$lQTL_o_z9bC2`1#Q|z3)!;Dj4yMI=dLkBgdjoiAc z8K>C!G}ft;8(MKLJM$1{Fyd5}`<>)4#Ho*U>g10RM?;)Bh!cqYgX7wgo?9S}!H zoZg5ND)-ySk%)5-)~ScG*YXE!2;O#`VdK%@)@>{eI8sjhTYn6gtR}R^;Ublj$?(%<_wK7~fL*XRv2J5i$=$m-$7goWs#x^TWD1 zKlJA+{Gm?f{M7b&zHZkqo-6h%z}2t6cE0%j43J)6h3zojUm$#5sw4sw1B**iT-CIGeFfot%$2-LdZy@;M^++sNmjuf;lb zax27HwqlF6`0YZy!_cQr)u@v#Qui$xaXB|rke5T}y7%?7#p9FL8SlZya_+@<^7(`J zpsJB^Ztz|ULBCa8L$q1)IWHRDziLL$=X3O01L~ZN9K!ctq)JYDdFqAB6TIJpk(ftC zsbk((J6R0JzM39!0#tHCe2z7D_O8eyK*lk+Hsbs?YUKISA18}q^Llr(K12GWlg~lt zy7$%ES3O;)&u!uJdx%3`1fRL?eWlNhkaviCUG`{=JjC6f9XMZxZq;dUDC+XvuaMjH z3+G&6-vYFm=ojm^u7*F;F+b$r;ZF+Isgo^#4S(3LVE4L==MVWd{gQRL^@~1pe(^c3 z`Cj2U&pLI?wN9Uz59h`5??i2l?2%7Qyv)PUt$hYpgU=UHzjZbIL0<)#{&C&=YV9-p zX@v9jt;hKa(a4_X>(+jyd7UrLih-N3Uq@3XcinOo%s#Q7jnbpUB3JEnC81*Wu-65W4M;p*snRLe z?)}<7Le^t%BJ(l$0P5mAzXPAw!cX#G_?&}v>g4M1xxm-_TyE*gr-gey4adrS4BhhC z;Qh#l^-qC6i{Jx!KK#kXI(4$OuUWdj+SmEAN8F}gc)lCD?0a49OPsIP$h*LOe(bF^ zvgi3~EbljT8K<$f8+8r!n&;)n|2Fe1sAHZ*Rrov`^-qF7><4)U^iRC{MSk1T({;{| z^c~NQ$73|IN52}4RIPdp-MSjge(`+CbAmX%!d3E4_wzL}O2#pCD~`ci;E(U~73k{M zrykG4PDxg+wr;5L)8Qara7w-GJ3+?J3l8e! z8Y${a!NIQz`$lu$v-(p%nsI827ea=e-5U0^AZKrodK0WuCwHG}#_6UfSou7w$NnJW z#LE3H&i_Kj>CSar$$^M70r>?qF?{aiteHFB5ksx(Bz=MO zK_~Zvz6|Ts$x=_(=fTHA@cBdRC%+GUCDy5v_rvFsYSN#3XSw;4RY4`MaQEkhW~irv zyc+zVb>Ad`_p)cddI{BMzZBeW=+fs3D*N>Y>VFOXAYVc513gpLC8(39S$evzA#dqB z|6W*C$(7vu71L0YdEotSF#9#lYo49(lk>a?`c|w{C+EXw=JSX9{EQD)$(7ynvEQ)% z-LhinR$c}(pCr^j4E|)n2QuepIo7F@rGKvPQJ(!eA$*^&Ud^RHhHm9!F#F~6eC37s zK3{qGo@(e;9D^_6T+D`#hdt+6kkg^BLS2|=K`wyLufv}S@Mi`3LhcV8^P*EHCtLc> z`BB`?MbNU|D)|Zbe(k7c`C{p2VFq7AU6*{FuWwptKF`;vs?rxjxBM~qp?`M1f*Sbl zS6(A|zoAQ?8@c+mtHbfAMU`iZe+;>}?eT91>5o`HjylMN9i=|b{4w`Il{)!o|CXs= z?D#+|h%@(%(*H~{!qFhEjmV!KZ}=oL#$ILM+`LMysH1>8M*Hr#Azq@>*Ul< zQtyLx>f~1tCmH!%N1XQ9PhNq1RIF1c*Fu~}kk3@aX(9Iu@>Jy073peVEkY9qn73C zeO4Zl{Z^{u`jCd_kFdS@#J=ZBbM?z}*606nWrBQfo`_V*X%)7gse5jfTp!*=wc;7N zbv5|m(4>^Z#S=sd_GP%osnlI155wnJ_g(KR^DuPsuOTu{dgV8YhUp{4{s(P|>wCVf zeg0U*H&bQUb@43y7lGs_^Z|IiK2A9YDbcT$2v_bug;LlR@ zg**;=8rG?k$HAZJ$a}h5U3sUUQ_1fAQeGTv)nn>bT?U^;zb3$+6t_6`!);Wu`}qp0 zA^$$u)Xg|1*Z*hdE2y3Ce%YUv{u;WKkHPHMq0W=lqgU668G+C5ShQk-e6F4-tCGib zk^1xIk23h)V(8>O-Aj@?+^31g*1nl~lPAnLJ;kc7Z4&)<^cLjG(5GOXI(eBDCqw_% z%4eqjJK`*n`(6A6^bfI4o!tFJGoK8d`NSj6O6(`sL_R54r%qmke15>^6?oM}UITrW zjN_`C+!^`&cy;@WZ7%Q9Lym|eO*hXKXZ(g#y(@dZAb$dVH`b|>tNKSR*-+-NUPsQo z_((o~@E+XjAUB3?^*Mc*ls>q4Tbk5=z&dquM&l7brk`J;)5j$Ed>#AA2jKHhSf@_D za<}xk8}hz<%FUngeu1V-Ew|(awX6C z(Wzs8^rP^H=e?X;zw!cJ*2wR<_p6Kio?z;hj|R6#zy7zKFWKk+|MmZh1Ff!E?^E%M zJQwmkm=L6rS3x)Hu{XuPdpC4)H~Bslh@;~^`MygC(#Sm?Ju&C|z{S?O5l7Ckp(Bnq zro5@==q1?4`%GPNi-UZuJpYY2%;%Q6kPqWn^%$Sys_WsB_i|Iu=ga%7B>vtncgmkT x)i3z`xBah$*@mNXVgCY~>!SJ(e1BDNFsCI6nwhr{?!ZXww6@O09 z0f%1H%ET!QKPc-~N*oo8Q&m-_A9}5D6^c{msvJK!^kj+COV$lQ9A8{0UftqAH$OR! zql8NwXR*Y&M%E3MIKf}zj-7YlDnInK!ZYiaPU>8y zxLH3e^MgYp&eD#T*LbPx)3#X8+($*~yU^ftT_sxJ(1yR+^<~r-mg}tcIx*`UQe0=v zDUQ~tGxA@+b>cYNb~w=BfK4E0Hn`qY;CBvjT1{Xl<)B@UUD`d-)%4Zif1t3T>CXQ>|?dh%J$SHI_zT+|(X zF#Alu`eNFL;n2tnb;kTsUcvO(&yBov&P#Jzr#0#vyFMwDw?20sIQgmS>>o+tcv>4C zYkuO+i(k0A<(|{}H7+<}{UfP+{ov51-?5$>=DwjmzxZR-CzSIm&5hx-PHWT${Xlhlp=*t0z{!1Mb=Z)}gS@#X%s3iT>*?|rV zucgdazGJB2&;upTdRezu;#61WtB(V{Rd^j`zHab?L${MSSa%q49;~jvI{2aE)^Xl> zZOzw6$56weM@gK0>)JO7Y;rxZ{M^WDU_pb^9Og9TIn*3_%2QwOj66|V`Fua=Ia6|p z3w_dg3DsqP%_)x7axUziCx58~)US#i%bV@`_z&g?uu&PFkbRy`(;3 zE8{?epPUZjDAhyigMOetwBDBxyI&IPdfHKT+rX^4t*1 z^$9h1p}{FaavQF%=FphunDt5O7<+!^b!30dDUR00KYD)9AM6}DcE7MH!$MwI2S2QU z2FLYNE&3j#Uy};uCDZXgSd`I8vXjyirYq24E-{&dpX(g)aeagSeyR5AG(I{ z>MbT+Sie)D0}eemwdebBAN~}$VwLfC!|7}TMh2NSX#aOevt zPHDKYtg9_?va9xs>mR(^51k}DR5dRyM-@2W(1=q@*2#SXh4HCz-GieY=*xxgJ+C5e zLS%vi4t=id=leG~uV%8J3C>r0RvdlEfvzHPx_wkS`X4C+9dPI>5~uf1cdwZ`)b+&j zJkR;gfd)7Ij@T`D{xpZa;k+)-t@@zU+so&|7V8MD|yBJ||kJGxC@GB#u(h$ou5_ zuu@U=K|jzRUa`l42IuP~{Qj=?h%Ir@FXWGU@cN9{^-*}gIAO_O*9Uo_KO86X41Gtx zGNn!@`qfkFELQ3hw4lK`erxWR!uwBacYw?@>WI1`FKBSS-d6U3yxLkbWDcb+C(n{AtUSjOTxf8LqYlY)S$ZK?{kKZXqx91O4IkOB`r$x*qyk>b!)l*67#WQXjEgAC+oB zgLC{;*+;76m+BJN^@*8RD77*#aNbWt_SYObAbsB|brvgiwp?g%@+;5U)-`3G zaX+c2#8Cy;0S&|D4&Rh7A`=mKE z`i1j@^Ta&<4ChOJ@1x=+P6Lbb*XIlKiut8DzDq4=aJpWW!ueoqL)*Bl!2gSkh4h%3eWeXqa<=j+KFM{{VLi%xQW#7f@-E;P9K@jg^QY;tAa z(GT>8))&{o=dsG?CC8D!zt9{S{mMxHVqc5Q6g!Qrm*smXzDCwOKXhZ^4<7vdw&#vL z<5Rgy~o6{LtmXC;$1*mD>xR@qn z8x3i?J#%P**O;x$;SoPHIQgk7W&V!%p_>@L9!2&Q5=Sh@$!ldnrwJ#&aq@GV=FnG4 zoc&`)?>!h_WQ*mzp1j6|1}DEmz1d%LXyQqWyez4+Sl&-&7Z)0w{MO0)sOHeQQs;Ep zw^*ron7c3 zg_B>XJ^O18OT*AJm-VY;({r)DiQ-e6JrGoUTXk_osXLd})pO=`QmumhYS0!G-Q7oZ=+N z=S$6@bEM8PCpPBAaNqe|E9VE?To2umN&bdI?=|(+^+CUohv7>11Nx-AR5RWO`lU5- z-7l8mQG#XP z%hp@%>R0Bce(34KpMJ``7W=`Wtw}3#e>mF5UUA68N%2b9+J7~#sc-qAUl$%;tzR=E z4mfl@iqqU%CF|ajINQt2Yl|QHP2t_>^_2q-y+`6KmUVAPoJ;BTk^}vk@VQ0iwa5X7 zZYpt>$-3ndr>l95eaC^8{JaeFTDi#&4*jjfd1!Fw%j=gr-UPO0p7Wkv9dx0=_4Tq> zTQ#-7q1$X2c;}BhOTAlp?rpJr-m~jl(BRrnsRwFvzcq)J=e)gczx~0xYYM!zZ25cz z(_QHIg_B=~{QgmM==FwoD5{LJAu3Mr^&AUYEbTM2;b!tT+{7`uRN_3kZtLo$tDmr^ zv8DA9b-e{WRXE-MvfdQOaOlZ~AA8XYNu9-VUK{$j(BO3co5Wvp=uT4SwYZC$XmIjtlF$B{Llf8aiP^7` zi(~gIIf?x>r#M=pUmwYQ$T{P_th~~NmUG7YJdnozno}IDalSCGm|x=ayI9cRbp7s5 z?5{a=7pXJmM=aep;nhJr7ptG}z>VLQ*9K7>!=Y)ubbavL*cg2-3a_)E(NDU6gRUH;OC z4sv^3a<$Ncr<*?VFw&~fW*Om z>dAiS_ptB&%Maa9_Cx#Vb=U!iE|NI8a=&cZ59PgZzXN@p@Wl8IiSco3{ov38CC-M0 z+tsH(ukys|a}lcc%-DPj8eI3sOKCGy_lY_5jHWXt)T;WimoLv9u{1Nj`TIrLG9GkoqN#WUB=vBgqenXN}z(BQ@|)Mga<8xB2${9p7wl{#-?%Xxis zvkNVj;;$S%l=IQt#4&oF)VbO$uFv^wIZi?+7g{XEQHfV^9L-G}qpOo&k&S+!Kg19F zV)rZW8?LkFv`%aEYrNFwKDK<{UHLBb{ldxb6L}xh9J)a2gMNvnyn3M$>(4dN;Knal^H%aV9D2Ug`9Z1A1h!P4k#CH! zpv6-B!qSyzNQ3YpRb%&v=5~@ z#nHOPSP#X69@=ex!+$C{c_`6aOlY&n*EINlXAb?? zP2pVwaz6Hs%JW?;&9e%==|F><{)BSn=M>GMXYskP56)Tf?ykc#yhChh{v#LIE;Kmx zd3u26&~WJSyO{7YHs2vt>;Q!=m+{^aBy2o3mV*9_f6w-Za6gh^@7yr`KbDYcJ;NO!Hu7qoJ)Bb z4n0HaGe`1zIy$cz$rd#FY5YPZ^6v#ThaM++#mPL2r8-}nHr9d$*M7kn!|tK|7!Hkj zE|C39iQdoQ4i>anuCv#L;wa6{ew5ajpZBFcYoqF;>fdHTi=};r9v#eaG&f<4_M|@O z2l`{Us-5UUgVXh^=D&=QC^ghX({9HK?nnOPzbv`2Vd@yRBBU47Y z(9!ceuJahm$8ZzJ=sRSdF&|>Np0{UO&@xA~&$6BJbDZWTj?s@xeQ++ONB66}Ps)(ZwrwT{JgorQN2Ei8I1$$JRFQZ?)d6 z=|YR;_rZ7fIpFj@sP#mO^Q8BK#5o{w%zb=^n!3>C!p%Oe{;)CC#c=3J5~qi(%Z`c@ zdcBql-BY-^|IZi7=P}Kp$4i{|W!<`{IKkO4tVy zJ+Y>q*7UTIE_5T|`g&Q?;yX=V#GwzJKl18P7k}(MEAOLXX`VxW>g+<#uY8}?b-S|- z)yHrhhx7tI7k2A^wHaK#Zl0IPmiikhYH2~Y5pMj{b^QMg8V=oW!-)5`JoS`)iL6VC ziW3@flM8)mRGh@txA1;6H*t*KBXM5+=Q}e8yfxW=g)QZE=t!;$Etd8fTGor=7;fSi zU1!ACYnyI=#Ty}Y6l>N8vpNm2pux@c!D@Y_3FXtAEY9z3>#etshkoh$Pkk@|R1^~qS#*@Xtz zern4(gQ+ftLuW~S&@br=%}-wAzAkihzd}22pg4w`I7aW1e&PJ!JQ@DqJYRJG|NrCv zOa*qyc@S&nU)^z40~Z?HT)!Z-A$0JBw|_+SFuwok3;RYrOuYY~K4;7; zI)1FYPU+Vfx^-JZ+(6%FK|=Cvz@R1 E1xhBregFUf literal 0 HcmV?d00001 diff --git a/explohh/output.log b/explohh/output.log new file mode 100644 index 0000000..424025a --- /dev/null +++ b/explohh/output.log @@ -0,0 +1,1331 @@ +-------------------------------------------------------------------------- +[[48331,1],0]: A high-performance Open MPI point-to-point messaging module +was unable to find any relevant network interfaces: + +Module: OpenFabrics (openib) + Host: doubs + +Another transport will be used instead, although this may result in +lower performance. + +NOTE: You can disable this warning by setting the MCA parameter +btl_base_warn_component_unused to 0. +-------------------------------------------------------------------------- +contraste 100 + +échantillon 0 + +fraction volumique 0.95 + +itération 1 + +itération 2 + +itération 3 + +[[ 1.20506100e-02 3.03213076e-03 3.83812638e-06] + [ 3.03212249e-03 1.20599446e-02 4.29364467e-06] + [ 3.51664043e-06 4.16078647e-06 8.94335601e-03]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 1.35064818e-02 3.14347500e-03 -1.49845495e-04] + [ 3.14347500e-03 1.32187262e-02 -1.65261321e-05] + [ -1.49845495e-04 -1.65261321e-05 9.55219642e-03]] [[ 1.23190923e-02 3.05285562e-03 -7.30808103e-06] + [ 3.05285562e-03 1.22829176e-02 -4.48828346e-07] + [ -7.30808103e-06 -4.48828346e-07 9.08402652e-03]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 1.22394872e-02 3.03000033e-03 -2.65947608e-05] + [ 3.03000033e-03 1.21735126e-02 -1.08674746e-05] + [ -2.65947608e-05 -1.08674746e-05 9.02686508e-03]] +fraction volumique 0.9 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +[[ 1.33136694e-02 3.33772436e-03 1.37369178e-05] + [ 3.33748650e-03 1.31180086e-02 1.41407448e-05] + [ 1.43514295e-05 1.47674136e-05 9.66057262e-03]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.01961171 0.00524111 -0.00059529] + [ 0.00524111 0.01871487 -0.00021216] + [-0.00059529 -0.00021216 0.01241936]] [[ 1.37942800e-02 3.39733903e-03 -7.77396503e-05] + [ 3.39733903e-03 1.36667891e-02 -2.31619748e-05] + [ -7.77396503e-05 -2.31619748e-05 9.96940355e-03]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.01375099 0.00342112 -0.0001981 ] + [ 0.00342112 0.01357528 -0.00012 ] + [-0.0001981 -0.00012 0.0100388 ]] +fraction volumique 0.8 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +[[ 1.63607579e-02 4.14272969e-03 -1.68364665e-05] + [ 4.14254747e-03 1.62970266e-02 -1.12700674e-05] + [ -1.04604172e-05 -1.54370967e-05 1.14910949e-02]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.03793709 0.01340611 -0.00087197] + [ 0.01340611 0.04089012 -0.00059683] + [-0.00087197 -0.00059683 0.02358639]] [[ 1.76169608e-02 4.29588482e-03 -1.39948087e-04] + [ 4.29588482e-03 1.74484197e-02 9.78312491e-06] + [ -1.39948087e-04 9.78312491e-06 1.22651832e-02]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 1.77044440e-02 4.50613180e-03 -1.14974387e-04] + [ 4.50613180e-03 1.76591468e-02 7.87348666e-05] + [ -1.14974387e-04 7.87348666e-05 1.28366571e-02]] +fraction volumique 0.67 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +[[ 2.30288786e-02 5.93940160e-03 -9.82848421e-05] + [ 5.94019320e-03 2.32856721e-02 -8.35877856e-05] + [ -7.12323858e-05 -7.77374509e-05 1.53540969e-02]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.08778503 0.03291238 -0.00200441] + [ 0.03291238 0.0973406 -0.00062201] + [-0.00200441 -0.00062201 0.05109433]] [[ 2.50200542e-02 6.12461308e-03 -1.34478868e-04] + [ 6.12461308e-03 2.57013375e-02 2.47219073e-05] + [ -1.34478868e-04 2.47219073e-05 1.66894649e-02]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.02714007 0.00742558 -0.00019558] + [ 0.00742558 0.02866668 -0.00027959] + [-0.00019558 -0.00027959 0.02039034]] +fraction volumique 0.6 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +[[ 0.02914128 0.00770913 -0.00019001] + [ 0.00771316 0.03061247 -0.00014187] + [-0.00016039 -0.00014018 0.01899342]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 1.38953601e-01 5.73909614e-02 -4.53943016e-03] + [ 5.73909614e-02 1.55881092e-01 4.57579599e-05] + [ -4.53943016e-03 4.57579599e-05 7.61766122e-02]] [[ 3.02379095e-02 7.71102708e-03 -3.78007833e-04] + [ 7.71102708e-03 3.38941964e-02 -9.79799763e-05] + [ -3.78007833e-04 -9.79799763e-05 2.01015825e-02]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.03601354 0.01081144 -0.00118214] + [ 0.01081144 0.04048323 -0.00211822] + [-0.00118214 -0.00211822 0.03019989]] +fraction volumique 0.5 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +itération 8 + +[[ 0.04559025 0.01254356 -0.00032087] + [ 0.01257494 0.0523807 -0.00016416] + [-0.00021485 -0.00014642 0.02861612]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 2.20288457e-01 9.09922926e-02 -4.70864930e-03] + [ 9.09922926e-02 2.73508577e-01 1.34229628e-04] + [ -4.70864930e-03 1.34229628e-04 1.27274775e-01]] [[ 0.04064449 0.00962148 -0.00080135] + [ 0.00962148 0.04668952 -0.00091778] + [-0.00080135 -0.00091778 0.02592773]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.05744661 0.01730969 -0.00447371] + [ 0.01730969 0.06744436 -0.00682377] + [-0.00447371 -0.00682377 0.05074681]] +fraction volumique 0.4 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +itération 8 + +itération 9 + +[[ 0.09323755 0.02732066 0.00036177] + [ 0.02743382 0.11773368 0.00098446] + [ 0.00072443 0.00091199 0.05815082]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.31812223 0.13968403 -0.01236708] + [ 0.13968403 0.44293175 -0.00829174] + [-0.01236708 -0.00829174 0.20600573]] [[ 0.05489906 0.01206026 -0.00275523] + [ 0.01206026 0.06404606 -0.00186044] + [-0.00275523 -0.00186044 0.03041697]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.10481109 0.03480211 -0.01407209] + [ 0.03480211 0.14563794 -0.01661331] + [-0.01407209 -0.01661331 0.09754844]] +fraction volumique 0.3 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +[[ 0.25091143 0.06888563 0.00047704] + [ 0.06922942 0.2902174 0.00220874] + [ 0.00031787 0.00170907 0.14999365]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.44348825 0.14531303 -0.00705179] + [ 0.14531303 0.50616611 -0.00717347] + [-0.00705179 -0.00717347 0.27311641]] [[ 0.07825596 0.01588329 -0.00011448] + [ 0.01588329 0.0902405 0.00095932] + [-0.00011448 0.00095932 0.03877313]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.23886549 0.05510129 -0.01297212] + [ 0.05510129 0.27461161 -0.01494107] + [-0.01297212 -0.01494107 0.1719581 ]] +fraction volumique 0.2 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +[[ 0.50895212 0.13281598 -0.00202797] + [ 0.13316091 0.54830711 -0.00126715] + [-0.00200236 -0.00145824 0.31934938]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.56338248 0.14516971 -0.00541436] + [ 0.14516971 0.62460199 -0.00452366] + [-0.00541436 -0.00452366 0.34504133]] [[ 0.15309262 0.02418984 -0.00161604] + [ 0.02418984 0.17696687 -0.0018827 ] + [-0.00161604 -0.0018827 0.07182198]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.47221428 0.10978332 -0.0057239 ] + [ 0.10978332 0.55903856 -0.00280136] + [-0.0057239 -0.00280136 0.3040896 ]] +fraction volumique 0.1 + +itération 1 + +itération 2 + +itération 3 + +[[ 8.15107752e-01 2.06776506e-01 5.84764677e-04] + [ 2.06728494e-01 8.13921686e-01 -7.95501907e-04] + [ 7.27659571e-04 -5.93043836e-04 5.57293893e-01]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.80387696 0.20230302 0.00469742] + [ 0.20230302 0.80215092 0.00308563] + [ 0.00469742 0.00308563 0.53719726]] [[ 0.42059926 0.06959841 0.00077292] + [ 0.06959841 0.39244323 -0.00166296] + [ 0.00077292 -0.00166296 0.20360261]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.78637979 0.18934682 0.00383187] + [ 0.18934682 0.78498331 0.00271657] + [ 0.00383187 0.00271657 0.53785694]] +fraction volumique 0.05 + +itération 1 + +itération 2 + +[[ 0.9572323 0.24237474 0.0020749 ] + [ 0.24239508 0.97447351 0.00215128] + [ 0.00233705 0.00241706 0.6932377 ]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 0.94349294 0.23845454 -0.00683209] + [ 0.23845454 0.95268614 -0.00641909] + [-0.00683209 -0.00641909 0.67778623]] [[ 0.66222882 0.0931501 0.03344322] + [ 0.0931501 0.64077709 0.03823771] + [ 0.03344322 0.03823771 0.33048857]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 0.93273587 0.23113421 -0.01119715] + [ 0.23113421 0.93883768 -0.01177015] + [-0.01119715 -0.01177015 0.67576004]] +contraste 0.01 + +échantillon 0 + +fraction volumique 0.95 + +itération 1 + +itération 2 + +[[ 9.64195180e+01 2.43644685e+01 6.61818526e-02] + [ 2.43626404e+01 9.60844295e+01 4.38560125e-02] + [ 6.37947300e-02 4.10741491e-02 6.90661613e+01]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 94.74352214 23.63572275 0.41641093] + [ 23.63572275 94.23454631 0.86987355] + [ 0.41641093 0.86987355 66.9508794 ]] [[ 75.46696365 16.17647343 -0.30970396] + [ 16.17647343 74.91938068 -2.35717516] + [ -0.30970396 -2.35717516 42.78186134]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 94.82305357 23.34949709 0.71685442] + [ 23.34949709 94.29804716 1.23217834] + [ 0.71685442 1.23217834 67.69247701]] +fraction volumique 0.9 + +itération 1 + +itération 2 + +itération 3 + +[[ 8.10651438e+01 2.05507986e+01 6.64357251e-02] + [ 2.05496938e+01 8.05093989e+01 3.57070319e-02] + [ 1.00783375e-01 4.24596847e-02 5.50106752e+01]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 78.7473458 19.98366634 1.07514869] + [ 19.98366634 77.9063362 1.89273322] + [ 1.07514869 1.89273322 53.27918455]] [[ 41.93737993 6.92654194 0.6075019 ] + [ 6.92654194 41.88182203 -1.36158401] + [ 0.6075019 -1.36158401 18.77238793]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 77.48953869 19.14340926 1.64602683] + [ 19.14340926 75.84639641 2.32871395] + [ 1.64602683 2.32871395 53.59983199]] +fraction volumique 0.8 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +[[ 51.17726365 13.31262957 -0.16674085] + [ 13.30559131 51.88425186 -0.25461485] + [ -0.17957313 -0.27542769 32.17639091]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 51.48282313 13.85684036 1.11033002] + [ 13.85684036 55.41939346 1.01387292] + [ 1.11033002 1.01387292 34.76142169]] [[ 16.91586692 3.21657036 -0.16853747] + [ 3.21657036 19.0887927 -0.13630693] + [ -0.16853747 -0.13630693 8.57675477]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 44.94823386 11.48757872 1.3353861 ] + [ 11.48757872 51.49803225 1.07060757] + [ 1.3353861 1.07060757 32.83459011]] +fraction volumique 0.67 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +itération 8 + +[[ 18.13813791 4.79706159 0.17955058] + [ 4.78406627 17.6269658 0.17962449] + [ 0.19680614 0.1720494 10.56651341]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 38.36172107 12.41797 -1.89120798] + [ 12.41797 38.35402929 -1.34374734] + [ -1.89120798 -1.34374734 23.81104818]] [[ 7.36881196 1.41866689 0.07894991] + [ 1.41866689 7.95015448 -0.0633751 ] + [ 0.07894991 -0.0633751 4.01885852]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 20.80201856 5.91416931 -2.48607981] + [ 5.91416931 20.20210652 -2.69595604] + [ -2.48607981 -2.69595604 14.93819981]] +fraction volumique 0.6 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +itération 8 + +itération 9 + +[[ 9.94533947 2.58387348 0.08956304] + [ 2.5805794 9.61453303 0.11879433] + [ 0.10402415 0.12566766 5.57532818]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 33.49677798 12.0379258 -0.16448221] + [ 12.0379258 34.27881881 -0.63386792] + [ -0.16448221 -0.63386792 18.94249126]] [[ 5.67547011 1.26571623 0.11948652] + [ 1.26571623 6.18933683 -0.1743594 ] + [ 0.11948652 -0.1743594 3.29823516]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 12.48569206 3.56017345 -0.42806788] + [ 3.56017345 13.9100877 -1.3478446 ] + [ -0.42806788 -1.3478446 8.67132599]] +fraction volumique 0.5 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +itération 8 + +[[ 4.87131757 1.24376115 0.06080735] + [ 1.24054817 4.72284505 0.05659635] + [ 0.07209931 0.06876481 2.84624979]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 23.33356463 8.65832668 -0.20232172] + [ 8.65832668 23.43339591 -0.16107321] + [ -0.20232172 -0.16107321 12.32349202]] [[ 4.36266096 0.92973496 -0.06212624] + [ 0.92973496 4.21314111 -0.11885985] + [-0.06212624 -0.11885985 2.54880754]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 6.65573865 1.58386997 -0.49239928] + [ 1.58386997 6.46541963 -0.57064421] + [-0.49239928 -0.57064421 4.34394107]] +fraction volumique 0.4 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +itération 6 + +itération 7 + +[[ 2.94315701 0.76597308 0.01591055] + [ 0.76575 3.00146059 0.01802857] + [ 0.02148001 0.02356865 1.8849416 ]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 13.54450592 5.14908136 -0.45661682] + [ 5.14908136 14.52403835 -0.56925408] + [ -0.45661682 -0.56925408 7.53521918]] [[ 3.11814416 0.69780213 -0.09448432] + [ 0.69780213 3.11328819 -0.13406135] + [-0.09448432 -0.13406135 1.95453335]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 3.88680228 0.99962568 -0.3121757 ] + [ 0.99962568 3.95691427 -0.4756405 ] + [-0.3121757 -0.4756405 2.63003775]] +fraction volumique 0.3 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +itération 5 + +[[ 2.08501102 0.53786325 0.00680349] + [ 0.53772345 2.09146673 0.00441852] + [ 0.00847065 0.00617137 1.41956661]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 6.71346705 2.31407929 -0.2371023 ] + [ 2.31407929 7.02480578 -0.17787526] + [-0.2371023 -0.17787526 4.0924178 ]] [[ 2.2111521 0.52494885 -0.04596848] + [ 0.52494885 2.23471853 -0.05389568] + [-0.04596848 -0.05389568 1.50326509]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 2.29471977 0.56594921 -0.09542988] + [ 0.56594921 2.35970617 -0.0964012 ] + [-0.09542988 -0.0964012 1.69871715]] +fraction volumique 0.2 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +[[ 1.6133298 0.40978823 0.00220145] + [ 0.40987102 1.62492241 0.00174642] + [ 0.00273227 0.00241142 1.14004998]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 3.7018822 1.2839443 -0.01294962] + [ 1.2839443 3.99655444 -0.09780641] + [-0.01294962 -0.09780641 2.3628141 ]] [[ 1.67428481 0.41347534 -0.0105579 ] + [ 0.41347534 1.70561132 -0.02033437] + [-0.0105579 -0.02033437 1.18582854]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 1.68582413 0.41821164 -0.01514296] + [ 0.41821164 1.74027567 -0.0263997 ] + [-0.01514296 -0.0263997 1.24341215]] +fraction volumique 0.1 + +itération 1 + +itération 2 + +itération 3 + +itération 4 + +[[ 1.31826828e+00 3.32297425e-01 1.64923506e-04] + [ 3.32299978e-01 1.32177239e+00 3.63845605e-05] + [ 4.48565323e-04 4.76870828e-04 9.63029525e-01]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 1.82513995 0.59020382 -0.00802537] + [ 0.59020382 2.04003882 -0.04418683] + [-0.00802537 -0.04418683 1.28889564]] [[ 1.34174534 0.33581633 -0.00469455] + [ 0.33581633 1.36265134 -0.00904841] + [-0.00469455 -0.00904841 0.98187662]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 1.33495227 0.33459582 -0.00542238] + [ 0.33459582 1.36471052 -0.01106408] + [-0.00542238 -0.01106408 0.99397336]] +fraction volumique 0.05 + +itération 1 + +itération 2 + +itération 3 + +[[ 1.20152506e+00 3.01870316e-01 -4.99462212e-04] + [ 3.01869230e-01 1.20248510e+00 -5.15056346e-04] + [ -4.79331909e-04 -4.68926219e-04 8.91304681e-01]] +Solving linear variational problem. +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +[[ 1.39939393 0.40098653 -0.0137074 ] + [ 0.40098653 1.46677039 -0.02616372] + [-0.0137074 -0.02616372 0.99440416]] [[ 1.21323813 0.30286247 -0.00135932] + [ 0.30286247 1.22029599 -0.00211928] + [-0.00135932 -0.00211928 0.89842494]] +Solving linear variational problem. +Solving linear variational problem. +Solving linear variational problem. +[[ 1.21143102 0.30168995 -0.00364948] + [ 0.30168995 1.22253109 -0.00514989] + [-0.00364948 -0.00514989 0.90090052]] -- GitLab