2D Gaussian quadratures: Triangular example

Contents

Compute the 2D Gauss points on the reference element

N=2; %order of the Gaussian quadrature
[w,ptGaussRef]=gaussValues2DTriang(N);
% this Matlab function is defined on the slide num. 15  of document T3-MN.

Define the shape functions and their derivatives for the reference element

We use Matlab implicit function definition:

% Shape functions
Psi1=@(x,y) 1-x-y;
Psi2=@(x,y) x;
Psi3=@(x,y) y;

% Gradient matrix
Jacb =[-1, 1, 0; -1, 0,1];

Define the function $$f(x,y)$ and the domain $$\Omega_k$

% The present function
F =@(x,y) 2-x-2*y;
%F =@(x,y) (x.^2)- y.^2;
% The domain vertices
v1=[0,0];
v2=[1,1/2];
v3=[0,1];

Compute the corresponding Gaussian points on the domain

% evaluate Shape functions on Gaussian refrence points
xx = ptGaussRef(:,1);
yy = ptGaussRef(:,2);
evalPsi1 = Psi1(xx,yy);
evalPsi2 = Psi2(xx,yy);
evalPsi3 = Psi3(xx,yy);
% according formula on slide num. 18, of document T3-MN
ptGaussDomain = evalPsi1*v1+evalPsi2*v2+evalPsi3*v3;

Compute the Jacobian terms

% vertex matrix
v = [v1;v2;v3];
% evaluate Jacobian contribution for each point
for i=1:size(xx,1)
    evalDetJacb(i) = det(Jacb*v);
end

Compute the integral value according Gauss formula

%evaluate the function on the domain points
evalF=F(ptGaussDomain(:,1),ptGaussDomain(:,2));

% according formula on slide num. 20, of document T3-MN
suma=0;
for i=1:size(ptGaussDomain,1)
    suma=suma+w(i)*evalF(i)*evalDetJacb(i);
end
integ = suma

integ =

   3.3333e-01

(c)Numerical Factory 2018


Published with MATLAB® R2017b