Skip to main content
Computational Science

Return to Question

Improve code
Source Link
nicoguaro ♦
nicoguaro ♦
  • 8.6k
  • 6
  • 26
  • 51
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
%%% correspondingCalculate tonormal theoutward elementpoining E1normalized vector
xi1       delta_x = t(edge_numberx_node_edge(1):,1end);  
nu1        = t(edge_number-x_node_edge(1):,21);  
P1        delta_y = elemT3_Polyy_node_edge(xi1:,nu1end); 
N_1        = P1{1,1}                ; % [N]
N_dot_1    = [P1{1-y_node_edge(:,2}; P1{1,3}]     ; % [N,xi ); N,nu]

% Jacobian matrix for the% referenceTangentiels element
J1        t1 = (N_dot_1*[x_node_1' y_node_1'])' ;
j1         =delta_x inv(J1'); % (J')^-1
N_dotp_1   = (t2 j1= *delta_y N_dot_1);  

% 
    % LegendreNormal polynomialsoutward p=1poining normalized vector 
% corresponding to the elementnv E2= [t2;-t1]/sqrt((t1^2)+(t2^2));

xi2    % Calculate the trace =: t(edge_number(2),1);Length of 
nu2the edge  
    trace_edge  = tsqrt(edge_number(2),2delta_x);
P2        ^2 =+ elemT3_Poly(xi2,nu2delta_y); 
N_2        = P2{1,1}                ; % [N]
N_dot_2    = [P2{1,2}; P2{1,3}]     ; % [N,xi ^2); N,nu]    

    % 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

%-------------------------------------------------------------------------- %--------------------------------------------------------------------------

end

end

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

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
added 4112 characters in body
Source Link
salah eddine Hebaz
salah eddine Hebaz
  • 21
  • 3
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

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

Formatting
Source Link
nicoguaro ♦
nicoguaro ♦
  • 8.6k
  • 6
  • 26
  • 51

I want to solve the 2D poisson problem using the interior penalty discontinuous galerkin methods (IPDGGalerkin 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]$

$$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  ?

Thank you

I want to solve the 2D poisson problem using the interior penalty discontinuous galerkin methods (IPDG 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  ?

Thank you

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?

Source Link
salah eddine Hebaz
salah eddine Hebaz
  • 21
  • 3