BIM2Aadvdiff

% -*- texinfo -*-
%
% @deftypefn {Function File} @
% {@var{A}} = BIM2Aadvdiff (@var{mesh}, @var{alpha}, @var{gamma}, @var{eta}, @var{beta})
%
% Builds the Scharfetter-Gummel matrix for the 
% discretization of the LHS
% of the equation
% @iftex 
% @tex
% $ -div ( \alpha  \gamma  ( \eta \vect{\nabla} u - \vect{beta} u )) = f $
% @end tex 
% @end iftex 
% @ifnottex
% -div (@var{alpha} * @var{gamma} (@var{eta} grad u - @var{beta} u )) = f
% @end ifnottex
% 
% where: 
% @itemize @minus
% @item @var{alpha}: element-wise constant scalar function
% @item @var{eta}, @var{gamma}: piecewise linear conforming 
% scalar functions
% @item @var{beta}: element-wise constant vector function
% @end itemize
%
% Instead of passing the vector field @var{beta} directly
% one can pass a piecewise linear conforming scalar function
% @var{phi} as the last input.  In such case @var{beta} = grad @var{phi}
% is assumed.  If @var{phi} is a single scalar value @var{beta}
% is assumed to be 0 in the whole domain.
% 
% Example:
% @example
% mesh = MSH2Mstructmesh([0:1/3:1],[0:1/3:1],1,1:4);
% mesh = BIM2Cmeshproperties(mesh);
% x = mesh.p(1,:)';
% Dnodes = BIM2Cunknownsonside(mesh,[2,4]);
% Nnodes = columns(mesh.p); Nelements = columns(mesh.t);
% Varnodes = setdiff(1:Nnodes,Dnodes);
% alpha  = ones(Nelements,1); eta = .1*ones(Nnodes,1);
% beta   = [ones(1,Nelements);zeros(1,Nelements)];
% gamma  = ones(Nnodes,1);
% f      = BIM2Arhs(mesh,ones(Nnodes,1),ones(Nelements,1));
% S = BIM2Aadvdiff(mesh,alpha,gamma,eta,beta);
% u = zeros(Nnodes,1);
% u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes);
% uex = x - (exp(10*x)-1)/(exp(10)-1);
% assert(u,uex,1e-7)
% @end example
%
% @seealso{BIM2Arhs, BIM2Areaction, BIM2Cmeshproperties}
% @end deftypefn