clear;

// Finite Elements Method in 2D-dimension space for elliptic PDE's - Delta u(x) + lambda u(x) = f(x). Delta: Laplacian operator
// Domain: 2D unit-disk


// Parameters

h = 0.1 // Step size for numerical integration
bound = "Dirichlet"  // Boundary conditions (Dirichlet or Neumann)
lambda = 0

disp("===================================")
disp("2D Finite Element Method")
disp(["Boundary condition:", bound])
disp("===================================")

// Functions and vectors

function z=Test(x)
    z = 1-x(1)-x(2) + (1e-12)
endfunction

function z=f(x)
    // Second member function
    // Inputs:
    // - x: Space variable, size: (2,1)
    z = 4 + (1e-12)*x(1)
endfunction


function z=u_ex(x)
    // Exact solution
    
    // Inputs:
    // - x: Space variable, size: (2,1)
    z = 1 - x(1).^2 - x(2).^2
endfunction




function z=grad(g,x)
    // Gradient of a function with evaluation at point x
    // Inputs:
    // - g: Function with 2D input variable
    // - x: Space variable (evaluation point), size: (2,1)
    z(1) = (g([x(1)+h,x(2)]) - g([x(1)-h,x(2)]))/(2*h)
    z(2) = (g([x(1),x(2)+h]) - g([x(1),x(2)-h]))/(2*h)
endfunction

function I = Int_Delta(g, X1, X2, X3)
    // Integration of function g on triangle with verticies X1, X2 and X3. Numerical approximation with existing numerical integration
    // Inputs:
    // - g: Function with 2D input variable
    // - X1, X2, X3: Vertices of the triangle. Size: (1,2)

    function zz = gg(x1,x2)
        // Adapted function for numerical integration
        x = [x1, x2]
        zz = g(x)
    endfunction
    
    I = int2d([X1(1);X2(1);X3(1)], [X1(2);X2(2);X3(2)], gg)
endfunction

function aaa = aa(u, v, X1, X2, X3)
    // Bilinear form in the weak formulation of the PDE on the triangle with vertices X1, X2 and X3
    // Inputs:
    // - u, v: Functions with 2D inputs and scalar outputs
    // - X1, X2, X3: Vertices of the triangle. Shape: (2,1)
    function z = g(x)
        // Function which is integrated on the triangle
        //z = (u([x(1)+h,x(2)]) - u([x(1)-h,x(2)]))/(2*h)
        z = sum(grad(u,x).*grad(v,x)) + lambda*u(x)*v(x)
    endfunction
    aaa = Int_Delta(g, X1, X2, X3)
endfunction

function LLL = LL(v, X1, X2, X3)
    // Linear form in the weak formulation of the PDE on the triangle with vertices X1, X2 and X3
    // Inputs:
    // - v: Function with 2D input and scalar output
    // - X1, X2, X3: Vertices of the triangle. Shape: (2,1)
    function z = g(x)
        // Function which is integrated on the triangle
        z = f(x)*v(x)
    endfunction
    LLL = Int_Delta(g, X1, X2, X3)
endfunction








function [C,V] = FE(J,K)
    // Gives mesh (cells and vertices) for finite elements
    // Inputs:
    // - J: Number of verticies at the edge of the domain
    // - K: Number of layers of points inside the unit disk
    V = zeros(J*K+1, 3)
    for k=0:K-1
        J0 = (1-(-1)^k)/4 + int(k/2)
        V(1+k*J:1+(k+1)*J-1,1) =  (1-k/K)*cos(2*%pi*(J0:1:J+J0-1)/J)'
        V(1+k*J:1+(k+1)*J-1,2) =  (1-k/K)*sin(2*%pi*(J0:1:J+J0-1)/J)'
    end

    V((J)*K+1,1) = 0
    V((J)*K+1,2) = 0
    
    V(1:J,3) = ones(J,1) // 1 for belong to bound
    
    
    C1 = []
    C2 = []
    C3 = []
    for k = 0:K-2
        C1 = [C1 ; [(J*k+1:J*k+J)']]
        C2 = [C2 ; [(J*k+2:J*k+J)';J*k+1]]
        C3 = [C3 ; [(J*k+J+1:J*k+2*J)']]

        C1 = [C1 ; [(J*k+J+1:J*k+2*J)']]
        C2 = [C2 ; [(J*k+J+2:J*k+2*J)';J*k+J+1]]
        C3 = [C3 ; [(J*k+2:J*k+J)';J*k+1]]
    end
    
    C1 = [C1 ; [(J*K-J+1:J*K)']]
    C2 = [C2 ; [(J*K-J+2:J*K)';J*K-J+1]]
    C3 = [C3 ; [(J*K+1)*ones(J,1)]]
    
    C = [C1,C2,C3]

endfunction



       
function [h_mesh, L2_error, H1_error] = FEM_Error(J, K)
    // Finite Elements Method to solve elliptic PDE on the unit disk in two-dimensional space, with homogeneous Dirichlet or Neumann boundary conditions
    // Inputs:
    // - J: Number of verticies at the edge of the domain
    // - K: Number of layers of points inside the unit disk
    // Returns biggest Lebesgue measure cell of the mesh and L2/H1 errors between exact and approximated solution.
    
    disp("J = "+string(J)+" - K = "+string(K))
    
    if bound == "Dirichlet" then
        A = zeros(J*K+1-J,J*K+1-J)
        B = zeros(J*K+1-J,1)
    elseif bound == "Neumann" then
        A = zeros(J*K+1,J*K+1)
        B = zeros(J*K+1,1)
    end
    
    
    //disp('Mesh creation.')
    [C, V] = FE(J,K)
    
    //disp('Matrices constructions.')
    for k = 1:size(C)(1)
        //disp("Loading [%]")
        //disp(int(1000*k/size(C)(1))/10)
        X1 = V(C(k,1),1:2)
        X2 = V(C(k,2),1:2)
        X3 = V(C(k,3),1:2)
        M = [[1,X1];[1,X2];[1,X3]]
        
        for betaa = 1:3
            i = C(k, betaa)
            e_beta = eye(3,3)(:,betaa)
            a_beta = M\e_beta
            
            function z = phi_beta(x)
                // Function of local basis with index alpha
                // Inputs
                // - x: Space variable. Shape: (2,1)
                z = a_beta(1) + a_beta(2)*x(1) + a_beta(3)*x(2)
            endfunction
            
            if bound == "Dirichlet" then
                if (V(i,3) == 0) then
                    B(i-J) = B(i-J) + LL(phi_beta, X1, X2, X3)
                end
            elseif bound == "Neumann" then
                B(i) = B(i) + LL(phi_beta, X1, X2, X3)
            end   
                
            for alphaa = 1:3 
                j = C(k, alphaa)
                e_alpha = eye(3,3)(:,alphaa)
                a_alpha = M\e_alpha      
                 
                function z = phi_alpha(x)
                    // Function of local basis with index alpha
                    // Inputs
                    // - x: Space variable. Shape: (2,1)
                    z = a_alpha(1) + a_alpha(2)*x(1) + a_alpha(3)*x(2)
                endfunction
                               
                if bound == "Dirichlet" then
                    if (V(i,3) == 0) && (V(j,3) == 0) then
                        A(i-J,j-J) = A(i-J,j-J) + aa(phi_alpha, phi_beta, X1, X2, X3)
                    end
                elseif bound == "Neumann" then
                    A(i,j) = A(i,j) + aa(phi_alpha, phi_beta, X1, X2, X3)
                end             
            end
        end
    end
    
    U = A\B
    
    if bound == "Dirichlet" then
        U = [zeros(J,1);U]
    end
    
    L2_error = 0
    H1_error = 0
    h_mesh = 0
    
    for k = 1:size(C)(1)
        X1 = V(C(k,1),1:2)
        X2 = V(C(k,2),1:2)
        X3 = V(C(k,3),1:2)
        M = [[1,X1];[1,X2];[1,X3]]
        a = M\[U(C(k,1));U(C(k,2));U(C(k,3))]
        
        funcprot(0)
        function z = diff(x)
            z = (a(1) + a(2)*x(1) + a(3)*x(2) - u_ex(x))^2
        endfunction
        
        function z = diff_grad(x)
            //disp("diff:",a(2:3)'-2*x)
            z = (a(2) + 2*x(1))^2 + (a(3) + 2*x(2))^2
        endfunction
        
        //hh = (1/2)*abs(det([X2-X1;X3-X1]))
        hh = max(norm(X1-X2),norm(X2-X3),norm(X3-X1))
        
        if hh > h_mesh then
            h_mesh = hh
        end
        
        Int_diff = Int_Delta(diff, X1, X2, X3)
        Int_diff_grad = Int_Delta(diff_grad, X1, X2, X3)
        
        //disp("k:", k, "diff:",Int_diff, "diff_grad", Int_diff_grad)
        
        L2_error = L2_error + Int_diff
        H1_error = H1_error + Int_diff + Int_diff_grad
    end
    
    L2_error = sqrt(L2_error)
    H1_error = sqrt(H1_error)
    
endfunction

function Plot_Convergence(save_fig)
    // Plots a convergence curve of FEM method.
    // Inputs:
    // - save_fig: Saves the figure or not. Boolean
    
    NN1 = 3; NN2 = 12; stride = 1
    NN = flipdim(NN1:stride:NN2,2)
    H_MESH = []; L2_ERROR = []; H1_ERROR = []
    
    for n = NN
        [h_mesh, L2_error, H1_error] = FEM_Error(2*n, n)
        H_MESH = [H_MESH, h_mesh]
        L2_ERROR = [L2_ERROR, L2_error]
        H1_ERROR = [H1_ERROR, H1_error]
    end
    
    [ord_L2,ord2_L2] = reglin(log(H_MESH),log(L2_ERROR))
    [ord_H1,ord2_H1] = reglin(log(H_MESH),log(H1_ERROR))
    
    
    scf()
    plot("ll",H_MESH,L2_ERROR,"rs",'thickness',2)
    plot("ll",H_MESH,H1_ERROR,"gs",'thickness',2)
    xlabel("$\delta(T)$")
    ylabel("Error")
    legend("$|| u_h - u ||_{L^2(\Omega)}$", "$|| u_h - u ||_{H^1(\Omega)}$")
    title("Convergence curve - order [L2]:"+string(int(100*ord_L2)/100)+" - order [H1]:"+string(int(100*ord_H1)/100))
    gca().grid=[1 1 1]
    
    if save_fig == %T then
        filename="C:\Users\Utilisateur\Documents\Scripts_Scilab\Finite_Elements_2D\Figures_Convergence\Finite_Element_Method_2D_NN1="+string(NN1)+"_NN2="+string(NN2)+"_stride="+string(stride)+"_bound="+bound;
        xs2pdf(0,filename);
        xs2pdf(gcf(),filename);
    end
    
endfunction




