From 25047a446e9e711ac9dca0c3aa835f982f1771f5 Mon Sep 17 00:00:00 2001 From: Jeremy Bleyer Date: Tue, 19 May 2020 11:11:42 +0200 Subject: [PATCH] Increase to 256 resolution for FE with Mumps --- Reference.py | 52 ++++++++++++++++++++++++++-------------------------- execute.py | 2 +- 2 files changed, 27 insertions(+), 27 deletions(-) diff --git a/Reference.py b/Reference.py index 25ac533..4b7b2d1 100644 --- a/Reference.py +++ b/Reference.py @@ -27,12 +27,12 @@ def reference(tabC,res1,c,ec,fr): 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): @@ -40,17 +40,17 @@ def reference(tabC,res1,c,ec,fr): # for i in range(int(i0*res/k),int((i0+1)*res/k)): # for j in range (int(j0*res/k),int((j0+1)*res/k)): # Cv[i,j]=tabC[i0,j0] - + # plt.imshow(tabC[:,:,0,0]) # plt.show() # plt.close() # plt.imshow(Cv[:,:,0,0]) # plt.show() # plt.close() - + #Cvar est la transformation du tableau Cv en une fonction définie sur 'mesh' Cvar= LocalC(tabC,mesh) - + # class used to define the periodic boundary map class PeriodicBoundary(SubDomain): def __init__(self, vertices, tolerance=DOLFIN_EPS): @@ -63,7 +63,7 @@ def reference(tabC,res1,c,ec,fr): # 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 @@ -71,7 +71,7 @@ def reference(tabC,res1,c,ec,fr): 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]) @@ -82,49 +82,49 @@ def reference(tabC,res1,c,ec,fr): 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' + + +#la fonction reference_problem prend en argument C, donnant le tenseur d'élasticité en un point de 'mesh', et renvoie le Chom du domaine. k étant le nombre de phases différentes sur le domaine, qui correspond à la taille du tableau Cn pris en argument par la fonction 'reference' def reference_problem(C, mesh,k): def eps(v): return as_vector([v[0].dx(0), v[1].dx(1), (v[0].dx(1)+v[1].dx(0))/sqrt(2)]) - + def sigma(v, Eps): return dot(C, eps(v)+Eps) - - + + Ve = VectorElement("CG", mesh.ufl_cell(), 2) Re = VectorElement("R", mesh.ufl_cell(), 0) vertices = np.array([[0, 0], [1, 0], [1, 1], [0, 1]]) W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10)) V = VectorFunctionSpace(mesh,'CG',degree=2,dim=3) - + v_,lamb_ = TestFunctions(W) dv, dlamb = TrialFunctions(W) w = Function(W) - + Eps = Constant((0,)*3) F = inner(sigma(dv, Eps), eps(v_))*dx a, L = lhs(F), rhs(F) a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx - + Chom = np.zeros((3, 3)) for (j, case) in enumerate(["Exx", "Eyy", "Exy"]): ee = np.zeros((3,)) ee[j] = 1 Eps.assign(Constant(ee)) - solve(a == L, w, [])#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"}) + solve(a == L, w, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"}) (v, lamb) = split(w) ffile = XDMFFile('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Strain EF '+case+'.xdmf') ffile.parameters["functions_share_mesh"] = True vf = Function(V, name='Strain-'+case) vf.assign(project(eps(v),V)) ffile.write(vf,0.) - + plt.figure() pp=plot(eps(v)[0]) plt.colorbar(pp) @@ -143,15 +143,15 @@ def reference(tabC,res1,c,ec,fr): plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf') #plt.show() plt.close() - + Sigma = np.zeros((3,)) for kk in range(3): Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx) Chom[j, :] = Sigma - + ffile.close() return Chom - + Chom = reference_problem(Cvar, mesh,k) CVoigt = np.array([[0., 0., 0.], [0., 0., 0.], @@ -159,14 +159,14 @@ def reference(tabC,res1,c,ec,fr): 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/execute.py b/execute.py index b2bd231..ed475b9 100644 --- a/execute.py +++ b/execute.py @@ -65,7 +65,7 @@ for c in [100,0.01]: if im[ii,jj]==1: tC[ii,jj]=C[1] - ChomEF,CVoigt = reference(tC,128,c,i,f) + ChomEF,CVoigt = reference(tC,256,c,i,f) fi['ChomEF'][j]=ChomEF fi['CVoigt'][j]=CVoigt -- GitLab