# Interior penalty discontinuous Galerkin Matlab implementation

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

• What are the problems that you are facing with triangles? – nicoguaro Jun 29 '17 at 15:49
• Did you already check this repository? – nicoguaro Jun 29 '17 at 15:59
• In fact, all of the major finite element libraries should have example programs that show how to do this in ways that are independent of whether you have quadrilaterals or triangles. – Wolfgang Bangerth Jun 29 '17 at 17:05
• Hi, Thank you for your response. I did check this matlab implementation but it does not include the IPDG methods. I tried to follow B. Riviere's book. It worked for quads but not triagles and i think i have a problem in the edges integrations. Here's the function that i've written to computes the flux terms, if you can help me correct it please. Thank you. Bests, Hebaz – salah eddine Hebaz Jul 7 '17 at 9:11

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.

• Hi , thank you very much for responding and for your help. The book is great i'm looking wright now to correct my code. And Indeed i will move to a free finite element library, it saves alot of trouble. Bests, Hebaz – salah eddine Hebaz Jul 12 '17 at 13:24
• I would like to develop a modified version of IPDG, which would modify the DG piecewise polynomial function space. Is it easy to do such change in the finite element library? Thanks! – Michael Feb 24 '19 at 22:03