www.pudn.com > 87361026FEA.rar > quadrature.m, change:2002-10-01,size:9368b
function [W,Q] = quadrature( quadorder, qt, sdim )
% The function quadrature returns a n x 1 column vector W of quadrature
% weights and a n x dim matrix of quadrature points, where n is the
% number of quadrature points. The function is called as follows:
%
% [W,Q]=quadrature( nint, type, dim )
%
% nint is the quadrature order, type is the type of quadrature
% (i.e. gaussian, triangular, etc.. ) and dim is the number of spacial
% dimentions of the problem. The default for type is GAUSS and the
% default for dim is unity.
%
% wrQ=quadrature(nint,'TRIANGULAR',2);itten by Jack Chessa
% j-chessa@northwestern.edu
% Department of Mechanical Engineering
% Northwestern University
if ( nargin < 3 ) % set default arguments
if ( strcmp(qt,'GAUSS') == 1 )
dim = 1;
else
dim = 2;
end
end
if ( nargin < 2 )
type = 'GAUSS';
end
if ( strcmp(qt,'GAUSS') == 1 )
if ( quadorder > 8 ) % check for valid quadrature order
disp('Order of quadrature to high for Gaussian Quadrature');
quadorder =8;
end
quadpoint=zeros(quadorder^sdim ,sdim);
quadweight=zeros(quadorder^sdim,1);
r1pt=zeros(quadorder,1); r1wt=zeros(quadorder,1);
switch ( quadorder )
case 1
r1pt(1) = 0.000000000000000;
r1wt(1) = 2.000000000000000;
case 2
r1pt(1) = 0.577350269189626;
r1pt(2) =-0.577350269189626;
r1wt(1) = 1.000000000000000;
r1wt(2) = 1.000000000000000;
case 3
r1pt(1) = 0.774596669241483;
r1pt(2) =-0.774596669241483;
r1pt(3) = 0.000000000000000;
r1wt(1) = 0.555555555555556;
r1wt(2) = 0.555555555555556;
r1wt(3) = 0.888888888888889;
case 4
r1pt(1) = 0.861134311594053;
r1pt(2) =-0.861134311594053;
r1pt(3) = 0.339981043584856;
r1pt(4) =-0.339981043584856;
r1wt(1) = 0.347854845137454;
r1wt(2) = 0.347854845137454;
r1wt(3) = 0.652145154862546;
r1wt(4) = 0.652145154862546;
case 5
r1pt(1) = 0.906179845938664;
r1pt(2) =-0.906179845938664;
r1pt(3) = 0.538469310105683;
r1pt(4) =-0.538469310105683;
r1pt(5) = 0.000000000000000;
r1wt(1) = 0.236926885056189;
r1wt(2) = 0.236926885056189;
r1wt(3) = 0.478628670499366;
r1wt(4) = 0.478628670499366;
r1wt(5) = 0.568888888888889;
case 6
r1pt(1) = 0.932469514203152;
r1pt(2) =-0.932469514203152;
r1pt(3) = 0.661209386466265;
r1pt(4) =-0.661209386466265;
r1pt(5) = 0.238619186003152;
r1pt(6) =-0.238619186003152;
r1wt(1) = 0.171324492379170;
r1wt(2) = 0.171324492379170;
r1wt(3) = 0.360761573048139;
r1wt(4) = 0.360761573048139;
r1wt(5) = 0.467913934572691;
r1wt(6) = 0.467913934572691;
case 7
r1pt(1) = 0.949107912342759;
r1pt(2) = -0.949107912342759;
r1pt(3) = 0.741531185599394;
r1pt(4) = -0.741531185599394;
r1pt(5) = 0.405845151377397;
r1pt(6) = -0.405845151377397;
r1pt(7) = 0.000000000000000;
r1wt(1) = 0.129484966168870;
r1wt(2) = 0.129484966168870;
r1wt(3) = 0.279705391489277;
r1wt(4) = 0.279705391489277;
r1wt(5) = 0.381830050505119;
r1wt(6) = 0.381830050505119;
r1wt(7) = 0.417959183673469;
case 8
r1pt(1) = 0.960289856497536;
r1pt(2) = -0.960289856497536;
r1pt(3) = 0.796666477413627;
r1pt(4) = -0.796666477413627;
r1pt(5) = 0.525532409916329;
r1pt(6) = -0.525532409916329;
r1pt(7) = 0.183434642495650;
r1pt(8) = -0.183434642495650;
r1wt(1) = 0.101228536290376;
r1wt(2) = 0.101228536290376;
r1wt(3) = 0.222381034453374;
r1wt(4) = 0.222381034453374;
r1wt(5) = 0.313706645877887;
r1wt(6) = 0.313706645877887;
r1wt(7) = 0.362683783378362;
r1wt(8) = 0.362683783378362;
otherwise
disp('Order of quadrature to high for Gaussian Quadrature');
end % end of quadorder switch
n=1;
if ( sdim == 1 )
for i = 1:quadorder
quadpoint(n,:) = [ r1pt(i) ];
quadweight(n) = r1wt(i);
n = n+1;
end
elseif ( sdim == 2 )
for i = 1:quadorder
for j = 1:quadorder
quadpoint(n,:) = [ r1pt(i), r1pt(j)];
quadweight(n) = r1wt(i)*r1wt(j);
n = n+1;
end
end
else % sdim == 3
for i = 1:quadorder
for j = 1:quadorder
for k = 1:quadorder
quadpoint(n,:) = [ r1pt(i), r1pt(j), r1pt(k) ];
quadweight(n) = r1wt(i)*r1wt(j)*r1wt(k);
n = n+1;
end
end
end
end
Q=quadpoint;
W=quadweight;
% END OF GAUSSIAN QUADRATURE DEFINITION
elseif ( strcmp(qt,'TRIANGULAR') == 1 )
if ( sdim == 3 ) %%% TETRAHEDRA
if ( quadorder ~= 1 & quadorder ~= 2 & quadorder ~= 3 )
% check for valid quadrature order
disp('Incorect quadrature order for triangular quadrature');
quadorder = 1;
end
if ( quadorder == 1 )
quadpoint = [ 0.25 0.25 0.25 ];
quadweight = 1;
elseif ( quadorder == 2 )
quadpoint = [ 0.58541020 0.13819660 0.13819660;
0.13819660 0.58541020 0.13819660;
0.13819660 0.13819660 0.58541020;
0.13819660 0.13819660 0.13819660];
quadweight = [1; 1; 1; 1]/4;
elseif ( quadorder == 3 )
quadpoint = [ 0.25 0.25 0.25;
1/2 1/6 1/6;
1/6 1/2 1/6;
1/6 1/6 1/2;
1/6 1/6 1/6];
quadweight = [-4/5 9/20 9/20 9/20 9/20]';
end
Q=quadpoint;
W=quadweight/6;
else %%% TRIANGLES
if ( quadorder > 7 ) % check for valid quadrature order
disp('Quadrature order too high for triangular quadrature');
quadorder = 1;
end
if ( quadorder == 1 ) % set quad points and quadweights
quadpoint = [ 0.3333333333333, 0.3333333333333 ];
quadweight = 1;
elseif ( quadorder == 2 )
quadpoint = zeros( 3, 2 );
quadweight = zeros( 3, 1 );
quadpoint(1,:) = [ 0.1666666666667, 0.1666666666667 ];
quadpoint(2,:) = [ 0.6666666666667, 0.1666666666667 ];
quadpoint(3,:) = [ 0.1666666666667, 0.6666666666667 ];
quadweight(1) = 0.3333333333333;
quadweight(2) = 0.3333333333333;
quadweight(3) = 0.3333333333333;
elseif ( quadorder <= 5 )
quadpoint = zeros( 7, 2 );
quadweight = zeros( 7, 1 );
quadpoint(1,:) = [ 0.1012865073235, 0.1012865073235 ];
quadpoint(2,:) = [ 0.7974269853531, 0.1012865073235 ];
quadpoint(3,:) = [ 0.1012865073235, 0.7974269853531 ];
quadpoint(4,:) = [ 0.4701420641051, 0.0597158717898 ];
quadpoint(5,:) = [ 0.4701420641051, 0.4701420641051 ];
quadpoint(6,:) = [ 0.0597158717898, 0.4701420641051 ];
quadpoint(7,:) = [ 0.3333333333333, 0.3333333333333 ];
quadweight(1) = 0.1259391805448;
quadweight(2) = 0.1259391805448;
quadweight(3) = 0.1259391805448;
quadweight(4) = 0.1323941527885;
quadweight(5) = 0.1323941527885;
quadweight(6) = 0.1323941527885;
quadweight(7) = 0.2250000000000;
else
quadpoint = zeros( 13, 2 );
quadweight = zeros( 13, 1 );
quadpoint(1 ,:) = [ 0.0651301029022, 0.0651301029022 ];
quadpoint(2 ,:) = [ 0.8697397941956, 0.0651301029022 ];
quadpoint(3 ,:) = [ 0.0651301029022, 0.8697397941956 ];
quadpoint(4 ,:) = [ 0.3128654960049, 0.0486903154253 ];
quadpoint(5 ,:) = [ 0.6384441885698, 0.3128654960049 ];
quadpoint(6 ,:) = [ 0.0486903154253, 0.6384441885698 ];
quadpoint(7 ,:) = [ 0.6384441885698, 0.0486903154253 ];
quadpoint(8 ,:) = [ 0.3128654960049, 0.6384441885698 ];
quadpoint(9 ,:) = [ 0.0486903154253, 0.3128654960049 ];
quadpoint(10,:) = [ 0.2603459660790, 0.2603459660790 ];
quadpoint(11,:) = [ 0.4793080678419, 0.2603459660790 ];
quadpoint(12,:) = [ 0.2603459660790, 0.4793080678419 ];
quadpoint(13,:) = [ 0.3333333333333, 0.3333333333333 ];
quadweight(1 ) = 0.0533472356088;
quadweight(2 ) = 0.0533472356088;
quadweight(3 ) = 0.0533472356088;
quadweight(4 ) = 0.0771137608903;
quadweight(5 ) = 0.0771137608903;
quadweight(6 ) = 0.0771137608903;
quadweight(7 ) = 0.0771137608903;
quadweight(8 ) = 0.0771137608903;
quadweight(9 ) = 0.0771137608903;
quadweight(10) = 0.1756152576332;
quadweight(11) = 0.1756152576332;
quadweight(12) = 0.1756152576332;
quadweight(13) =-0.1495700444677;
end
Q=quadpoint;
W=quadweight/2;
end
end % end of TRIANGULAR initialization
% END OF FUNCTION