general/geometry/cubequadrature.m

Go to the documentation of this file.

function res = cubequadrature(poldeg, func, varargin)
%function res = cubequadrature(poldeg, func, varargin)
% integration of function func over reference cube == unit cube.
% by Gaussian quadrature exactly integrating polynoms of poldeg.
%
% arguments:
%  poldeg:  degree of polynomial the quadrature rule exactly approximates
%           (0-23)
%  dim:     dimension of cube on which we want to integrate (>=2)
%  func:    is a function getting a local coordinate vector (1d) and giving a
%           (vectorial or scalar) result
%  varargin: optional further arguments for function
%
% return values:
%  res:     result of quadrature
%

[points, weights]= get_quadrature_weights_1d(poldeg);

n = length(points);

weightss = weights';
pointss  = points;

dim = 2;

for i = 1:dim-1
  m = size(pointss,1);
  weightss = kron(weightss, weights);
  weightss = reshape(weightss, m*n, 1);
  pointss  = [repmat(pointss, n, 1), reshape(repmat(points,1,m)',n*m,1)];
end

res = quadrature(weightss,pointss,func,varargin{:});