Reference.py 6.61 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
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



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



#cf 'Contenu des fichiers' à propos de 'reference'
def reference(tabC,res1,c,ec,fr):
    #res=256
    k,ll,mm,nn=tabC.shape
    mesh = UnitSquareMesh(res1, res1, "crossed")

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

Jeremy BLEYER's avatar
Jeremy BLEYER committed
31 32
        def value_shape(self):
            return (3, 3)
33 34 35



Jeremy BLEYER's avatar
Jeremy BLEYER committed
36 37 38 39 40 41 42
#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]
43

Jeremy BLEYER's avatar
Jeremy BLEYER committed
44 45 46 47 48 49
    # plt.imshow(tabC[:,:,0,0])
    # plt.show()
    # plt.close()
    # plt.imshow(Cv[:,:,0,0])
    # plt.show()
    # plt.close()
50

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

Jeremy BLEYER's avatar
Jeremy BLEYER committed
54 55 56 57 58 59 60 61 62 63 64 65
    # 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
66

Jeremy BLEYER's avatar
Jeremy BLEYER committed
67 68 69 70 71 72 73
        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)
74

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

86 87 88


#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
89 90 91 92 93 94
    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)])
95

Jeremy BLEYER's avatar
Jeremy BLEYER committed
96 97
        def sigma(v, Eps):
            return dot(C, eps(v)+Eps)
98 99


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

Jeremy BLEYER's avatar
Jeremy BLEYER committed
106 107 108
        v_,lamb_ = TestFunctions(W)
        dv, dlamb = TrialFunctions(W)
        w = Function(W)
109

Jeremy BLEYER's avatar
Jeremy BLEYER committed
110 111 112 113
        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
114

Jeremy BLEYER's avatar
Jeremy BLEYER committed
115 116 117 118 119
        Chom = np.zeros((3, 3))
        for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
            ee = np.zeros((3,))
            ee[j] = 1
            Eps.assign(Constant(ee))
120
            solve(a == L, w, solver_parameters={"linear_solver": "mumps"})#, solver_parameters={"linear_solver": "bicgstab", "preconditioner":"ilu"})
Jeremy BLEYER's avatar
Jeremy BLEYER committed
121 122 123 124 125 126
            (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.)
127

Jeremy BLEYER's avatar
Jeremy BLEYER committed
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
            plt.figure()
            pp=plot(eps(v)[0])
            plt.colorbar(pp)
            plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 1'+str(j+1)+'.pdf')
            #plt.show()
            plt.close()
            plt.figure()
            pp=plot(eps(v)[1])
            plt.colorbar(pp)
            plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 2'+str(j+1)+'.pdf')
            #plt.show()
            plt.close()
            plt.figure()
            pp=plot(eps(v)[2])
            plt.colorbar(pp)
            plt.savefig('c='+str(c)+'/ech='+str(ec)+'/frac='+str(fr)+'/Aloc EF 3'+str(j+1)+'.pdf')
            #plt.show()
            plt.close()
146

Jeremy BLEYER's avatar
Jeremy BLEYER committed
147 148 149 150
            Sigma = np.zeros((3,))
            for kk in range(3):
                Sigma[kk] = assemble(sigma(v, Eps)[kk]*dx)
            Chom[j, :] = Sigma
151

Jeremy BLEYER's avatar
Jeremy BLEYER committed
152 153
        ffile.close()
        return Chom
154

Jeremy BLEYER's avatar
Jeremy BLEYER committed
155 156 157 158 159 160 161
    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]
162

Jeremy BLEYER's avatar
Jeremy BLEYER committed
163
    CVoigt=CVoigt/k/k
164

Jeremy BLEYER's avatar
Jeremy BLEYER committed
165 166 167 168 169 170
    print(Chom)
    return Chom,CVoigt




171

Jeremy BLEYER's avatar
Jeremy BLEYER committed
172