2D Gaussian quadratures: Quadrilateral Example

Contents

Example of 2D integration

Let's consider the function

$$ F(x,y)= x^2 - y^2$ defined on the quadrilateral $$\Omega_k$ with vertices

$$v_1=(0,0), \quad v_2=(5,-1), \quad v_3=(4,5), \quad v_4=(1,4)$.

We want to compute

$$\int_{\Omega_k} F(x,y) dx dy$.

Below, you have an slide showing the change of variables needed to relate the reference quadrilateral [-1,1]x[-1,1] with a general one.

Compute the 2D Gauss points on the reference element

First we compute the appropriate Gauss points in the reference quadrilateral.

We can use a Gauss quadrature using only N=2 in this example, because $$F(x,y)$ is a polynomial function of degree less than 3 in each variable.

N=2; %order of the Gaussian quadrature
[w,ptGaussRef]=gaussValues2DQuad(N);
% Draw Gauss points in the reference quadrilateral
plotRectangle([-1 -1],[1,-1],[1,1],[-1,1]);
hold on;
plot(ptGaussRef(:,1),ptGaussRef(:,2),'ro');
hold off;
%

A first example on the reference quadrilateral

Consider our function $$F(x)=x^2y^2$ and we code in Matlab the above change of variable formulas:

N=2;
[w,ptGaussRef]=gaussValues2DQuad(N)
F=@(x,y) (x.^2).*y.^2;
Fxy=F(ptGaussRef(:,1),ptGaussRef(:,2));
Integ=sum(w'.*Fxy)
w =

     1     1     1     1


ptGaussRef =

   -0.5774   -0.5774
   -0.5774    0.5774
    0.5774   -0.5774
    0.5774    0.5774


Integ =

    0.4444

General case: Define the function $$F(x,y)$ and the domain $$\Omega_k$

Now consider our function $$F(x)=x^2-y^2$ 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];

Compute the change of variables functions

As it is shown in the previous slide, the needed functions to define the change of variables$\psi_i(x,y)$ are known as shape functions. These functions are associated to the respective vertex of the quadrilateral and they satisfy

$$\psi_i(v_j)= 0 \quad$ if $$\quad i \neq j$ and

$$\psi_i(v_j)= 1 \quad$ if $$\quad i = j$

% 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

The change of variables from the reference quadrilateral to a genral one is:

$$(x,y)=\psi_1(\xi, \eta) v_1+\psi_2(\xi, \eta) v_2+\psi_3(\xi, \eta)
v_3+\psi_4(\xi, \eta) v_4$$

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*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 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

Exercise 1: Build the integQuad function

From the previous example, 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

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

Application: Integration over a mesh

Let's consider the integral of $$F(x,y)= 2 e^{x+y^2} y$ in a rectangular domain [xmin,xmax]x[ymin,ymax]

We want to show how to approximate the integral value using a mesh subdivision of the domain.

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
[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
numElem =

    16


numVert =

    25

Mesh definition

vertex: is a numPx3 matrix, where numP is the number of grid points. For each row:

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

N=2; %number of Gauss points = NxN
integralTot=0;
for i=1:numElem %compute the integral on each element
    v1=[vertex(elem(i,1),1),vertex(elem(i,1),2)];
    v2=[vertex(elem(i,2),1),vertex(elem(i,2),2)];
    v3=[vertex(elem(i,3),1),vertex(elem(i,3),2)];
    v4=[vertex(elem(i,4),1),vertex(elem(i,4),2)];
    vertices=[v1;v2;v3;v4];
    elemInteg=integQuad(F,vertices,N);
    integralTot=integralTot+elemInteg;
end
actualIntegVal= (exp(1)-1)^2 %exact known value
errorInt=abs(actualIntegVal-integralTot) %absolute error

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
actualIntegVal =

    2.9525


errorInt =

   2.9420e-04

Exercise 2:

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

(c) Numerical Factory 2018