# -*- coding: utf-8 -*-
"""
Created on Thu Jan 29 15:21:13 2015
Listing 13.4 : Equation de la chaleur, résolution par 
éléments finis

@author: ribot/grivet
"""

from pylab import *

#   discrétisation en espace

xx=array([0, 0.025, 0.05, 0.125, 0.175, 0.2, 0.25, 0.3, 0.4, 0.425, 0.5, 0.575,\
 0.6, 0.675, 0.725, 0.75, 0.8, 0.85, 0.95, 0.975, 1.0])
xx = xx.transpose()
npt = len(xx); nx=npt-2
hx=xx[1:npt]-xx[0:npt -1]
xxint=xx[1:nx+1]

#  matrice de discrétisation

K1D = zeros((nx,nx)); K1D[0,0] = 1.0/hx[0]
M1D = zeros((nx,nx)); M1D[0,0] = hx[0]/3.0
for i in range(1,nx):
    K1D[i-1,i-1] += 1.0/hx[i]; K1D[i,i] = 1.0/hx[i]
    K1D[i-1,i] = -1.0/hx[i]; K1D[i,i-1] = K1D[i-1,i]
    M1D[i-1,i-1] += hx[i]/3.0; M1D[i,i] = hx[i]/3.0
    M1D[i-1,i] = hx[i]/6.0; M1D[i,i-1] = M1D[i-1,i]
K1D[nx-1,nx-1] += 1.0/hx[nx]
M1D[nx-1,nx-1] += hx[nx]/3.0

#discrétisation en temps

Tfin = 1.0; dt = 0.9*min(hx)**2/10.0
ntps = int(Tfin/dt)
u0 = u = hstack((xxint[xxint <= 0.5],1.0 - xxint[xxint > 0.5]))
u = copy(u0)
for t in range(ntps):
    v = solve(M1D,dot(K1D,u)); u = u - dt*v;
uu = hstack([0,u,0]); uu0 = hstack([0,u0,0]) 
figure(1)
plot(xx,uu0)
figure(2)
plot(xx,uu*1e4)