2D Integrals : Gaussian Quadrilateral Example

Table of Contents

A first example on the reference quadrilateral Compute the 2D Gauss points on the reference element General case: A quadrilateral domain Example: Compute the corresponding Gaussian points on the domain Compute the Jacobian terms Compute the integral value according Gauss formula Exercise 1: Build the integQuad function Application: Integration over a mesh Structure of the Mesh Additional Functions Gauss points in 2D Quadrilateral integration function: integQuad Exercise 2: numerical integration error

A first example on the reference quadrilateral

Let's consider the function

defined on the reference quadrilateral

with the vertices

We want to compute numerically the double integral

by using the 2D Gauss formula which is similar to the 1D case but it involves the 2D Gauss points

with weights

. where N is the total number of Gauss points.

Compute the 2D Gauss points on the reference element

First we compute the appropriate Gauss points in the reference quadrilateral. They are obtained from the 1D Gauss points in

in both axes x and y and are generated by the Matlab function gaussValues2DQuad which can be found at the end of this script.

Let N denote the number of 1D Gauss points in each of the axes x,y. Then the number of 2D Gauss points in

is N^2. Notice that if

is a polynomial function of degree n < (2·N+1) in each of the variables, then the 2D Gauss integration formula returns an exact value of the integral. Since in our case

, to get an exact result it is sufficient to consider the case

N = 2; %order of the Gaussian 1D quadrature in each axis
[w,ptGaussRef] = gaussValues2DQuad(N);
fprintf('N = %d, w= %f %f %f %f\n', 2*N,w); 
N = 4, w= 1.000000 1.000000 1.000000 1.000000
ptGaussRef % exact values +/- (1/sqrt(3))
ptGaussRef = 4×2
   -0.5774   -0.5774
   -0.5774    0.5774
    0.5774   -0.5774
    0.5774    0.5774

grafical representation

% Display Gauss points in the reference quadrilateral [-1,1]x[-1,1]

vertices = [-1 -1; 1,-1; 1,1; -1,1; -1 -1];

plot(vertices(:,1), vertices(:,2));

hold on;

plot(ptGaussRef(:,1),ptGaussRef(:,2),'ro');

hold off;

F = @(x,y) (x.^2).*y.^2;

Fxy = F(ptGaussRef(:,1),ptGaussRef(:,2));

Integ = w*Fxy

Integ = 0.4444

General case: A quadrilateral domain

Let's consider the function

and the quadrilateral

with the vertices

We want to compute the following integral by using the formula of change of variables

where

is the absolute value of the determinat of the Jacobian matrix of the change of coordinates.

Below, you have an slide showing the change of variables needed to relate the reference quadrilateral

with a general one. The

are known as shape functions associated to each vertice

and are explained below:

where

are functions associated to the respective vertex of the refrence quadrilateral and satisfy

and

Example:

Now consider our function

and we code in Matlab the above change of variable formulas.

% The present function
F = @(x,y) (x.^2)- y.^2;
% The domain vertices
v1 = [0,0];
v2 = [5,-1];
v3 = [4,5];
v4 = [1,4];
% We use Matlab *implicit function* definition:
%
% Shape functions
Psi1 = @(xi,eta)(1-xi).*(1-eta)/4;
Psi2 = @(xi,eta)(1+xi).*(1-eta)/4;
Psi3 = @(xi,eta)(1+xi).*(1+eta)/4;
Psi4 = @(xi,eta)(1-xi).*(1+eta)/4;
% Shape function derivatives
dPsi11 = @(xi,eta) -(1-eta)/4;
dPsi21 = @(xi,eta) (1-eta)/4;
dPsi31 = @(xi,eta) (1+eta)/4;
dPsi41 = @(xi,eta) -(1+eta)/4;
dPsi12 = @(xi,eta) -(1-xi)/4;
dPsi22 = @(xi,eta) -(1+xi)/4;
dPsi32 = @(xi,eta) (1+xi)/4;
dPsi42 = @(xi,eta) (1-xi)/4;
% Gradient matrix
Jacb = @(xi,eta) [dPsi11(xi,eta), dPsi21(xi,eta),dPsi31(xi,eta),dPsi41(xi,eta);...
              dPsi12(xi,eta), dPsi22(xi,eta),dPsi32(xi,eta),dPsi42(xi,eta)];

Compute the corresponding Gaussian points on the domain

As we know from above, the change of variables from the reference quadrilateral to a general one is:

evaluate Shape functions on Gaussian reference points

xi = ptGaussRef(:,1);

eta = ptGaussRef(:,2);

evalPsi1 = Psi1(xi,eta);

evalPsi2 = Psi2(xi,eta);

evalPsi3 = Psi3(xi,eta);

evalPsi4 = Psi4(xi,eta);

% change these points to the new quadrilateral(from the change of variables function)

ptGaussDomain = evalPsi1*v1+evalPsi2*v2+evalPsi3*v3+evalPsi4*v4;

% Draw Gauss points in the present domain

figure()

vertices = [v1;v2;v3;v4;v1];

plot(vertices(:,1), vertices(:,2));

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 Gauss point 
for i = 1:size(xi,1)
    evalDetJacb(i) = abs(det(Jacb(xi(i),eta(i))*v));
end

Compute the integral value according Gauss formula

%evaluate the function on the domain points
evalF = F(ptGaussDomain(:,1),ptGaussDomain(:,2));
 
% Finally, apply Gauss integration formula
suma = 0;
for i = 1:size(ptGaussDomain,1)
    suma = suma+w(i)*evalF(i)*evalDetJacb(i);
end
integ = suma
integ = 60.0000

Warning 1. As we know from the theory, the area of any elementary domain in the plane (x,y) and, in particular, of a quadrilateral

, can be calculated via the double integral

Note that if we apply the above MATLAB code to this case (or the case of any constant function) by writing the subintegral function

F = @(x,y) 1
F = function_handle with value:
    @(x,y)1

F = @(x,y) 1.
F = function_handle with value:
    @(x,y)1. 

the program will return an error message: "Index exceeds the number of array elements".

One of the simplest solutions would be to write instead, for example

F = @(x,y) ones(size(x));
%or
F = @(x,y) x./x 
F = function_handle with value:
    @(x,y)x./x

but then we must ensure that the function F will not be evaluated at points with zero x-coordinate.

The definitive solution is

F = @(x,y) (x.^2+y.^2+1)./(x.^2+y.^2+1)
F = function_handle with value:
    @(x,y)(x.^2+y.^2+1)./(x.^2+y.^2+1)

Warning 2. In the above list v=[v1; v2; v3; v4] the vertices v1,...,v4 of the quadrilateral

must be arranged in the counterclockwise or clockwise order around

. Namely, for our choice

v1=[0,0]; v2=[5,-1]; v3=[4,5]; v4=[1,4];

we can write

v=[v1; v2; v3; v4] or v=[v4; v3; v2; v1] or v=[v2; v3; v4, v1]

(this will give the same value of the integral),

but we cannot write

v=[v1; v3; v2; v4] or v=[v4; v2; v1; v3],

as it will lead to a wrong value of the determinant evalDetJacb(i) of the Jacobian matrix.

Exercise 1: Build the integQuad function

From the previous example, one can build a function that returns the value of the integral of an inline function defined on a quadrilateral domain (see at the end of this document).

function integralValue=integQuad(F,vertices,N)

where

F: is an inline function previously defined
vertices: are the coordinates of the quadrilateral vertices.
N: order of Gauss integration in each of the axis x, y.

(Hint: find the answer at the end of this document)

Application: Integration over a mesh

Let's consider the integral of

in a general rectangular domain

, in particular we know that

Instead of integrating only in one quadrilateral, we want to show how to approximate the integral value using a mesh subdivision of the domain

The final integral will be the sum of the integrals in each small quadrilateral. Therefore, we increase the precission on the final value.

F = @(x,y) 2*exp(x+y.^2).*y; %present function

xmin = 0; % rectangle dimensions

xmax = 1;

ymin = 0;

ymax = 1;

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

% In order to create a rectangular mesh in Matlab

% we can make the following procedure

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

numDiv = 4; %number of subdivisions

hx = (xmax-xmin)/numDiv;

hy = (ymax-ymin)/numDiv;

xi = xmin:hx:xmax; %all points in the x-axis

eta = ymin:hy:ymax; %all points in the y-axis

[X,Y] = meshgrid(xi,eta); % create a grid of points

Z = F(X,Y); % function values

surf(X,Y,Z); % optional: is just to visualize the result

title('Plot of F(x,y) in the mesh domain')

Hint. Here we plot the surface as the approximation of the continous function using the mesh

[elem,vertex] = surf2patch(X,Y,Z); % the main variables
numElem = size(elem,1); %total number of elements
numVert = size(vertex,1); % total number of vertices
fprintf('num Elements =%d, num vertices = %d \n',numElem,numVert)
num Elements =16, num vertices = 25 

Fig. Visualization of the numbering for vertex and elements for a 4x4 mesh.

Structure of the Mesh

In general a mesh is defined by two tables: vertex and elem.

vertex: is a matrix with 3 columns and numVert rows, where numVert is the number of grid points.

For each row:

the two first coordinates are the (x,y) grid values
the third coordinate is the function value: z=F(x,y)

elem: is a numElem×4 matrix, where numElem is the number of elements For each row:

the four numbers are the number of the vertices defining this element

N = 2; %number of Gauss 2D points = NxN 
integralTot=0;
for i = 1:numElem %compute the integral on each element
    n1 = elem(i,1); %first vertex of the i-th element
    v1 = [vertex(n1,1),vertex(n1,2)];
    n2 = elem(i,2); %second vertex of the i-th element
    v2 = [vertex(n2,1),vertex(n2,2)];
    n3 = elem(i,3); %third vertex of the i-th element
    v3 = [vertex(n3,1),vertex(n3,2)];
    n4 = elem(i,4); %fourth vertex of the i-th element
    v4 = [vertex(n4,1),vertex(n4,2)];
    vertices = [v1;v2;v3;v4];
    elemInteg = integQuad(F,vertices,N);
    integralTot = integralTot+elemInteg;   
end 
actualIntegVal = (exp(1)-1)^2 %exact known value
actualIntegVal = 2.9525
errorInt=abs(actualIntegVal-integralTot) %absolute error
errorInt = 2.9420e-04

Additional Functions

Gauss points in 2D

They are the same as in the 1D version, but now using the points in both x and y axis.

function [w2D,pt2D] = gaussValues2DQuad(n)    
    switch (n)
        case 1
            w = 2; pt = 0;
        case 2
            w = [1,1]; pt = [-1/sqrt(3), 1/sqrt(3)];
        case 3
            w = [5/9, 8/9, 5/9]; pt = [-sqrt(3/5), 0, sqrt(3/5)];
        case 4
            w = [(18+sqrt(30))/36, (18+sqrt(30))/36, (18-sqrt(30))/36, (18-sqrt(30))/36]; 
            pt = [-sqrt(3/7-2/7*(sqrt(6/5))), sqrt(3/7-2/7*(sqrt(6/5))),...
                -sqrt(3/7+2/7*(sqrt(6/5))), sqrt(3/7+2/7*(sqrt(6/5)))];
        case 5
            w = [(322+13*sqrt(70))/900, (322+13*sqrt(70))/900, 128/225,...
                (322-13*sqrt(70))/900,(322-13*sqrt(70))/900]; 
            pt = [-1/3*sqrt(5-2*sqrt(10/7)), 1/3*sqrt(5-2*sqrt(10/7)),0,...
                -1/3*sqrt(5+2*sqrt(10/7)), 1/3*sqrt(5+2*sqrt(10/7))];
        otherwise
            error('No data are defined for this value');
    end
    pt2D = [];
    w2D = [];
    for i = 1:n
        for j = 1:n
            pt2D = [pt2D; [pt(i),pt(j)]];
            w2D = [w2D,w(i)*w(j)];
        end
    end    
end

Quadrilateral integration function: integQuad

Copy this function in a new file named integQuad.m

function valInteg = integQuad(F,vertices,N)
    [w,ptGaussRef] = gaussValues2DQuad(N);
    % 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)]; 
    % evaluate Shape functions on Gaussian reference points
    xi = ptGaussRef(:,1);
    eta = ptGaussRef(:,2);
    evalPsi1 = Psi1(xi,eta);
    evalPsi2 = Psi2(xi,eta);
    evalPsi3 = Psi3(xi,eta);
    evalPsi4 = Psi4(xi,eta);
    % from the change of variables function 
    ptGaussDomain = [evalPsi1,evalPsi2,evalPsi3,evalPsi4]*vertices;
    % evaluate Jacobian contribution for each point 
    for i = 1:size(xi,1)
        evalDetJacb(i) = abs(det(Jacb(xi(i),eta(i))*vertices));
    end
    %evaluate the function on the domain points
    evalF = F(ptGaussDomain(:,1),ptGaussDomain(:,2));
    % Finally, apply Gauss formula
    suma = 0;
    for i = 1:size(ptGaussDomain,1)
        suma = suma+w(i)*evalF(i)*evalDetJacb(i);
    end
    valInteg = suma;
end

Exercise 2: numerical integration error

Increase the number of subdivisions and see how the error evolves.