add explo2 explohh

explo2/Champs_bis.py

+
+
+
+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
+
+
+
+
+
+

explo2/Homog_hierar.py

+
+
+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
+
+
+
+    
+    
+
+
+
+
+

explo2/Reference.py

+
+
+
+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
+