# -*- coding: utf-8 -*-
"""
Created on Tue Jan 27 14:45:13 2015

Listing 13.2 : équation de Poisson, résolution par
éléments finis
-u(x)'' = x(1 - x)

@author: ribot/grivet
"""

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
from numpy.linalg import solve

xx=np.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])
npt = len(xx); nx=npt-2; 
hx=xx[1:npt]-xx[0:npt -1]
xxint=xx[1:nx+1]

#  Matrice de discrétisation

K1D = np.zeros((nx,nx)); F = np.zeros(nx)
K1D[0,0]=1/hx[0]
#uu = lambda x : x*(1-x)*(x-xx[0])/hx[0]
def fn0(x,a,b):
    return x*(1-x)*(x-a)/b
a = xx[0];b = hx[0]
F[0],err = quad(fn0,xx[0],xxint[0],args = (a,b))
def fnim1(x,a,b):
    return x*(1-x)*(a-x)/b
def fni(x,a,b):
    return x*(1-x)*(x-a)/b
for i in range(1,nx):
    K1D[i-1,i-1] = K1D[i-1,i-1] + 1/hx[i] 
    K1D[i,i] = 1/hx[i]
    K1D[i-1,i] = -1/hx[i] 
    K1D[i,i-1] = K1D[i-1,i]
    a = xxint[i]; b = hx[i];
    F[i-1],err = F[i-1] + quad(fnim1,xxint[i-1],xxint[i],args = (a,b))
    a = xxint[i-1]; b = hx[i]
    F[i],err = quad(fni,xxint[i -1],xxint[i], args = (a,b))
K1D[nx-1,nx-1] = K1D[nx-1,nx-1]+1/hx[nx]
def fnn(x,a,b):
    return x*(1-x)*(a-x)/b
a = xx[nx+1]; b = hx[nx]
F[nx-1] = F[nx-1] + quad(fnn,xxint[nx-1] ,xx[nx+1], args = (a,b))[0];

#Résolution du système

u = solve(K1D,F)
uu = np.hstack([0,u,0]);
plt.plot(xx,uu);