// TP Eléments Finis 1D
// H. Haddar 

// Résolution par EF P2 de l'equation
// - d/dx (k d/dx) u + alpha u = f
// u(1)=u(0) = 0
//
// et etude de la convergence 
//
clear; // suppression de toutes les variables
close();

// Boucle sur le nombre de points de discrétisations
// N = 2^(3+p) , p=1, ..., 6
N = 9; Np = 9;
h = zeros(1,Np);
errL2 = zeros(1,Np);
errH1 = zeros(1,Np);
for p=1:Np

N=2^(p+1) + 1; // Nb de points de discrétisation

x=linspace(0,1,N); // points x_i répartis de manière uniforme
eps = 0.0/N; // eps doit être < Dx/2!
veps = eps*(2*rand(x)-1); // nombres aléatoires entre -eps et eps
x(2:N-1) = x(2:N-1)+veps(2:N-1); // On perturbe le maillage uniforme 

hx = x(2:N) - x(1:N-1) ; // les h_i
h(p) = max(hx); 
k = ones(hx); // le coefficient de conductivité
alpha = 0; // le coefficient alpha

xex=linspace(0,1,10000); 
//uexex = sin(4*%pi * xex); 
uexex = xex.*((0.5-xex).^2).*(xex<0.5);

//============ EF P2 ==================================

Nddl = (N-2) + (N-1) ;
xddl(2:2:2*(N-2)) = x(2:N-1);
xddl(1:2:2*(N-2)+1) = (x(1:N-1)+x(2:N))/2.;
//================ Solution exacte
//Exemple 1
//uex = sin(4*%pi * xddl); 
//f = (16*%pi^2+alpha)*uex;

// Exemple 2
uex = xddl.*((0.5-xddl).^2).*(xddl<0.5);
f = (-6*(xddl)+2*ones(xddl)).*(xddl<0.5)+alpha*uex;
//
 

//================ Marices de masse et de rigidité
Kjj = zeros(1,Nddl);
Kjj(2:2:2*(N-2)) = (7.0/3.0)*(k(1:N-2)./hx(1:N-2) + k(2:N-1)./hx(2:N-1)); 
Kjj(1:2:2*(N-2)+1)= (16.0/3.0)*(k(1:N-1)./hx(1:N-1))

Kjjp1=zeros(1,Nddl-1);
Kjjp1(2:2:2*(N-2)) = - (8.0/3.0)*k(2:N-1)./hx(2:N-1);
Kjjp1(1:2:2*(N-2)) = - (8.0/3.0)*k(1:N-2)./hx(1:N-2); 

Kjjp2=zeros(1,Nddl-2);
Kjjp2(2:2:2*(N-3)) = (1.0/3.0)*k(2:N-2)./hx(2:N-2);
Kjjp2(1:2:2*(N-2)) = 0; 

Kh = diag(sparse(Kjj)) + diag(sparse(Kjjp1),1) + diag(sparse(Kjjp1),-1) ...
                       + diag(sparse(Kjjp2),2) + diag(sparse(Kjjp2),-2); //Kh

Mjj = zeros(1,Nddl);
Mjj(2:2:2*(N-2)) = (2.0/15.0)*(hx(1:N-2) + hx(2:N-1)); 
Mjj(1:2:2*(N-2)+1)= (8.0/15.0)*(hx(1:N-1))

Mjjp1=zeros(1,Nddl-1);
Mjjp1(2:2:2*(N-2)) = (1.0/15.)*hx(2:N-1);
Mjjp1(1:2:2*(N-2)) = (1.0/15.)*hx(1:N-2); 

Mjjp2=zeros(1,Nddl-2);
Mjjp2(2:2:2*(N-3)) = -(1.0/30.)*hx(2:N-2);
Mjjp2(1:2:2*(N-2)) = 0; 

Mh = diag(sparse(Mjj)) + diag(sparse(Mjjp1),1) + diag(sparse(Mjjp1),-1) ...
                       + diag(sparse(Mjjp2),2) + diag(sparse(Mjjp2),-2); //Mh

Ah=Kh+alpha*Mh; // matrice du systeme

Fh=(Mh*f')'; // Le second membre 

// ================ Résolution du système linéaire 

spcho = chfact(Ah); // factorisation de Cholesky
Uh=chsolve(spcho,Fh); // resolution du systeme

clf()
plot(xex,uexex,"b");
plot2d(xddl, [Uh' ], leg="Uh @ U exact");
xtitle(["Pas de discrétisation = ", string(h(p))])
xpause(1000000);
scf
// ======== On cacule l'erreur 
errL2(p) = ((Uh-uex)*Mh*(Uh-uex)');
errH1(p) = ((Uh-uex)*Kh*(Uh-uex)');
end
//
err1 = sqrt(errH1);
err0=sqrt(errL2);

clf();
dN=log2(h.^1);
dN1=log2(h.^(1.5));
dN3=log2(h.^3);
dN4=log2(h.^(2.5));
dN2=log2(h.^2);
dN5=log2(h.^(3.5));
dN6=log2(h.^4);
logh = log2(h); logerr1= log2(err1);logerr0=log2(err0) ;

//plot2d(logh,[logerr1'    dN4' dN3' dN5']);
plot2d(logh,[logerr1'    dN' dN1' dN2']);

xtitle("courbe de convergence Kh")
scf
//plot2d(logh,[logerr0'  dN3' dN5' dN6'  ]);
plot2d(logh,[logerr0'    dN' dN1' dN2']);

xtitle("courbe de convergence Mh")

