function dg_poisson1D()

    dg_params = struct;
    
    % dG method parameters
    dg_params.N   = 8;    % mesh elements
    dg_params.K   = 2;    % method degree
    dg_params.Np  = 10;   % plot evaluation points
    dg_params.eta = dg_params.K*dg_params.K*dg_params.N; % dG penalization parameter

    dg_main(dg_params);

end

% Function that evaluates the basis functions and their derivatives
% in a specific point. It returns two arrays with the values.
function [phi, dphi] = basis(center, point, h, degree)
    for ii = 0:degree
        phi(ii+1) = ((point-center)/h)^ii;
        if ii == 0
            dphi(ii+1) = 0;
        else
            dphi(ii+1) = (ii/h)*((point-center)/h)^(ii-1);
        end
    end
    phi = phi';
    dphi = dphi';
end

% Function that returns the quadrature nodes and their weights
% on the 1D reference element. Uses the Golub-Welsh algorithm.
function [nodes, weights] = gauss_quadrature(doe)
    doe = doe+1;
    if rem(doe,2) == 0
        num_nodes = floor((doe+1)/2);
    else
        num_nodes = floor((doe+2)/2);
    end
    
    if num_nodes == 1
        nodes = 0;
        weights = 1;
        return;
    end
    M = zeros(num_nodes, num_nodes);
    for ii = 1:num_nodes-1
        v = sqrt(1/(4-1/((ii)^2)));
        M(ii+1,ii) = v;
        M(ii,ii+1) = v;
    end

    [d,v] = eig(M);
    nodes = diag(v)';
    weights = (d(1,:).^2);   
end

% Load function of the problem
function l = f_load(x)
    l = pi*pi*sin(pi*x);
end

function hd = howmany_dofs(K,D)
    hd = nchoosek(K+D, D);
end

function dg_main(dg_params)
    [nodes,weights] = gauss_quadrature(2*dg_params.K);
    num_local_dofs = howmany_dofs(dg_params.K, 1);
    num_total_dofs = num_local_dofs * (dg_params.N);
    
    % initialize global matrix (A) and global rhs (b)
    A = sparse(num_total_dofs, num_total_dofs);
    f = zeros(num_total_dofs,1);
    
    % compute mesh 'h'
    h = 1/dg_params.N;

    % assemble local stiffness matrices
    for ii = 1:dg_params.N
        xT = (ii-0.5) * h;
        offset = num_local_dofs * (ii-1);
        
        % map quadrature to the current element
        lweights = weights*h;
        all_ones = ones(size(nodes));
        lnodes = h*( 0.5*(nodes) + ((ii-0.5)*all_ones));
        
        % initialize local matrices
        AT = zeros(num_local_dofs, num_local_dofs);
        bT = zeros(num_local_dofs, 1);
        
        % compute local contributions
        for in = 1:length(lnodes)
            [phi, dphi] = basis(xT, lnodes(in), h, dg_params.K);    
            AT = AT + lweights(in) * dphi * dphi';
            bT = bT +  lweights(in) * f_load(lnodes(in)) * phi;
        end
        
        % copy in the global matrix
        for ib = 1:length(dphi)
            for jb = 1:length(dphi)
                A(offset+ib, offset+jb) = A(offset+ib, offset+jb) + AT(ib,jb);
            end
            f(offset+ib) = f(offset+ib) + bT(ib);
        end
    end
    
    % assemble face contributions
    for ii = 1:dg_params.N-1
        xF = ii * h;
        xT1 = (ii-0.5)*h;
        xT2 = (ii+0.5)*h;
        
        % integrating on faces in 1D means just evaluating the function
        % on a point, don't need to use quadratures
        [phiT1, dphiT1] = basis(xT1, xF, h, dg_params.K);
        [phiT2, dphiT2] = basis(xT2, xF, h, dg_params.K);
        
        % compute offsets in the global matrix
        offsetT1 = num_local_dofs * (ii-1);
        offsetT2 = num_local_dofs * ii;
        
        % compute local contributions
        AF11 = dg_params.eta * (phiT1 * phiT1');
        AF11 = AF11 - 0.5*(phiT1 * dphiT1' + dphiT1 * phiT1');
        
        AF12 = - dg_params.eta * phiT1 * phiT2';
        AF12 = AF12 - 0.5*(phiT1 * dphiT2' - dphiT1 * phiT2');
        
        AF21 = - dg_params.eta * phiT2 * phiT1';
        AF21 = AF21 - 0.5*(-phiT2 * dphiT1' + dphiT2 * phiT1');
        
        AF22 = dg_params.eta * phiT2 * phiT2';
        AF22 = AF22 - 0.5*(-phiT2 * dphiT2' - dphiT2 * phiT2');

        % copy in the global matrix
        for ib = 1:num_local_dofs
           for jb = 1:num_local_dofs
               A(offsetT1+ib, offsetT1+jb) = A(offsetT1+ib, offsetT1+jb) + AF11(ib,jb);
               A(offsetT1+ib, offsetT2+jb) = A(offsetT1+ib, offsetT2+jb) + AF12(ib,jb);
               A(offsetT2+ib, offsetT1+jb) = A(offsetT2+ib, offsetT1+jb) + AF21(ib,jb);
               A(offsetT2+ib, offsetT2+jb) = A(offsetT2+ib, offsetT2+jb) + AF22(ib,jb);
           end
        end
    end
    
    % Boundary terms
    boundary_nodes = [0,1];
    boundary_normals = [-1,1];
    boundary_elements = [1,dg_params.N];
    
    for ibf = 1:2
        xF = boundary_nodes(ibf);
        nF = boundary_normals(ibf);
        ii = boundary_elements(ibf);
        xT = (ii-0.5) * h;
        offsetT = num_local_dofs * (ii-1);
        [phi, dphi] = basis(xT, xF, h, dg_params.K);
        AF = dg_params.eta * phi * phi' - nF * (phi * dphi' + dphi * phi');
        for ib = 1:num_local_dofs
            for jb = 1:num_local_dofs
                A(offsetT+ib, offsetT+jb) = A(offsetT+ib, offsetT+jb) + AF(ib, jb);
            end
        end
    end
    
    % Solve the linear system
    u = A\f;
    
    % Postprocess
    for ii = 1:dg_params.N
        xT = (ii-0.5) * h;
        offset = num_local_dofs * ii;
        
        plotnodes = zeros(1, dg_params.Np);
        
        % Choose the plot nodes inside the element
        for ipn = 1:dg_params.Np
            plotnodes(ipn) = h * (ii-1) + ( h / (dg_params.Np - 1.) ) * (ipn-1);
        end
    
        % Extract the DoFs of the current element
        local_x = u(num_local_dofs*(ii-1)+1:num_local_dofs*(ii));
    
        % Evaluate against the basis
        for ipn = 1:dg_params.Np
            [phi, dphi] = basis(xT, plotnodes(ipn), h, dg_params.K);
            pn((ii-1) * dg_params.Np + ipn) = plotnodes(ipn);
            sol((ii-1) * dg_params.Np + ipn) = phi'*local_x;
        end;
    end

    % Plot
    plot(pn, sol);
end