Homog_hierar.py 5.59 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 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160


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,CVoigt):
    
    # K0=np.array([[0.5,-0.5, 0.        ],
    #    [-0.5, 0.5, 0.        ],
    #    [0.        , 0.        , 1.]])
    # J0=np.array([[0.5, 0.5, 0.        ],
    #    [0.5, 0.5, 0.        ],
    #    [0.        , 0.        , 0.]])
    
    pp,qq,rr,ss=tabC.shape
    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,j,l,k]
                        gamma12[n,3*i+l,3*j+k]=t12[9-n,i,j,l,k]
    
    
    
    
    
    #la fonction auxiliary_problem calcule les A locaux : l'argument principal est C_list, la liste des tenseurs d'élasticité, et la fonction renvoie Aloc, un tableau à trois dimensions. Aloc[i] correspond à la moyenne du tenseur de localisation sur la case i, calculé par homogénéisation périodique
    def auxiliary_problem(C_list, n,gamma,K):
        
        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.]])
        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)])
        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_))        
        
        return Chom,E
            
        
#la fonction tC prend en argument un tableau de type tCA et renvoie le tableau de tenseurs d'élasticité qui sera celui de l'étape d'après
    def tC(tC_,tA,n,gamma,K):
        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=auxiliary_problem(C_list, n,gamma,K)
                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])
                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])
                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])
                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])
                    
        return t1
        
  
    

#la fonction Chom prend en argument le tableau initial tC0 des tenseurs d'élasticité associés à la grille image de taille n, et appelle les fonctions décrites plus haut pour calculer les nouveaux tableaux de tenseurs d'élasticité associés à chaque étape. La fonction enregistre ces tableaux à chaque étape dans un fichier. Le nombre d'étapes 'e' doit être précisé pour des méthodes 10 ou 11.   

    def Chom(tC0,n,K_):
        # 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)])
        for i in range(n):
            t2=tC(t,tA,n-i,gamma,K)
            t=t2
        return(K,t[0,0],tA)
       
    
    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=Chom(tabC,p,C_hom)
        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(C0,C_hom1)
    
        
    return C_hom_list,tA_list