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

// Résolution par EF P1 de l'equation
// - d/dx (k d/dx) u + alpha u = f
// u(1)=u(0) = 0

// Etude de la convergence 
//
clear; // suppression de toutes les variables
g1=-sqrt(3/5)
g2=0
g3=sqrt(3/5)
w1=5/9
w2=8/9
w3=5/9
// Boucle sur le nombre de points de discrétisations
// N = 2^(3+p) , p=1, ..., 6
N = 9; Np = 7;
h = zeros(1,Np);
errL2 = zeros(1,Np);
errH1 = zeros(1,Np);
errL2ap = zeros(1,Np);
errH1ap = zeros(1,Np);
for p=1:Np

N=2^(p+3) + 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,4000); 
//uexex = sin(4*%pi * xex); 
uexex = xex.*((0.5-xex).^2).*(xex<0.5);
//uexex= (1-xex).*xex;

//============ EF P1 ==================================

Nddl = N-2 ;
xddl = x(2:N-1);

//================ 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;
//
// == Exemple 3
//uex = (1-xddl).*xddl;
//f = 2*ones(xddl)+alpha*uex;
//

//================ Marices de masse et de rigidité

Kjj = k(1:N-2)./hx(1:N-2) + k(2:N-1)./hx(2:N-1); // la digonale de Kh
Kjjp = - k(2:N-2)./hx(2:N-2); // les sur et sous diagonales de Kh

Kh = diag(sparse(Kjj)) + diag(sparse(Kjjp),1) + diag(sparse(Kjjp),-1); //Kh

Mjj  = (hx(1:N-2) + hx(2:N-1))/3.;// la digonale de Mh
Mjjp = hx(2:N-2)/6.;// les sur et sous diagonales de Mh

Mh = diag(sparse(Mjj)) + diag(sparse(Mjjp),1) + diag(sparse(Mjjp),-1); //Kh

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

Fh=(hx(1:N-2) + hx(2:N-1)).*f /2.; // le second membre avec formule des trapèzes
//Fh = (Mh*f')'; // façon plus précise pour discrétiser 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 
errL2ap(p) = ((Uh-uex)*Mh*(Uh-uex)');
errH1ap(p) = ((Uh-uex)*Kh*(Uh-uex)');
errH1(p)=0
errL2(p)=0
// par méthode exacte : avec une formule de quadrature exacte jusqu au degré 5



xg1(1:N-1)=0.5*(x(1:N-1)+x(2:N))+0.5*g1*hx(1:N-1)
up1(1:N-1)=([Uh(1:N-2) ,0]-[0 ,Uh(1:N-2)])./hx(1:N-1)
u1(1:N-1)=[Uh(1:N-2),0].*(xg1(1:N-1)-x(1:N-1))+[0,Uh(1:N-2)].*(x(2:N)-xg1(1:N-1))
//exemple 1
//fp1=(4*%pi*cos(4*%pi*xg1)-up1).^2
//f1=(sin(4*%pi*xg1)-u1./hx).^2
//exemple 2
fp1=((0.5-xg1).*(0.5-3*xg1).*(xg1<0.5)-up1).^2
f1=(xg1.*((0.5-xg1).^2).*(xg1<0.5)-u1./hx).^2
//exemple 3
//fp1=(1-2*xg1-up1).^2
//f1=((1-xg1).*xg1-u1./hx).^2
xg2(1:N-1)=0.5*(x(1:N-1)+x(2:N))+0.5*g2*hx(1:N-1)
up2(1:N-1)=([Uh(1:N-2),0]-[0 ,Uh(1:N-2)])./hx(1:N-1)
u2(1:N-1)=[Uh(1:N-2),0].*(xg2(1:N-1)-x(1:N-1))+[0 ,Uh(1:N-2)].*(x(2:N)-xg2(1:N-1))
//exemple 1
//fp2=(4*%pi*cos(4*%pi*xg2)-up2).^2
//f2=(sin(4*%pi*xg2)-u2./hx).^2
//exemple 2
fp2=((0.5-xg2).*(0.5-3*xg2).*(xg2<0.5)-up2).^2
f2=(xg2.*((0.5-xg2).^2).*(xg2<0.5)-u2./hx).^2
//exemple 3
//fp2=(1-2*xg2-up2).^2
//f2=((1-xg2).*xg2-u2./hx).^2
xg3(1:N-1)=0.5*(x(1:N-1)+x(2:N))+0.5*g3*hx(1:N-1)
up3(1:N-1)=([Uh(1:N-2),0]-[0,Uh(1:N-2)])./hx(1:N-1)
u3(1:N-1)=[Uh(1:N-2),0].*(xg3(1:N-1)-x(1:N-1))+[0 ,Uh(1:N-2)].*(x(2:N)-xg3(1:N-1))
//exemple 1
//fp3=(4*%pi*cos(4*%pi*xg3)-up3).^2
//f3=(sin(4*%pi*xg3)-u3./hx).^2
//exemple 2
fp3=((0.5-xg3).*(0.5-3*xg3).*(xg3<0.5)-up3).^2
f3=(xg3.*((0.5-xg3).^2).*(xg3<0.5)-u3./hx).^2
//exemple 3
//fp3=(1-2*xg3-up3).^2
//f3=((1-xg3).*xg3-u3./hx).^2
errH1(p)=0.5*(w1*fp1+w2*fp2+w3*fp3)*hx'
errL2(p)=0.5*(w1*f1+w2*f2+w3*f3)*hx'
end
//
err1 = sqrt(errH1);
err0= sqrt(errL2);
err1ap = sqrt(errH1ap);
err0ap= sqrt(errL2ap);
clf();
dN=log2(h);
dN12=log2(h.^0.5);
dN14=log2(h.^0.25);
dN2=log2(h.^2);
logh = log2(h); logerr1= log2(err1);logerr0=log2(err0) ;
logerr1ap= log2(err1ap);logerr0ap=log2(err0ap) ;
plot2d(logh,[logerr1ap' logerr0ap' dN' dN2']);
xtitle("courbe de convergence normes approchées")
xclick()
scf
plot2d(logh,[logerr1ap'  dN' dN2' ]);
xtitle("courbe de convergence Kh h h^2 ")
xclick
scf
plot2d(logh,[logerr0ap'  dN' dN2' ]);
xtitle("courbe de convergence Mh h h^2 ")
xclick()
scf
plot2d(logh,[logerr1' logerr0' dN' dN2']);
xtitle("courbe de convergence")
xclick()
scf
plot2d(logh,[logerr1'  dN' dN2' ]);
xtitle("courbe de convergence err1 h h^2 ")
xclick
scf
plot2d(logh,[logerr0'  dN' dN2' ]);
xtitle("courbe de convergence err0 h h^2 ")
xclick()
//xclick()
//pente= (logerr(2:Np)-logerr(1:Np-1))./(logh(2:Np)-logh(1:Np-1));
//scf()
//plot(pente);
//xtitle("pentes de la courbe de convergence")
