Commit 25047a44 authored by Jeremy BLEYER's avatar Jeremy BLEYER

Increase to 256 resolution for FE with Mumps

parent 2897e76a
...@@ -27,12 +27,12 @@ def reference(tabC,res1,c,ec,fr): ...@@ -27,12 +27,12 @@ def reference(tabC,res1,c,ec,fr):
p = Cell(mesh, ufc_cell.index).midpoint() p = Cell(mesh, ufc_cell.index).midpoint()
i, j = int(p[0]*(self.Nx)), int(p[1]*(self.Ny)) i, j = int(p[0]*(self.Nx)), int(p[1]*(self.Ny))
value[:] = self.Cij[j, i,:,:].flatten() value[:] = self.Cij[j, i,:,:].flatten()
def value_shape(self): def value_shape(self):
return (3, 3) return (3, 3)
#Ici on crée, à partir du tableau des tenseurs d'élasticité correspondant au domaine dont on veut calculer le Chom, un tableau Cv représentant le même Cn mais adapté à une résolution plus fine 'res' #Ici on crée, à partir du tableau des tenseurs d'élasticité correspondant au domaine dont on veut calculer le Chom, un tableau Cv représentant le même Cn mais adapté à une résolution plus fine 'res'
# Cv = np.array([[np.zeros((3,3)) for i in range(res)] for j in range(res)]) # Cv = np.array([[np.zeros((3,3)) for i in range(res)] for j in range(res)])
# for i0 in range(k): # for i0 in range(k):
...@@ -40,17 +40,17 @@ def reference(tabC,res1,c,ec,fr): ...@@ -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 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)): # for j in range (int(j0*res/k),int((j0+1)*res/k)):
# Cv[i,j]=tabC[i0,j0] # Cv[i,j]=tabC[i0,j0]
# plt.imshow(tabC[:,:,0,0]) # plt.imshow(tabC[:,:,0,0])
# plt.show() # plt.show()
# plt.close() # plt.close()
# plt.imshow(Cv[:,:,0,0]) # plt.imshow(Cv[:,:,0,0])
# plt.show() # plt.show()
# plt.close() # plt.close()
#Cvar est la transformation du tableau Cv en une fonction définie sur 'mesh' #Cvar est la transformation du tableau Cv en une fonction définie sur 'mesh'
Cvar= LocalC(tabC,mesh) Cvar= LocalC(tabC,mesh)
# class used to define the periodic boundary map # class used to define the periodic boundary map
class PeriodicBoundary(SubDomain): class PeriodicBoundary(SubDomain):
def __init__(self, vertices, tolerance=DOLFIN_EPS): def __init__(self, vertices, tolerance=DOLFIN_EPS):
...@@ -63,7 +63,7 @@ def reference(tabC,res1,c,ec,fr): ...@@ -63,7 +63,7 @@ def reference(tabC,res1,c,ec,fr):
# check if UC vertices form indeed a parallelogram # 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[3, :] - self.a1) <= self.tol
assert np.linalg.norm(self.vv[2, :]-self.vv[1, :] - self.a2) <= self.tol assert np.linalg.norm(self.vv[2, :]-self.vv[1, :] - self.a2) <= self.tol
def inside(self, x, on_boundary): def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the # return True if on left or bottom boundary AND NOT on one of the
# bottom-right or top-left vertices # bottom-right or top-left vertices
...@@ -71,7 +71,7 @@ def reference(tabC,res1,c,ec,fr): ...@@ -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 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 (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) (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): 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 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[0] = x[0] - (self.a1[0]+self.a2[0])
...@@ -82,49 +82,49 @@ def reference(tabC,res1,c,ec,fr): ...@@ -82,49 +82,49 @@ def reference(tabC,res1,c,ec,fr):
else: # should be on top boundary else: # should be on top boundary
y[0] = x[0] - self.a2[0] y[0] = x[0] - self.a2[0]
y[1] = x[1] - self.a2[1] 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 reference_problem(C, mesh,k):
def eps(v): def eps(v):
return as_vector([v[0].dx(0), return as_vector([v[0].dx(0),
v[1].dx(1), v[1].dx(1),
(v[0].dx(1)+v[1].dx(0))/sqrt(2)]) (v[0].dx(1)+v[1].dx(0))/sqrt(2)])
def sigma(v, Eps): def sigma(v, Eps):
return dot(C, eps(v)+Eps) return dot(C, eps(v)+Eps)
Ve = VectorElement("CG", mesh.ufl_cell(), 2) Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0) Re = VectorElement("R", mesh.ufl_cell(), 0)
vertices = np.array([[0, 0], [1, 0], [1, 1], [0, 1]]) vertices = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10)) W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))
V = VectorFunctionSpace(mesh,'CG',degree=2,dim=3) V = VectorFunctionSpace(mesh,'CG',degree=2,dim=3)
v_,lamb_ = TestFunctions(W) v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W) dv, dlamb = TrialFunctions(W)
w = Function(W) w = Function(W)
Eps = Constant((0,)*3) Eps = Constant((0,)*3)
F = inner(sigma(dv, Eps), eps(v_))*dx F = inner(sigma(dv, Eps), eps(v_))*dx
a, L = lhs(F), rhs(F) a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
Chom = np.zeros((3, 3)) Chom = np.zeros((3, 3))
for (j, case) in enumerate(["Exx", "Eyy", "Exy"]): for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
ee = np.zeros((3,)) ee = np.zeros((3,))
ee[j] = 1 ee[j] = 1
Eps.assign(Constant(ee)) Eps.assign(Constant(ee))
solve(a == L, w, [])#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"}) solve(a == L, w, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"})
(v, lamb) = split(w) (v, lamb) = split(w)
ffile = XDMFFile('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Strain EF '+case+'.xdmf') ffile = XDMFFile('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Strain EF '+case+'.xdmf')
ffile.parameters["functions_share_mesh"] = True ffile.parameters["functions_share_mesh"] = True
vf = Function(V, name='Strain-'+case) vf = Function(V, name='Strain-'+case)
vf.assign(project(eps(v),V)) vf.assign(project(eps(v),V))
ffile.write(vf,0.) ffile.write(vf,0.)
plt.figure() plt.figure()
pp=plot(eps(v)[0]) pp=plot(eps(v)[0])
plt.colorbar(pp) plt.colorbar(pp)
...@@ -143,15 +143,15 @@ def reference(tabC,res1,c,ec,fr): ...@@ -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.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf')
#plt.show() #plt.show()
plt.close() plt.close()
Sigma = np.zeros((3,)) Sigma = np.zeros((3,))
for kk in range(3): for kk in range(3):
Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx) Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx)
Chom[j, :] = Sigma Chom[j, :] = Sigma
ffile.close() ffile.close()
return Chom return Chom
Chom = reference_problem(Cvar, mesh,k) Chom = reference_problem(Cvar, mesh,k)
CVoigt = np.array([[0., 0., 0.], CVoigt = np.array([[0., 0., 0.],
[0., 0., 0.], [0., 0., 0.],
...@@ -159,14 +159,14 @@ def reference(tabC,res1,c,ec,fr): ...@@ -159,14 +159,14 @@ def reference(tabC,res1,c,ec,fr):
for i in range(k): for i in range(k):
for j in range(k): for j in range(k):
CVoigt=CVoigt+tabC[i,j] CVoigt=CVoigt+tabC[i,j]
CVoigt=CVoigt/k/k CVoigt=CVoigt/k/k
print(Chom) print(Chom)
return Chom,CVoigt return Chom,CVoigt
...@@ -65,7 +65,7 @@ for c in [100,0.01]: ...@@ -65,7 +65,7 @@ for c in [100,0.01]:
if im[ii,jj]==1: if im[ii,jj]==1:
tC[ii,jj]=C[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['ChomEF'][j]=ChomEF
fi['CVoigt'][j]=CVoigt fi['CVoigt'][j]=CVoigt
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment