function [D_face] = DiffusionCoefficient(D,I,J,DiffusionDiscrModel,varargin)

% Evaluates the diffusion coeff. (function) at the edges of the control
% volumes using the input variables. DiffusionDiscrModel could be either
% one of the 'HarmonicAverage' , 'ArithmeticAverage' methods.
%
% Two important notes: 
%   1- The independent variables should be defined inside the domain only.
%   2- The function should be defined using Element-By-Element notations.
%
% For example, if diffusion coeff. 'D(a,x) = a.*x' depends on species 'A'
% and x-coordinate, you should first restrict them to the nodes inside
% of the domain: 
%      a = A(InsideOldMask);   x = PosX(InsideOldMask)
% and then send them to this function:
%     [D_face]= DiffusionCoefficient(D,I,J,DiffCoefApproxMeth,a,x)

% evaluate diffusion coeff. function at the cell centers
D_cell = D(varargin{:});

% evaluate the diffusion coeff. at the edges of the control volumes.
if strcmp(DiffusionDiscrModel,'HarmonicAverage')
  D_face = 2*( D_cell(I,:).*D_cell(J,:) ./ ( D_cell(I,:)+D_cell(J,:) + 1e-16 ) );
  
else                        % 'ArithmeticAverage'
  D_face = 1/2 *( D_cell(I,:)+D_cell(J,:) );
  
end


end % end function 'DiffusionCoefficient'