2D Gaussian quadratures: Quadrilateral example

Contents

Compute the 2D Gauss points on the reference element

N=2; %order of the Gaussian quadrature
[w,ptGaussRef]=gaussValues2DQuad(N);
% Draw Gauss points in the reference element
plotRectangle([-1 -1],[1,-1],[1,1],[-1,1]);
hold on;
plot(ptGaussRef(:,1),ptGaussRef(:,2),'ro');
hold off;
% 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).*(1-y)/4;
Psi2=@(x,y)(1+x).*(1-y)/4;
Psi3=@(x,y)(1+x).*(1+y)/4;
Psi4=@(x,y)(1-x).*(1+y)/4;
% Shape function derivatives
dPsi11=@(x,y) -(1-y)/4;
dPsi21=@(x,y) (1-y)/4;
dPsi31=@(x,y) (1+y)/4;
dPsi41=@(x,y) -(1+y)/4;
dPsi12=@(x,y) -(1-x)/4;
dPsi22=@(x,y) -(1+x)/4;
dPsi32=@(x,y) (1+x)/4;
dPsi42=@(x,y) (1-x)/4;
% Gradient matrix
Jacb =@(x,y) [dPsi11(x,y), dPsi21(x,y),dPsi31(x,y),dPsi41(x,y);...
              dPsi12(x,y), dPsi22(x,y),dPsi32(x,y),dPsi42(x,y)];

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

% The present function
F =@(x,y) (x.^2)- y.^2;
refElem=0; % =1 if the function is one of the shape functions (expressed in
          % in the reference element)
% The domain vertices
v1=[0,0];
v2=[5,-1];
v3=[4,5];
v4=[1,4];

Compute the corresponding Gaussian points on the domain

% evaluate Shape functions on Gaussian reference points
xx = ptGaussRef(:,1);
yy = ptGaussRef(:,2);
evalPsi1 = Psi1(xx,yy);
evalPsi2 = Psi2(xx,yy);
evalPsi3 = Psi3(xx,yy);
evalPsi4 = Psi4(xx,yy);
% according formula on slide num. 18, of document T3-MN
ptGaussDomain = evalPsi1*v1+evalPsi2*v2+evalPsi3*v3+evalPsi4*v4;
% Draw Gauss points in the present domain
figure()
plotRectangle(v1,v2,v3,v4);
hold on;
plot(ptGaussDomain(:,1),ptGaussDomain(:,2),'ro');
hold off;

Compute the Jacobian terms

vertex matrix

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

Compute the integral value according Gauss formula

%evaluate the function on the domain points
if(refElem ~= 1)
    evalF=F(ptGaussDomain(:,1),ptGaussDomain(:,2));
else
    evalF=F(ptGaussRef(:,1),ptGaussRef(:,2));
end
% Hint ---> if the function F is one the shape functions, we will use instead
% evalF=F(ptGaussRef(:,1),ptGaussRef(:,2));
% because we don't know the expression for shape functions for arbitary
% quadrilaterals but we know that they are coincidents, after change of variables,
% with the formula of the reference shape functions evaluated at ptGaussRef
%

% Finally, according to 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 =

   6.0000e+01

(c)Numerical Factory 2018


Published with MATLAB® R2017b