Uscharfettergummel3

% -*- texinfo -*-
%
% @deftypefn {Function File} @
% {@var{S}} = Uscharfettergummel3 (@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 
% @ifinfo
% -div (@var{alpha} * @var{gamma} (@var{eta} grad u - @var{beta} u )) = f
% @end ifinfo
% 
% where: 
% @itemize @minus
% @item @var{alpha} is an element-wise constant scalar function
% @item @var{eta}, @var{gamma} are piecewise linear conforming 
% scalar functions
% @item @var{beta}  is an 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.p,mesh.e,mesh.t] = Ustructmesh([0:1/3:1],[0:1/3:1],1,1:4);
% mesh = Umeshproperties(mesh);
% x = mesh.p(1,:)';
% Dnodes = Unodesonside(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      = Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1));
% S = Uscharfettergummel3(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{Ucomplap, Ucompconst, Ucompmass2, Uscharfettergummel}
% @end deftypefn