% ----------------------------------------------------------------------- %
% Methodes de discretisation en fluides                                   %
% Laboratoire I - Code de départ                                          %
% M.A. Habisreutinger                                                     %
% Jeudi 14 mars 2024                                                      %
% ----------------------------------------------------------------------- %

function Lab01()

    %  1.   Parameters

    %  1.1  Domain
    limits = [-1, +1] ;
    
    %  1.2  Mesh
    npoints = 64 ;

    %  1.3  Diffusion coefficient
    nu = 1e+00 ;


    %  2.  Computation

    %  2.1  Mesh points
    x = linspace(limits(1), limits(2), npoints)' ;

    %  2.2  Mesh normals
    normals = zeros(npoints, 1) ;
    normals(1) = -1 ;
    normals(end) = +1 ;
    
    %  2.3  Mesh size
    h = (limits(2)-limits(1))/(npoints-1) ;

    %  2.4  Analytical solution
    k = 2*pi/(limits(2)-limits(1)) ;
    u = cos(k*x+pi/4) ;
    u_x = -k*sin(k*x+pi/4) ;

    %  2.5  Mass matrix (colocation method)
    e = ones(npoints, 1) ;
    M = spdiags(e, 0, npoints, npoints) ;
    
    %  Correction for boundary conditions
    M(1, 1) = 0 ;
    M(end, end) = 0 ;

    %  2.6  Laplacian matrix (central finite differences - 2nd order)
    laplacian = DifferentiationMatrix(2, -1:1, npoints, h) ;
    %         = (1/h^2)*spdiags([+1*e, -2*e, +1*e], -1:1, npoints, npoints) ;
    
    %  Correction for Neumann boundary conditions
    taylor = TaylorCoefficients(0:+2, h) ;
    laplacian(1, 1) = +taylor(2, 1) ; % = +3/(2*h)
    laplacian(1, 2) = +taylor(2, 2) ; % = -4/(2*h)
    laplacian(1, 3) = +taylor(2, 3) ; % = +1/(2*h)

    taylor = TaylorCoefficients(-2:0, h) ;
    laplacian(end, end-2) = -taylor(2, 1) ; % = -1/(2*h)
    laplacian(end, end-1) = -taylor(2, 2) ; % = +4/(2*h)
    laplacian(end, end)   = -taylor(2, 3) ;   % = -3/(2*h)
        
    %  2.7  Discretization matrix
    A = -nu*laplacian ;
    
    %  2.8  Forcing term
    forcing = M*norm(k)^2*nu*u ;
    
    %  2.9  Neumann boundary conditions
    neumann = x == limits(2) ;
    forcing = forcing+neumann.*(nu*u_x.*normals) ;
    
    %  2.10 Dirichlet boundary conditions
    dirichlet = x == limits(1) ;
    
    Dh = distribution(dirichlet, true) ;
    Lh = distribution(dirichlet, false) ;
    
    %  2.11 System of equations
    Ah = Lh*A*Lh+Dh ;
    bh = Lh*(forcing-A*Dh*u)+Dh*u ;

    %  2.12 Numerical solution (Ah*uh = bh)
    uh = Ah\bh ;


    %  3.   Data analysis

    %  3.1  Error norms
    eh = norm(M*(uh-u))/norm(M*u) ;
    disp(['Error norm : ', num2str(eh)])

    %  3.2  Solutions plot
    figure(1) ; set(gca, 'FontName', 'Lucida Sans', 'FontSize', 14) ;
    plot(x, uh, 'bo-', x, u, 'k-') ; 
    xlabel('x') ;
    ylabel('u') ;
    title('Solutions of the problem') ;
    legend('Numerical solution', 'Exact solution') ;
    grid on ;
    axis([limits(1), limits(2), -1, +1]) ;

    %  3.3  Matrix plot
    figure(2) ;
    spy(Ah) ;
    xlabel('column index') ;
    ylabel('row index') ; 
    title('Structure of the disretisation matrix') ;

end

function w = distribution(set, mode)

    if mode
        w = spdiags(set, 0, numel(set), numel(set)) ;
    else
        w = spdiags(not(set), 0, numel(set), numel(set)) ;
    end

end

function w = TaylorCoefficients(stencil, h)
    
    nstencil = numel(stencil) ;
    T = zeros(nstencil, nstencil) ;
    
    for ip = 0 : nstencil-1
        for jp = 0 : nstencil-1 
            T(ip+1, jp+1) = (stencil(ip+1)*h)^jp/factorial(jp) ;
        end
    end
    
    w = inv(T) ;
    
end

function w = DifferentiationMatrix(order, stencil, npoints, h)

    e = ones(npoints, 1) ;
    taylor = TaylorCoefficients(stencil, h) ;
    
    diagonals = zeros(npoints, numel(stencil)) ;
    for idiag = 1 : numel(stencil)
        diagonals(:, idiag) = e*taylor(order+1, idiag) ;
    end
    
    w = spdiags(diagonals, stencil, npoints, npoints) ;
    
end