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