
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
//Application du raffinement du maillage
//=========================================================================================

//=========================================================================================
//Saisie des données:
//      - f : Fonction "second memmbre"
//      - Lambda : Fonction "paramètre"
//      - k : Degré des éléments finis k, égal à 1 ou 2
//      - N : Nombre de listes de points
//      - List1,...,ListN : Liste des points du maillage dans ]0,1[, qui sont des
//VECTEURS COLONNE, List1 étant le moins raffiné, ListN le plus raffiné, N est le nombre de
//raffinements
//=========================================================================================

function y=f(x)
    if (x<0.5) then
        y = sin(1./x);
    else
        y=0;
    end
endfunction

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

k = 1;

N = 7;

List1 = [0.2;0.4;0.6;0.8];
List2 = [0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9];
List3 = [(0.05:0.05:0.45)';(0.5:0.1:0.9)'];
List4 = [(0.025:0.025:0.475)';(0.5:0.1:0.9)'];
List5 = [(0.0125:0.0125:0.2375)';(0.25:0.025:0.475)';(0.5:0.1:0.9)'];
List5 = [(0.00625:0.00625:0.24375)';(0.25:0.025:0.475)';(0.5:0.1:0.9)'];
List6 = [(0.003125:0.003125:0.246875)';(0.25:0.025:0.475)';(0.5:0.1:0.9)'];
List7 = [(0.0015625:0.0015625:0.2484375)';(0.25:0.025:0.475)';(0.5:0.1:0.9)'];

//=========================================================================================
//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 de la liste des listes
//=========================================================================================

LIST = [List1;List2;List3;List4;List5;List6;List7];
Long = [0;length(List1);length(List2);length(List3);length(List4);length(List5);length(List6);length(List7)];

for jj=1:N

imin = sum(Long(1:jj))+1;
imax = sum(Long(1:jj+1));    
List = LIST(imin:imax);

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

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

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

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);
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 \text{ avec maillage raffiné}$');
plot(Sol2(1,:),Sol2(2,:),'color',color(255-35*(jj-1),0,35*(jj-1)),'thickness',1);

end //end de la boucle "for" balayant les listes de maillages raffinés

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