function I=MailleInt(p,t,d)
%MESHINTEGRATE Integrate values on a mesh
%   I=MESHINTEGRATE(P,T,D) computes the integral I over the mesh
%   given by P and T, with values (for each point) in D. D is of size
%   1-by-np, where np is the number of points in P (=size(P,2)).
%   The elements are considered to be linear.
%
%   I=MESHINTEGRATE(P,D) assumes T=[1,2,3,... (NP-1) ; 2,3,4,... NP],
%   (where np=size(p,2)), i.e., that the mesh is a line and that the
%   points in P are sorted.
%
%   I=MESHINTEGRATE(P) calls MESHINTEGRATE(P(1,:),P(2,:)).
%
%   This function is useful for computing integrals along cross-sections
%   plotted with postcrossplot, in which case P, T and D are extracted from
%   the output when 'outtype' is 'postdata'.
%
%   Examples:
%
%   Line integral in 2D:
%   --------------------
%   % Just set up a problem:
%   clear fem
%   fem.geom = circ2+rect2;
%   fem.mesh = meshinit(fem);
%   fem.shape = 2; fem.equ.c = 1; fem.equ.f = 1; fem.bnd.h = 1;
%   fem.xmesh = meshextend(fem);
%   fem.sol = femlin(fem);
%
%   % Make a cross-section plot, with output being a postdata structure
%   pd = postcrossplot(fem,1,[0 1;0 1],'lindata','u','npoints',100,'outtype','postdata');
%
%   % Call meshintegrate:
%   I = meshintegrate(pd.p);
%
%   Line integral in 3D:
%   --------------------
%   % Just set up a problem:
%   clear fem, fem.geom = block3;
%   fem.mesh = meshinit(fem,'hmax',0.15);
%   fem.shape = 2;
%   fem.equ.c = 1; fem.equ.f = 1;
%   fem.bnd.h = {1 1 0 0 1 1};
%   fem.xmesh = meshextend(fem);
%   fem.sol = femlin(fem);
%
%   % Make cross-section plot:
%   pd = postcrossplot(fem,1,[0 1;0 1;0 1],'lindata','u','npoints',100,'outtype','postdata');
%
%   % Call meshintegrate:
%   I = meshintegrate(pd.p)
%
%   Surface integral in 3D using the same problem as above:
%   -------------------------------------------------------
%   pd = postcrossplot(fem,2,[0 0 0;0 1 0;1 0 1]','surfdata','u','outtype','postdata');
%   I = meshintegrate(pd.p, pd.t, pd.d)
%
%   This method only works for lines and surfaces actually
%   intersecting the geometry. For plots along geometry boundaries
%   or edges (or 1-D subdomains), better results are achieved using
%   POSTINT.
%
%   See also POSTCROSSPLOT, POSTINT

%   Copyright (c) 1998-2008 by COMSOL AB
%   $Revision: 1.9.2.1 $  $Date: 2008/10/30 16:25:03 $

np=size(p,2);
if nargin==1
  d=p(2,:);
  p=p(1,:);
  t=[1:(np-1) ; 2:np];
elseif nargin==2
  d=t;
  t=[1:(np-1) ; 2:np];
end

for i=1:(size(t,2)./(3^size(p,1)));
    tri=(i-1).*(3.^size(p,1))+([1:(3^size(p,1))]);
    vol=l_area(p,t(:,tri));
    I(i)=sum(mean(reshape(d(t(:,tri)),size(t(:,tri))),1).*vol);
end

    
function A=l_area(p,t)
edim=size(t,1)-1;
switch edim
case 3
  i0 = t(1,:);
  i1 = t(2,:);
  i2 = t(3,:);
  i3 = t(4,:);
  x0 = p(1,i0);
  y0 = p(2,i0);
  z0 = p(3,i0);
  x1 = p(1,i1);
  y1 = p(2,i1);
  z1 = p(3,i1);
  x2 = p(1,i2);
  y2 = p(2,i2);
  z2 = p(3,i2);
  x3 = p(1,i3);
  y3 = p(2,i3);
  z3 = p(3,i3);
  dx0 = x1-x0;
  dy0 = y1-y0;
  dz0 = z1-z0;
  dx1 = x2-x1;
  dy1 = y2-y1;
  dz1 = z2-z1;
  dx5 = x3-x1;
  dy5 = y3-y1;
  dz5 = z3-z1;
  A = abs((dx0.*(dy1.*dz5-dz1.*dy5)+dy0.*(dz1.*dx5-dx1.*dz5)+dz0.*(dx1.*dy5-dy1.*dx5))/6);
case 2
  a=sqrt(sum((p(:,t(1,:))-p(:,t(2,:))).^2,1));
  b=sqrt(sum((p(:,t(2,:))-p(:,t(3,:))).^2,1));
  c=sqrt(sum((p(:,t(3,:))-p(:,t(1,:))).^2,1));
  s=(a+b+c)/2;
  A=sqrt(s.*(s-a).*(s-b).*(s-c));
case 1
  A=sqrt(sum((p(:,t(1,:))-p(:,t(2,:))).^2,1));
otherwise
  error('Not implemented.')
end