clear();
//=========================================================================================
//Mise en oeuvre des éléments finis Pk de Lagrange en une dimension
//Résolution du problème aux limites
// -u'' + Lambda*u = f sur ]0,1[
//u(0) = u(1) = 0
//Etude de la convergence en fonction du pas de subdivision de ]0,1[, UNIFORME
//=========================================================================================

//=========================================================================================
//Saisie des données:
//      - f : Fonction "second memmbre"
//      - Lambda : Fonction "paramètre"
//      - uex : Solution exacte
//      - duex : Dérivée de la solution exacte
//      - k : Degré des éléments finis k, égal à 1 ou 2
//      - H : Liste des pas de subdivision de ]0,1[, VECTEUR COLONNE
//=========================================================================================

function y=f(x)
    y=6*abs(x-0.5);
endfunction

function y=Lambda(x)
    y=0;
endfunction

function y=uex(x)
    y=-abs(x-0.5)^3+1/8;
endfunction

function y=duex(x)
    y=-3*abs(x-0.5).*(x-0.5);
endfunction

k = 1;
H = 2.^(-[5:0.2:7])';

//=========================================================================================
//Fonction de création du maillage
//=========================================================================================

function [Nint,Next,Nodes,Cells]=createMesh(List,k)
    //Créée le maillage en fonction des points et du degré k
    Np = length(List);
    select k
    case 1 then
        Nint = Np;
        Next = 2;
        Nodes = [List;0;1];
        Cells = [[Np+1;(1:Np)'],[(1:Np)';Np+2]];
    case 2 then
        Nint = 2*Np+1;
        Next = 2;
        Nodes = [List;0.5*([0;List] + [List;1]);0;1];
        Cells = [[2*Np+2;(1:Np)'],[(1:Np)';2*Np+3],(Np+1:2*Np+1)'];
    end
endfunction

//=========================================================================================
//Saisie des fonctions de base et de leurs dérivées
//=========================================================================================

function S = RefBasisFct(x,k)
    //Fonctions de base pour l'élément de référence
    select k
    case 1 then
        S = [1-x;x];
    case 2 then
        S = [2*(x-0.5).*(x-1);2*x.*(x-0.5);4*x.*(1-x)];
    end
endfunction

function DS = DerRefBasisFct(x,k)
    //Dérivées des fonctions de base pour l'élément de référence
    select k
    case 1 then
        DS = [0.*x-1;0.*x+1];
    case 2 then
        DS = [4*x-3;4*x-1;4-8*x];
    end
endfunction

//=========================================================================================
//Création des vacteurs d'erreur
//=========================================================================================

nh = length(H);
ErrL2 = zeros(nh,1); //Vecteur d'erreur en norme L2
ErrH1 = zeros(nh,1); //Vecteur d'erreur en norme H1


//=========================================================================================
//Calcul des vecteurs ErrL2 et ErrH1. L'étude de la convergence commence maintenant
//=========================================================================================
disp("Calcul des vecteurs d erreur.");

for hh=1:nh //On parcourt tous les pas de subdivision de H
disp("Progression: "+string(round(100*hh/nh))+"%");

h = H(hh);
List = [(h:h:1-h)']; //Subdivision uniforme

//=========================================================================================
//Création du maillage
//=========================================================================================

[Nint,Next,Nodes,Cells]=createMesh(List,k);

//=========================================================================================
//Calcul de l'application affine F_K
//=========================================================================================

function y=FK(x,K)
    //Application affine envoyant l'élément de référence [0,1] vers la maille K, K étant ici le
    //numéro de la maille, et x est dans [0,1]
    y = (Nodes(Cells(K,2)) - Nodes(Cells(K,1)))*x + Nodes(Cells(K,1));
endfunction

//=========================================================================================
//Construction des matrices élémentaires et du second membre élémentaire
//Utilise la méthode de quadrature de Gauss-Legendre à quatre points pour le calcul
//d'intégrales
//=========================================================================================

Q = 4;
xq = [-sqrt(3/7+2/7*sqrt(6/5)),-sqrt(3/7-2/7*sqrt(6/5)),sqrt(3/7-2/7*sqrt(6/5)),sqrt(3/7+2/7*sqrt(6/5))]; //Points de la quadrature sur ]-1,1[
wq = [(18-sqrt(30))/36,(18+sqrt(30))/36,(18+sqrt(30))/36,(18-sqrt(30))/36]; //Poids associés à chaque point de quadrature
xxq = 0.5*xq + 0.5; //Pour adapter les points de quadrature à l'élément de référence [0,1]
wwq = 0.5*wq; //Pour adapter les poids de quadrature à l'élémente de référence [0,1]

function y = m(alpha,betaa,K)
    //Matrice élémentaire associée à la maille K, renvoie l'élément d'indice (alpha,betaa), où
    //alpha,betaa sont compris entre 1 et k+1=size(Cells)(2)
    mesK = Nodes(Cells(K,2)) - Nodes(Cells(K,1)); //Mesure de K, qui correspond à la dérivée de F_K, en valeur absolue
    y = (1/mesK)*sum(wwq.*DerRefBasisFct(xxq,k)(alpha,:).*DerRefBasisFct(xxq,k)(betaa,:)) + mesK*sum(wwq.*Lambda(FK(xxq)).*RefBasisFct(xxq,k)(alpha,:).*RefBasisFct(xxq,k)(betaa,:)); //Approximation d'intégrale
endfunction

function y = b(betaa,K)
    //Second membre élémentaire associée à la maille K, renvoie l'élément d'indice alpha, où
    //alpha est compris entre 1 et k+1=size(Cells)(2)
    mesK = Nodes(Cells(K,2)) - Nodes(Cells(K,1)); //Mesure de K, qui correspond à la dérivée de F_K, en valeur absolue
    y = mesK*sum(wwq.*f(FK(xxq)).*RefBasisFct(xxq,k)(betaa,:)); //Approximation d'intégrale
endfunction

//=========================================================================================
//Construction de la matrice
//=========================================================================================

A = zeros(Nint,Nint);

for K = 1:size(Cells)(1); //Parcourt les éléments du maillage
    for betaa = 1:size(Cells)(2); //Parcourt les noeuds de la maille K
        i = Cells(K,betaa); //Indice global du premier noeud considéré
        for alpha = 1:size(Cells)(2); //Parcourt les noeuds de la maille K une seconde fois
            j = Cells(K,alpha); //Indice global du second noeud considéré
            if (i<=Nint)&&(j<=Nint) then
                A(i,j) = A(i,j) + m(alpha,betaa,K);
            end
        end
    end
end

//=========================================================================================
//Construction du second membre
//=========================================================================================

F = zeros(Nint,1);
for K = 1:size(Cells)(1); //Parcourt les éléments du maillage
    for betaa = 1:size(Cells)(2); //Parcourt les noeuds de la maille K
        i = Cells(K,betaa); //Indice global du noeud considéré
            if (i<=Nint) then
                F(i) = F(i) + b(betaa,K);
            end
    end
end

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

U = A\F; //Routine toute faite de réolution de système linéaire

//=========================================================================================
//Calcul de la solution et se sa dérivée
//=========================================================================================

Sol = [[0;Nodes(1:Nint);1],[0;U;0]]; //On rajoute les valeurs extrêmes
Sol2 = gsort(Sol','lc','i'); //Trie selon la valeur de chaque noeud, c'est une matrice à deux lignes

X = Sol2(1,:)'; //Axe des abscisses, VECTEUR COLONNE
Ucalc = Sol2(2,:)'; //Solution calculée, VECTEUR COLONNE

//=========================================================================================
//Calcul de l'erreur en norme L2 et H1
//=========================================================================================
errL2 = 0; //Erreur en norme L2
errH1 = 0; //Erreur en norme H1
X = Sol2(1,:)'; //Axe des abscisses, VECTEUR COLONNE
nx = length(X);
select k
case 1 then
for i=1:nx-1
    function y=uh(x)
        y = ((x-X(i))./(X(i+1)-X(i)))*(Ucalc(i+1)-Ucalc(i)) + Ucalc(i);
    endfunction
    
    function y=duh(x)
        y = 0.*x + (Ucalc(i+1)-Ucalc(i))/(X(i+1)-X(i));
    endfunction
    
    function y=err(x)
        y = abs(uex(x)-uh(x))^2;
    endfunction
    
    function y=derr(x)
        y = abs(duex(x)-duh(x))^2;
    endfunction
    
    errL2 = errL2 + intg(X(i),X(i+1),err);
    errH1 = errH1 + intg(X(i),X(i+1),err)+intg(X(i),X(i+1),derr);
end
case 2 then
for i=(1:0.5*(nx-3))
    function y=uh(x)
        //On procède par interpolation polynoiale de degré 2 sur le segment
        V = [[X(2*i+1)^2,X(2*i+1),1];[X(2*i+2)^2,X(2*i+2),1];[X(2*i+3)^2,X(2*i+3),1]];
        UU = [Ucalc(2*i+1);Ucalc(2*i+2);Ucalc(2*i+3)];
        AA = V\UU;
        y = AA(1)*x^2+AA(2)*x+AA(3);
    endfunction
    
    function y=duh(x)
        //On procède par interpolation polynoiale de degré 2 sur le segment puis on dérive
        V = [[X(2*i+1)^2,X(2*i+1),1];[X(2*i+2)^2,X(2*i+2),1];[X(2*i+2)^2,X(2*i+3),1]];
        UU = [Ucalc(2*i+1);Ucalc(2*i+2);Ucalc(2*i+3)];
        AA = V\UU;
        y = 2*AA(1)*x+AA(2);
    endfunction
    
    function y=err(x)
        y = abs(uex(x)-uh(x))^2;
    endfunction
    
    function y=derr(x)
        y = abs(duex(x)-duh(x))^2;
    endfunction
    
    errL2 = errL2 + intg(X(2*i+1),X(2*i+3),err);
    errH1 = errH1 + intg(X(2*i+1),X(2*i+3),err)+intg(X(2*i+1),X(2*i+3),derr);
end    
    
end //end du select k - case


ErrL2(hh) = sqrt(errL2);
ErrH1(hh) = sqrt(errH1);





end //end de la boucle for h=H

//=========================================================================================
//Calcul du taux de convergence
//=========================================================================================
[aL2,bL2,sig]=reglin(log(H)',log(ErrL2)'); //Fonction régression linéaire,
[aH1,bH1,sig]=reglin(log(H)',log(ErrH1)'); //routine toute faite

//=========================================================================================
//Affichage des courbes de convergence
//=========================================================================================
disp("Affichage des courbes de convergence.");

f1=scf(1);

subplot(1,2,1);
set(gca(),"grid",[1 1]);
xlabel('$log(h)$');
ylabel('$log||u-u_h||_{L^2}$');
title('$\text{Courbe de convergence en norme L^2 pour les éléments finis de Lagrange}$');
plot(log(H),log(ErrL2),'sred','thickness',2);
plot(log(H),aL2*log(H)+bL2,'--blue','thickness',1);
legend('Mesures avec les élements finis IP'+string(k)+' de Lagrange','Ordre de convergence observé: '+string(aL2));

subplot(1,2,2);
set(gca(),"grid",[1 1]);
xlabel('$log(h)$');
ylabel('$log||u-u_h||_{H^1}$');
title('$\text{Courbe de convergence en norme H^1 pour les éléments finis de Lagrange}$');
plot(log(H),log(ErrH1),'sred','thickness',2);
plot(log(H),aH1*log(H)+bH1,'--blue','thickness',1);
legend('Mesures avec les élements finis IP'+string(k)+' de Lagrange','Ordre de convergence observé: '+string(aH1));

disp("C est prêt !")

