Interior penalty discontinuous Galerkin Matlab implementation

Question

I want to solve the 2D poisson problem using the interior penalty discontinuous Galerkin methods:

−∇a(x)(∇u)=0 in Ω.

The variational formulation is such that :

$$a_{\epsilon}(u,v)=\sum_{K\in T_h}\int_K a\nabla u\cdot \nabla v-\sum_{e\in \Gamma_h\cup \Gamma_D}\int_e\{a\nabla u\cdot n_e\}[v]+ \sum_{e\in \Gamma_h\cup \Gamma_D}\int_e\epsilon\{a\nabla v\cdot n_e\}[u] + \sum_{e\in \Gamma_h\cup \Gamma_D}\frac{\sigma_0}{|e|}\int_e[u][v]$$

I managed to get the implementation right and to compute and assemble the flux terms for a quadrilateral element. Yet for the triangles, it gets a little bit hard and I dont know how to do it. Can someone help me please?

function Sigma = edgeT3_Lagrange(x_node_edge,y_node_edge,eps,alp,gam,bet,edge_data)

    %% Calculate normal outward poining normalized vector
    delta_x = x_node_edge(:,end)-x_node_edge(:,1); 
    delta_y = y_node_edge(:,end)-y_node_edge(:,1); 

    % Tangentiels 
    t1 = delta_x ; 
    t2 = delta_y ;  

    % Normal outward poining normalized vector 
    nv = [t2;-t1]/sqrt((t1^2)+(t2^2));

    % Calculate the trace : Length of the edge 
    trace_edge  = sqrt((delta_x)^2 + (delta_y)^2);

    % the average local mesh size 
    h_avg  = (h_E1 + h_E2)/2 ; /circumradius 
    h_avg_beta  = h_avg.^bet ; 

    %% Space allocation
    % Matrices Initiation 
    Sigma = [{zeros(DG_DATA.NN_elem)} {zeros(DG_DATA.NN_elem)} {zeros(DG_DATA.NN_elem)} {zeros(DG_DATA.NN_elem)}];

    %% Gauss points and weights
    N_GAUSS_POINT = 3; %total number of gauss points = (N_GAUSS_POINT)^2
    GAUSS_POINT   = [-1 1 0 ]*sqrt(3/5);
    GAUSS_WEIGHT  = [5 5 8 ]/9;

    %% Calculus Loop sur Sigma de s
    for ng=1:N_GAUSS_POINT

        %-----------------------------------------------------------------------
        %-----------------------------------------------------------------------
            s = GAUSS_POINT(ng);

            % Jacobian matrix for surface integral
            Jc     = trace_edge/2 ; 
            detJc  = det(Jc); 

        %-----------------------------------------------------------------------
        %-----------------------------------------------------------------------
        % Assignment of the edge
        %-------------------------------------
        %-------------------------------------
        % For     :      the edge {i} of E1 
        % Thereby >>     the edge {j} of E2 
        % And the integration is done on s
        %-------------------------------------
        %-------------------------------------

            S = (s+1)/2;    

        %   t  = [xi  nu]
            t  = [S        0;
                  1-S      S;
                   0       S;];

        %-----------------------------------------------------------------------
        %-----------------------------------------------------------------------

            % Legendre polynomials p=1
            % corresponding to the element E1 
            xi1        = t(edge_number(1),1);  
            nu1        = t(edge_number(1),2);  
            P1         = elemT3_Poly(xi1,nu1); 
            N_1        = P1{1,1}                ; % [N]
            N_dot_1    = [P1{1,2}; P1{1,3}]     ; % [N,xi ; N,nu]

            % Jacobian matrix for the reference element
            J1         = (N_dot_1*[x_node_1' y_node_1'])' ;
            j1         = inv(J1'); % (J')^-1
            N_dotp_1   = ( j1 * N_dot_1); 

            % 
            % Legendre polynomials p=1
            % corresponding to the element E2 

            xi2        = t(edge_number(2),1);  
            nu2        = t(edge_number(2),2);
            P2         = elemT3_Poly(xi2,nu2); 
            N_2        = P2{1,1}                ; % [N]
            N_dot_2    = [P2{1,2}; P2{1,3}]     ; % [N,xi ; N,nu]    

            % Jacobian matrix for the reference element
            J2         = (N_dot_2*[x_node_2' y_node_2'])';
            j2         = inv(J2'); % (J')^-1
            N_dotp_2   = ( j2 * N_dot_2) ; 

        %-----------------------------------------------------------------------
        %-----------------------------------------------------------------------

            % Flux expressions
            Sigma{1} = Sigma{1} + GAUSS_WEIGHT(ng)*( -0.5*N_1.'*(nv'*(N_dotp_1)) + (eps/2)*((N_dotp_1).'*nv)*N_1 + (alp/h_avg_beta)*(N_1.')*N_1 + (gam/h_avg_beta)*(N_dotp_1.')*N_dotp_1 )*detJc ; %B
            Sigma{2} = Sigma{2} + GAUSS_WEIGHT(ng)*( +0.5*N_2.'*(nv'*(N_dotp_2)) - (eps/2)*((N_dotp_2).'*nv)*N_2 + (alp/h_avg_beta)*(N_2.')*N_2 + (gam/h_avg_beta)*(N_dotp_2.')*N_dotp_2 )*detJc ; %C
            Sigma{3} = Sigma{3} + GAUSS_WEIGHT(ng)*( +0.5*N_2.'*(nv'*(N_dotp_1)) + (eps/2)*((N_dotp_2).'*nv)*N_1 - (alp/h_avg_beta)*(N_2.')*N_1 - (gam/h_avg_beta)*(N_dotp_2.')*N_dotp_1 )*detJc ; %D    
            Sigma{4} = Sigma{4} + GAUSS_WEIGHT(ng)*( -0.5*N_1.'*(nv'*(N_dotp_2)) - (eps/2)*((N_dotp_1).'*nv)*N_2 - (alp/h_avg_beta)*(N_1.')*N_2 - (gam/h_avg_beta)*(N_dotp_1.')*N_dotp_2 )*detJc ; %E

        %-----------------------------------------------------------------------
        %-----------------------------------------------------------------------

    end

end

Answer 1

If you are adamant on using MATLAB, chapter 14 of the following book walks through a 2D IPDG Poisson problem using piecewise linear basis functions on triangles:

The Finite Element Method: Theory, Implementation, and Applications

The MATLAB code used for the most examples in the book is freely available on github:

https://github.com/Jumziey/FEM/tree/master/MatlabCodeFEMBook

The IPDG specific code is in the script "dG1PoissonSolver2D.m "

As mentioned in the comments, it might be worth while to look into a finite element library, as you can probably have a much easier time extending your problem (3D, different elements, adaptivity, high order, and so on) than in your own MATLAB codes.