Reference.py 6.87 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2 3 4 5 6 7 8 9 10 11 12



from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter



Jeremy BLEYER's avatar
Jeremy BLEYER committed
13 14 15 16 17 18 19 20 21 22 23
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

Jeremy BLEYER's avatar
Jeremy BLEYER committed
24 25 26 27
#cf 'Contenu des fichiers' à propos de 'reference'
def reference(tabC,res1,c,ec,fr):
    #res=256
    k,ll,mm,nn=tabC.shape
Jeremy BLEYER's avatar
Jeremy BLEYER committed
28
    mesh = UnitSquareMesh(res1, res1, "left/right")
Jeremy BLEYER's avatar
Jeremy BLEYER committed
29 30 31 32 33 34 35 36 37 38 39 40

    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()
41

Jeremy BLEYER's avatar
Jeremy BLEYER committed
42 43
        def value_shape(self):
            return (3, 3)
44 45


Jeremy BLEYER's avatar
Jeremy BLEYER committed
46 47 48 49 50 51 52
#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]
53

Jeremy BLEYER's avatar
Jeremy BLEYER committed
54 55 56 57 58 59
    # plt.imshow(tabC[:,:,0,0])
    # plt.show()
    # plt.close()
    # plt.imshow(Cv[:,:,0,0])
    # plt.show()
    # plt.close()
60

Jeremy BLEYER's avatar
Jeremy BLEYER committed
61 62
#Cvar est la transformation du tableau Cv en une fonction définie sur 'mesh'
    Cvar= LocalC(tabC,mesh)
63

Jeremy BLEYER's avatar
Jeremy BLEYER committed
64 65 66 67 68 69 70 71 72 73 74 75
    # 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
76

Jeremy BLEYER's avatar
Jeremy BLEYER committed
77 78 79 80 81 82 83
        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)
84

Jeremy BLEYER's avatar
Jeremy BLEYER committed
85 86 87 88 89 90 91 92 93 94 95
        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]

96 97 98


#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'
Jeremy BLEYER's avatar
Jeremy BLEYER committed
99 100 101 102 103 104
    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)])
105

Jeremy BLEYER's avatar
Jeremy BLEYER committed
106 107
        def sigma(v, Eps):
            return dot(C, eps(v)+Eps)
108 109


Jeremy BLEYER's avatar
Jeremy BLEYER committed
110 111 112 113
        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))
Jeremy BLEYER's avatar
Jeremy BLEYER committed
114
        V = VectorFunctionSpace(mesh,'DG',degree=1,dim=3)
115

Jeremy BLEYER's avatar
Jeremy BLEYER committed
116 117 118
        v_,lamb_ = TestFunctions(W)
        dv, dlamb = TrialFunctions(W)
        w = Function(W)
119

Jeremy BLEYER's avatar
Jeremy BLEYER committed
120 121 122 123
        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
124

Jeremy BLEYER's avatar
Jeremy BLEYER committed
125 126
        Chom = np.zeros((3, 3))
        for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
Jeremy BLEYER's avatar
Jeremy BLEYER committed
127 128
            ffile = XDMFFile('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Strain_EF_{}.xdmf'.format(case))

Jeremy BLEYER's avatar
Jeremy BLEYER committed
129 130 131
            ee = np.zeros((3,))
            ee[j] = 1
            Eps.assign(Constant(ee))
132
            solve(a == L, w, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"})
Jeremy BLEYER's avatar
Jeremy BLEYER committed
133 134
            (v, lamb) = split(w)
            vf = Function(V, name='Strain-'+case)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155
            vf.assign(local_project(eps(v), V))
            ffile.write_checkpoint(vf, case, 0.)

            # plt.figure()
            # pp=plot(vf[0])
            # plt.colorbar(pp)
            # plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 1'+str(j+1)+'.pdf')
            # #plt.show()
            # plt.close()
            # plt.figure()
            # pp=plot(vf[1])
            # plt.colorbar(pp)
            # plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 2'+str(j+1)+'.pdf')
            # #plt.show()
            # plt.close()
            # plt.figure()
            # pp=plot(vf[2])
            # plt.colorbar(pp)
            # plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf')
            # #plt.show()
            # plt.close()
156

Jeremy BLEYER's avatar
Jeremy BLEYER committed
157 158 159 160
            Sigma = np.zeros((3,))
            for kk in range(3):
                Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx)
            Chom[j, :] = Sigma
161

Jeremy BLEYER's avatar
Jeremy BLEYER committed
162
            ffile.close()
Jeremy BLEYER's avatar
Jeremy BLEYER committed
163
        return Chom
164

Jeremy BLEYER's avatar
Jeremy BLEYER committed
165 166 167 168 169 170 171
    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]
172

Jeremy BLEYER's avatar
Jeremy BLEYER committed
173
    CVoigt=CVoigt/k/k
174

Jeremy BLEYER's avatar
Jeremy BLEYER committed
175 176 177 178 179 180
    print(Chom)
    return Chom,CVoigt




181

Jeremy BLEYER's avatar
Jeremy BLEYER committed
182