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? Commented Jun 29, 2017 at 15:49
• Did you already check this repository? Commented Jun 29, 2017 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. Commented Jun 29, 2017 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 Commented Jul 7, 2017 at 9:11

1 Answer

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 Commented Jul 12, 2017 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! Commented Feb 24, 2019 at 22:03