general/geometry/quadrature.m

function res = quadrature(weights,points,func,varargin)
%function res = quadrature(weights,points,func,varargin)
% 
% integration of function func by given quadrature. func 
% is a function getting a local coordinate vector and varargin
% and giving a (vectorial or scalar) result or a cell array of such.
% points is expected to be a npoints x dimpoint matrix
% result is the same size of f-evaluations.
% we emphasize, that also cell-array valued functions can be
% integrated (for use in rb-methods)

% Bernard Haasdonk 27.8.2009

npoints = length(weights);
if isempty(varargin)
  f = func(points(1,:));
  if ~iscell(f)
    res = weights(1)*f;
    for qp=2:npoints
      f = func(points(qp,:));
      res = res +weights(qp)*f;
    end;  
  else % iscell!!
    res = f;
    for q = 1:length(f(:));
      res{q} = weights(1)*f{q};
    end;
    for qp=2:npoints
      f = func(points(qp,:));
      for q = 1:length(f(:));
        res{q} = res{q} +weights(qp)*f{q};
      end;
    end;  
  end;
else % not isempty varargin: same but different...
  f = func(points(1,:),varargin{:});
  if ~iscell(f)
    res = weights(1)*f;
    for qp=2:npoints
      f = func(points(qp,:),varargin{:});
      res = res +weights(qp)*func(points(qp,:),varargin{:});
    end;
  else % iscell
    for q = 1:length(f(:));
      res{q} = weights(1)*f{q};
    end;
    for qp=2:npoints
      f = func(points(qp,:),varargin{:});
      for q = 1:length(f(:));
        res{q} = res{q} +weights(qp)*f{q};
      end;
    end;
  end;
end;

%| \docupdate