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
//=========================================================================================

//=========================================================================================
//Saisie des données:
//      - f : Fonction "second memmbre"
//      - Lambda : Fonction "paramètre"
//      - TraceSolEx : 0 ou 1, trace la solution exacte (1) ou pas (0)
//      - k : Degré des éléments finis k, égal à 1 ou 2
//      - List : Liste des points du maillage dans ]0,1[, qui est un VECTEUR COLONNE
//=========================================================================================

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

function y=Lambda(x)
    y=0.*x+0.01;
endfunction

TraceSolEx = 1;

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

k = 1;
List = [0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9];

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

Np = length(List);

function [Nint,Next,Nodes,Cells]=createMesh(List,k)
    //Créée le maillage en fonction des points et du degré k
    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

//=========================================================================================
//A partir de maintenant, on résoud le problème
//=========================================================================================

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

disp("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]

disp('Construction des matrices et second membres élémentaires.');

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
//=========================================================================================

disp("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
//=========================================================================================

disp("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
//=========================================================================================
disp("Résolution du système linéaire.");

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

//=========================================================================================
//Affichage de la solution
//=========================================================================================
disp("Affichage de la solution.");

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

f0=scf(1);
select TraceSolEx
case 0 then
set(gca(),"grid",[1 1]);
xlabel('$x$');
ylabel('$u(x)$');
title('$\text{Solutions du problème elliptique }-u''(x) + \lambda(x) u(x)=f(x) \text{ avec }u(0)=u(1)=0$');
plot(Sol2(1,:),Sol2(2,:),'color','red','thickness',2);
legend('Approchée (IP'+string(k)+')');

case 1 then
set(gca(),"grid",[1 1]);
xlabel('$x$');
ylabel('$u(x)$');
title('$\text{Solutions du problème elliptique }-u''(x) + \lambda(x) u(x)=f(x) \text{ avec }u(0)=u(1)=0$');
plot(Sol2(1,:),Sol2(2,:),'color','red','thickness',2);
XX = 0:0.001:1;
YY = [];
for x = XX
    YY = [YY;uex(x)];
end
disp("Tracé de la solution exacte.");
plot(XX',YY,'--black','thickness',2);
legend('Approchée (IP'+string(k)+')','Exacte');

end
disp("C est prêt !")
