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 vectordelta_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 weightsN_GAUSS_POINT = 3; %total number of gauss points = (N_GAUSS_POINT)^2GAUSS_POINT   = [-1 1 0 ]*sqrt(3/5);GAUSS_WEIGHT  = [5 5 8 ]/9;%% Calculus Loop sur Sigma de sfor ng=1:N_GAUSS_POINT%--------------------------------------------------------------------------%--------------------------------------------------------------------------    s = GAUSS_POINT(ng);% Jacobian matrix for surface integralJc     = 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 elementJ1         = (N_dot_1*[x_node_1' y_node_1'])' ;j1         = inv(J1'); % (J')^-1N_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 elementJ2         = (N_dot_2*[x_node_2' y_node_2'])';j2         = inv(J2'); % (J')^-1N_dotp_2   = ( j2 * N_dot_2) ; %--------------------------------------------------------------------------%--------------------------------------------------------------------------% Flux expressionsSigma{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 ; %BSigma{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 ; %CSigma{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%--------------------------------------------------------------------------%--------------------------------------------------------------------------endend
Interior penalty discontinuous Galerkin Matlab implementation
matlab;discontinuous galerkin