clear;
clf;

parm_N = 10;
parm_K = 3;
parm_h = 1/parm_N;

dg_num_local_dofs = parm_K+1;
dg_num_total_dofs = parm_N * dg_num_local_dofs;
dg_num_plot_nodes = parm_K+1;

parm_eta = 2;
parm_beta = 0.2;

function l = f_load(x)
    l = sin(%pi*x) * %pi^2;
    return;
endfunction

function l = u(x)
    l = sin(%pi*x);
    return;
endfunction

function [nodes, weights] = gauss_quadrature(doe)
    if modulo(doe,2) == 0 then
        num_nodes = int((doe+1)/2);
        doe = doe+1;
    else
        num_nodes = int((doe+2)/2);
    end

    if num_nodes == 1 then
        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] = spec(M);
    nodes = diag(v)';
    weights = (d(1,:).^2);   
    return;
endfunction

function [phi, dphi] = basis(centr, point, h, dgree)
    phi=zeros(dgree+1,1)
    dphi=zeros(dgree+1,1)
    
       phi(1)=1
       dphi(1)=0
    for i=2:dgree+1
       phi(i) = ((point-centr)/h)^(i-1)
       dphi(i) = (i-1)*(1/h)*((point-centr)/h)^(i-2)
    end
    
endfunction



function dg_main()
    
    [nodes,weights] = gauss_quadrature(2*parm_K);
    A = zeros(dg_num_total_dofs,dg_num_total_dofs);
    b = zeros(dg_num_total_dofs,1);
    
    for ii = 1:parm_N
        xT = (ii-0.5) * parm_h;
        offset = dg_num_local_dofs * (ii-1);
        
        lweights = weights*parm_h;
        all_ones = ones(length(nodes));
        lnodes = parm_h*( 0.5*(nodes) + ((ii-0.5)*all_ones));
        
        bs = parm_K+1;
        
        AT = zeros(bs, bs);
        bT = zeros(bs, 1);
        
        for in = 1:length(lnodes)
            [phi, dphi] = basis(xT, lnodes(in), parm_h, parm_K);

            for ib = 1:length(dphi)
                for jb = 1:length(dphi)
                    AT(ib, jb) = AT(ib, jb) + lweights(in) * dphi(ib) * dphi(jb);
                    AT(ib, jb) = AT(ib, jb) + parm_beta* phi(ib) * dphi(jb);
                end
                bT(ib) = bT(ib) +  lweights(in) * f_load(lnodes(in)) * phi(ib);
            end
        end
        
        for ib = 1:length(dphi)
            for jb = 1:length(dphi)
                A(offset+ib, offset+jb) = A(offset+ib, offset+jb) + AT(ib,jb);
            end
            b(offset+ib) = b(offset+ib) + bT(ib);
        end
        
    end
    
    for ii = 1:parm_N-1
        xF = ii * parm_h;
        xT1 = (ii-0.5)*parm_h;
        xT2 = (ii+0.5)*parm_h;
        
        [phiT1, dphiT1] = basis(xT1, xF, parm_h, parm_K);
        [phiT2, dphiT2] = basis(xT2, xF, parm_h, parm_K);
        
        offsetT1 = dg_num_local_dofs * (ii-1);
        offsetT2 = dg_num_local_dofs * ii;

        AF11 = parm_eta * phiT1 * phiT1';
        AF11 = AF11 - 0.5*(phiT1 * dphiT1' - dphiT1 * phiT1');
        
        AF12 = - parm_eta * phiT1 * phiT2';
        AF12 = AF12 - 0.5*(phiT1 * dphiT2' + dphiT1 * phiT2');
        
        AF21 = - parm_eta * phiT2 * phiT1';
        AF21 = AF21 - 0.5*(-phiT2 * dphiT1' - dphiT2 * phiT1');
        
        AF22 = parm_eta * phiT2 * phiT2';
        AF22 = AF22 - 0.5*(-phiT2 * dphiT2' + dphiT2 * phiT2');
        
        // Advection terms
        AF11 = AF11 - parm_beta/2*phiT1 * phiT1' + abs(parm_beta)/2*phiT1 * phiT1'
        AF12 = AF12 + parm_beta/2*phiT1 * phiT2' - abs(parm_beta)/2*phiT1 * phiT2'
        AF21 = AF21 - parm_beta/2*phiT2 * phiT1' - abs(parm_beta)/2*phiT2 * phiT1'
        AF22 = AF22 + parm_beta/2*phiT2 * phiT2' + abs(parm_beta)/2*phiT2 * phiT2'
        
         for ib = 1:dg_num_local_dofs
            for jb = 1:dg_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_nodes = [0,1];
    boundary_normals = [-1,1];
    boundary_elements = [1,parm_N];
    
    for ibf = 1:2
        xF = boundary_nodes(ibf);
        nF = boundary_normals(ibf);
        ii = boundary_elements(ibf);
        xT = (ii-0.5) * parm_h;
        offsetT = dg_num_local_dofs * (ii-1);
        [phi, dphi] = basis(xT, xF, parm_h, parm_K);
        
        AF = parm_eta*phi*phi'+nF*(dphi*phi'-phi*dphi') 
        
        // Advection term
        AF=AF +  1/2*(abs(parm_beta)-parm_beta)*phi*phi'

        for ib = 1:dg_num_local_dofs
            for jb = 1:dg_num_local_dofs
                A(offsetT+ib, offsetT+jb) = A(offsetT+ib, offsetT+jb) + AF(ib, jb);
            end
        end
    end
    
    x = A\b;
    
    for ii = 1:parm_N
        xT = (ii-0.5) * parm_h;
        offset = dg_num_local_dofs * ii;
        
        plotnodes = zeros(1, dg_num_plot_nodes);
        
        for ipn = 1:dg_num_plot_nodes
            plotnodes(ipn) = parm_h * (ii-1) + ( parm_h / (dg_num_plot_nodes - 1.) ) * (ipn-1);
        end
    
        local_x = x(dg_num_local_dofs*(ii-1)+1:dg_num_local_dofs*(ii));
    
        for ipn = 1:dg_num_plot_nodes
            [phi, dphi] = basis(xT, plotnodes(ipn), parm_h, parm_K);
            pn((ii-1) * dg_num_plot_nodes + ipn) = plotnodes(ipn);
            sol((ii-1) * dg_num_plot_nodes + ipn) = phi'*local_x;
        end;
    end

    disp(pn);
    disp(sol);

    plot(pn, sol);
    
endfunction

dg_main()
