Numerical Integration - Numerical Methods in Engineering with MATLAB - page 239

orderedina counterclockwise direction around the element. The determinantof
the Jacobian matrix isobtainedbycalling detJ ; mapping is performedby map .
The weights and the values of at the integrationpoints arecomputedby

gaussNodes listedinthe previous article (note that

ξ and η

ξ and η appear as s and t in

listing).

functionI=gaussQuad2(func,x,y,n)
% Gauss-Legendre quadrature over a quadrilateral.
%USAGE:I=gaussQuad2(func,x,y,n)
% INPUT:
%func=handleoffunctiontobeintegrated.
% x = [x1;x2;x3;x4] = x-coordinates of corners.
% y = [y1;y2;y3;y4] = y-coordinates of corners.
% n = order of integration
% OUTPUT:
% I

= integral

[t,A]=gaussNodes(n);I=0;
fori=1:n
forj=1:n
[xNode,yNode] = map(x,y,t(i),t(j));
z = feval(func,xNode,yNode);
detJ = jac(x,y,t(i),t(j));
I=I+A(i)*A(j)*detJ*z;

end

end

function detJ = jac(x,y,s,t)
% Computes determinant of Jacobian matrix.
J = zeros(2);
J(1,1)=-(1-t)*x(1)+(1-t)*x(2)...
+(1+t)*x(3)-(1+t)*x(4);
J(1,2)=-(1-t)*y(1)+(1-t)*y(2)...
+(1+t)*y(3)-(1+t)*y(4);
J(2,1)=-(1-s)*x(1)-(1+s)*x(2)...
+(1+s)*x(3)+(1-s)*x(4);
J(2,2)=-(1-s)*y(1)-(1+s)*y(2)...
+(1+s)*y(3)+(1-s)*y(4);
detJ = (J(1,1)*J(2,2) - J(1,2)*J(2,1))/16;