import numpy as np
import matplotlib.pyplot as plt

h=100
lamb=220 # plaque en aluminium
Text=20

Lx=10e-2 # taille du domaine de travail en mètre
Ly=10e-2

Nx=40  # nombre de points de la grille
Ny=40

dx=Lx/Nx # intervalle entre deux points
dy=Ly/Ny

precision = 1e-3# règle la précision du calcul : si la methode converge à ce
                 # nombre près, le programme arrête les itérations

T=np.zeros((Ny,Nx))       # tableaux à 2 dimensions vide au départ
T2=np.zeros((Ny,Nx))      # T température
B=np.zeros((Ny,Nx))       # jx et jy vecteur courant
jx=np.zeros((Ny,Nx))      # B : si B=1, la taempérature est imposée, B=0 vide
jy=np.zeros((Ny,Nx))

# conditions aux limites : plaque
T[:,:]=20
B[:,0]=1
T[:,0]=80

B[:,Nx-1]=2

B[0,:]=1
T[0,:]=60
B[Ny-1,:]=1
T[Ny-1,:]=50

#------------calcul du potentiel méthode de Gauss Seidel -------------------
w=2/(1+np.pi/max(Nx,Ny))        # paramètre optimal pour la convergence
emax=1000                       # on met une valeur arbitraire > précision
k=0                             # pour que le while ne stoppe pas tout de suite
                                # les itérations, emax sera l'écart max de
                                # potentiel entre deux itérations

while emax>precision:           # tant que la présicion n'est pas atteinte
                                # on met à zéro au départ pour calculer le max
    emax=0                      # d'écart sur l'ensemble de la grille

    for i in range(Ny):
        for j in range(Nx):

            if B[i,j]==0: # si on est dans le vide : on résout l'equ de Laplace
                          # méthode de Gauss -Seidel : converge + vite
                T2=T[i,j]
                T[i,j]=(1.-w)*T[i,j]+w*(T[i+1,j]+T[i-1,j]+T[i,j+1]+T[i,j-1])/4.
                emax=max(emax,abs(T2-T[i,j])) # on compare pour trouver
                                                   # l'écart maximum

        T[i,Nx-1]= (h*Text+lamb/dx*T[i,Nx-2])/(h+lamb/dx) # convection à droite

    k=k+1                 # on compte le nombre d'itérations
#---------------------------------------------------------------------------

print(str(k)+" itérations")

# -------------calcul des courants thermiques jx et jy--------------------------

for i in range(1,Ny-1):
    for j in range(1,Nx-1):

        jy[i,j]=-(T[i+1,j]-T[i-1,j])/(2*dy) # calcul du gradient selon Oy
        jx[i,j]=-(T[i,j+1]-T[i,j-1])/(2*dx) # selon Ox

#----------------------------------------------------------------------------

# tracé de la carte : T et vecteur courant j

x=np.linspace(0,Lx,Nx)
y=np.linspace(0,Ly,Ny)
X,Y=np.meshgrid(x,y)
plt.figure(figsize=(10,8))
cf=plt.contourf(X,Y,T,100,cmap='jet')
graph=plt.contour(X,Y,T,12,colors='black')
plt.clabel(graph,inline=1,fontsize=10,fmt='%3.2f') #affiche les valeurs de T
plt.colorbar(cf)
plt.xlabel("en m")
plt.ylabel("en m")
plt.title("Carte thermique")
plt.streamplot(X,Y,jx,jy,density=0.5,color='black')
plt.show()