function eta = effectiveness(u, a, p, t, phi)

% Purpose:    Computes the effectiveness factor,
%
%               eta = Int(kstar*u)/Int(1),
%
%             where Int(kstar*u) denotes the double integral
%             of kstar*u over a cross-section, Omega, of the
%             domain Vstar. Int(1) is thus the area of Omega.
%
% Input data: 'u' is an M x 1 matrix, containing the nodal values of
%             the FEM solution of the reaction-diffusion equation,
%
%               -divstar(Dstar*gradstar(u)) = -(phi^2)*kstar*u,
%
%             in Omega.
%
%             'a' is a text string used to compute the value of the
%             coefficient a = (phi^2)*kstar in Omega. If there is
%             only one subdomain ( = Omega), then 'a' is just a text
%             expression, e.g., a = '1+x+y'. If there are more than
%             one subdomain, 'a' is a sequence of such expressions,
%             separated by exclamation marks. Example: a = '0!1'
%             means that there are two subdomains and that a = 0,
%             in the first, and a = 1, in the second.
%
%             'p' is a 2 x M matrix, the "Point matrix", containing
%             the node coordinates. For 1 <= i <= M,
%             p(1,i) is the x coordinate of node number i,
%             p(2,i) is the y coordinate of node number i.
%
%             't' is a 4 x N matrix, the "Triangle matrix", containing
%             information about the triangles in the triangulation
%             of Omega. For 1 <= n <= N,
%             t(1,n), t(2,n), and t(3,n), are the node numbers of the
%             three corners in triangle number n (counted counter-
%             clockwise), and
%             t(4,n) is the number of the subdomain in which the
%             triangle lies.
%
%             'phi' is the Thiele modulus.
%
%             Note: 'u', 'a', 'p', and 't', should all be exported
%                   from the GUI 'pdetool' to the workspace.
%                   'u' is exported from the "Solve Menu".
%                   'a' is exported from the "PDE Menu".
%                   'p' and 't' are exported from the "Mesh Menu".

nominator = 0;    % initialize the nominator of eta
denominator = 0;  % initialize the denominator of eta

% Find the positions of '!' in 'a':
separators = findstr(a, '!');
  
% "Add" one '!' in the beginning and one '!' in the end:
separators = [0 separators length(a)+1];

N = size(t, 2);  % number of triangles

for n = 1:N % loop over all triangles

  % Compute the node numbers and coordinates
  % of the three corners in triangle number n:
  for index = 1:3
    node(index) = t(index, n);
    x(index) = p(1, node(index));
    y(index) = p(2, node(index));
  end

  % Compute the area of the triangle:
  area = ((x(2) - x(1))*(y(3) - y(1)) - ...
          (x(3) - x(1))*(y(2) - y(1)))/2;
  
  % Compute the coordinates of the centre of gravity:
  xmid = (x(1) + x(2) + x(3))/3;
  ymid = (y(1) + y(2) + y(3))/3;

  % Compute the value of 'u' in the centre of gravity:
  umid = (u(node(1)) + u(node(2)) + u(node(3)))/3; % since 'u' is linear
                                                   % on each triangle

  % Compute the number of the subdomain in which the triangle lies:
  subdomain = t(4, n);

  % Compute starting position and ending position in 'a'
  % for this subdomain:
  startpos = separators(subdomain) + 1;
  endpos = separators(subdomain+1) - 1;
  
  % Compute the value of 'kstar' in the centre of gravity:
  x = xmid;
  y = ymid;
  kstarmid = eval(a(startpos:endpos))/(phi^2);

  % Add contributions to the integrals from triangle n:
  nominator = nominator + kstarmid*umid*area;
  denominator = denominator + area;
  
end

eta = nominator/denominator;
