octave_packages/secs2d-0.0.8/Utilities/Ufvsgcurrent3.m

   1 function [Fx,Fy]=Ufvsgcurrent3(mesh,u,alpha,gamma,eta,beta);
   2
   3
   4   ## -*- texinfo -*-
   5   ##
   6   ## @deftypefn {Function File} {[@var{jx},@var{jy}]} = Ufvsgcurrent3 @
   7   ## (@var{mesh}, @var{u}, @var{alpha}, @var{gamma}, @var{eta}, @var{beta});
   8   ##
   9   ##
  10   ## Builds the Scharfetter-Gummel approximation of the vector
  11   ## field
  12   ##
  13   ## @iftex
  14   ## @tex
  15   ## $ \vect{J}(u) = \alpha \gamma (\eta\vect{\nabla}u-\vect{beta}u) $
  16   ## @end tex
  17   ## @end iftex
  18   ## @ifinfo
  19   ## J(@var{u}) = @var{alpha}* @var{gamma} * (@var{eta} * grad @var{u} - @var{beta} * @var{u}))
  20   ## @end ifinfo
  21   ##
  22   ## where:
  23   ## @itemize @minus
  24   ## @item @var{alpha} is an element-wise constant scalar function
  25   ## @item @var{eta}, @var{u}, @var{gamma} are piecewise linear
  26   ## conforming scalar functions
  27   ## @item @var{beta}  is an element-wise constant vector function
  28   ## @end itemize
  29   ##
  30   ## J(@var{u}) is an element-wise constant vector function
  31   ##
  32   ## Instead of passing the vector field @var{beta} directly
  33   ## one can pass a piecewise linear conforming scalar function
  34   ## @var{phi} as the last input.  In such case @var{beta} = grad @var{phi}
  35   ## is assumed.  If @var{phi} is a single scalar value @var{beta}
  36   ## is assumed to be 0 in the whole domain.
  37   ##
  38   ## @seealso{Uscharfettergummel3}
  39   ## @end deftypefn
  40
  41   %% This file is part of
  42   %%
  43   %%            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
  44   %%         -------------------------------------------------------------------
  45   %%            Copyright (C) 2004-2006  Carlo de Falco
  46   %%
  47   %%
  48   %%
  49   %%  SECS2D is free software; you can redistribute it and/or modify
  50   %%  it under the terms of the GNU General Public License as published by
  51   %%  the Free Software Foundation; either version 2 of the License, or
  52   %%  (at your option) any later version.
  53   %%
  54   %%  SECS2D is distributed in the hope that it will be useful,
  55   %%  but WITHOUT ANY WARRANTY; without even the implied warranty of
  56   %%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  57   %%  GNU General Public License for more details.
  58   %%
  59   %%  Youx should have received a copy of the GNU General Public License
  60   %%  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
  61
  62   Nelem  = columns(mesh.t);
  63   Nnodes = columns(mesh.p);
  64
  65   uloc      = u(mesh.t(1:3,:));
  66
  67   shgx = reshape(mesh.shg(1,:,:),3,Nelem);
  68   shgy = reshape(mesh.shg(2,:,:),3,Nelem);
  69
  70   x      = reshape(mesh.p(1,mesh.t(1:3,:)),3,[]);
  71   dx     = [ (x(3,:)-x(2,:)) ;
  72             (x(1,:)-x(3,:)) ;
  73             (x(2,:)-x(1,:)) ];
  74
  75   y      = reshape(mesh.p(2,mesh.t(1:3,:)),3,[]);
  76   dy     = [ (y(3,:)-y(2,:)) ;
  77             (y(1,:) -y(3,:)) ;
  78             (y(2,:) -y(1,:)) ];
  79
  80   ##ei     = sqrt( dx.^2 + dy.^2 );
  81
  82   if all(size(beta)==1)
  83     v12=0;v23=0;v31=0;
  84   elseif all(size(beta)==[2,Nelem])
  85     v23    = beta(1,:) .* dx(1,:) + beta(2,:) .* dy(1,:);
  86     v31    = beta(1,:) .* dx(2,:) + beta(2,:) .* dy(2,:);
  87     v12    = beta(1,:) .* dx(3,:) + beta(2,:) .* dy(3,:);
  88   elseif all(size(beta)==[Nnodes,1])
  89     betaloc = beta(mesh.t(1:3,:));
  90     v23    = betaloc(3,:)-betaloc(2,:);
  91     v31    = betaloc(1,:)-betaloc(3,:);
  92     v12    = betaloc(2,:)-betaloc(1,:);
  93   end
  94
  95   etaloc = eta(mesh.t(1:3,:));
  96
  97   eta23    = etaloc(3,:)-etaloc(2,:);
  98   eta31    = etaloc(1,:)-etaloc(3,:);
  99   eta12    = etaloc(2,:)-etaloc(1,:);
 100
 101   etalocm1 = Utemplogm(etaloc(2,:),etaloc(3,:));
 102   etalocm2 = Utemplogm(etaloc(3,:),etaloc(1,:));
 103   etalocm3 = Utemplogm(etaloc(1,:),etaloc(2,:));
 104
 105   gammaloc = gamma(mesh.t(1:3,:));
 106   geloc      = gammaloc.*etaloc;
 107
 108   gelocm1 = Utemplogm(geloc(2,:),geloc(3,:));
 109   gelocm2 = Utemplogm(geloc(3,:),geloc(1,:));
 110   gelocm3 = Utemplogm(geloc(1,:),geloc(2,:));
 111
 112   [bp23,bm23] = Ubern( (v23 - eta23)./etalocm1);
 113   [bp31,bm31] = Ubern( (v31 - eta31)./etalocm2);
 114   [bp12,bm12] = Ubern( (v12 - eta12)./etalocm3);
 115
 116   gfigfj = [ shgx(3,:) .* shgx(2,:) + shgy(3,:) .* shgy(2,:) ;
 117             shgx(1,:) .* shgx(3,:) + shgy(1,:) .* shgy(3,:) ;
 118             shgx(2,:) .* shgx(1,:) + shgy(2,:) .* shgy(1,:) ];
 119
 120   Fx = - alpha' .* ( gelocm1 .* etalocm1 .* dx(1,:) .*  ...
 121                     gfigfj(1,:) .* ...
 122                    ( bp23 .* uloc(3,:)./etaloc(3,:) -...
 123                     bm23 .* uloc(2,:)./etaloc(2,:)) +... %% 1
 124                     gelocm2 .* etalocm2 .* dx(2,:) .*  ...
 125                     gfigfj(2,:) .* ...
 126                    (bp31 .* uloc(1,:)./etaloc(1,:) -...
 127                     bm31 .* uloc(3,:)./etaloc(3,:)) +... %% 2
 128                     gelocm3 .* etalocm3 .* dx(3,:) .* ...
 129                     gfigfj(3,:) .* ...
 130                    (bp12 .* uloc(2,:)./etaloc(2,:) -...
 131                     bm12 .* uloc(1,:)./etaloc(1,:)) ... %% 3
 132                    );
 133
 134   Fy = - alpha' .* ( gelocm1 .* etalocm1 .* dy(1,:) .*  ...
 135                     gfigfj(1,:) .* ...
 136                    ( bp23 .* uloc(3,:)./etaloc(3,:) -...
 137                     bm23 .* uloc(2,:)./etaloc(2,:)) +... %% 1
 138                     gelocm2 .* etalocm2 .* dy(2,:) .*  ...
 139                     gfigfj(2,:) .* ...
 140                     (bp31 .* uloc(1,:)./etaloc(1,:) -...
 141                      bm31 .* uloc(3,:)./etaloc(3,:)) +... %% 2
 142                     gelocm3 .* etalocm3 .* dy(3,:) .* ...
 143                     gfigfj(3,:) .* ...
 144                     (bp12 .* uloc(2,:)./etaloc(2,:) -...
 145                      bm12 .* uloc(1,:)./etaloc(1,:)) ... %% 3
 146                     );