0
$\begingroup$

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
$\endgroup$
  • 2
    $\begingroup$ What are the problems that you are facing with triangles? $\endgroup$ – nicoguaro Jun 29 '17 at 15:49
  • $\begingroup$ Did you already check this repository? $\endgroup$ – nicoguaro Jun 29 '17 at 15:59
  • 1
    $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Jun 29 '17 at 17:05
  • $\begingroup$ 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 $\endgroup$ – salah eddine Hebaz Jul 7 '17 at 9:11
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – salah eddine Hebaz Jul 12 '17 at 13:24
  • $\begingroup$ 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! $\endgroup$ – Michael Feb 24 at 22:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.