octave_packages/secs2d-0.0.8/Utilities/Ufvsgcurrent.m

   1 function [jx,jy]=Ufvsgcurrent(omesh,n,psi,psith,coeffe);
   2
   3 % [jx,jy]=Udrawcurrent(omesh,n,psi,psith,coeffe);
   4
   5 % This file is part of
   6 %
   7 %            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
   8 %         -------------------------------------------------------------------
   9 %            Copyright (C) 2004-2006  Carlo de Falco
  10 %
  11 %
  12 %
  13 %  SECS2D is free software; you can redistribute it and/or modify
  14 %  it under the terms of the GNU General Public License as published by
  15 %  the Free Software Foundation; either version 2 of the License, or
  16 %  (at your option) any later version.
  17 %
  18 %  SECS2D is distributed in the hope that it will be useful,
  19 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 %  GNU General Public License for more details.
  22 %
  23 %  You should have received a copy of the GNU General Public License
  24 %  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
  25
  26 Nelem = size(omesh.t,2);
  27 jx = NaN*ones(Nelem,1);
  28 jy = jx;
  29
  30 for iel=1:Nelem
  31
  32         K = sum(omesh.wjacdet(:,iel));
  33
  34         dpsi1 = psi(omesh.t(3,iel))-psi(omesh.t(2,iel));
  35         dvth1 = psith(omesh.t(3,iel))-psith(omesh.t(2,iel));
  36         vthm1 = Utemplogm(psith(omesh.t(3,iel)),psith(omesh.t(2,iel)));
  37         [bp,bn] = Ubern((dpsi1-dvth1)/vthm1);
  38         t1    = omesh.p(:,omesh.t(3,iel))-omesh.p(:,omesh.t(2,iel));
  39         l1    = norm(t1,2);
  40         t1    = t1/l1;
  41         j1    = vthm1*(coeffe(iel)/l1) * ( n(omesh.t(3,iel)) * bp - ...
  42                 n(omesh.t(2,iel)) * bn) * t1;
  43
  44         gg1= omesh.shg(:,2,iel)'*omesh.shg(:,3,iel)*l1^2;
  45
  46         dpsi2 = psi(omesh.t(1,iel))-psi(omesh.t(3,iel));
  47         dvth2 = psith(omesh.t(1,iel))-psith(omesh.t(3,iel));
  48         vthm2 = Utemplogm(psith(omesh.t(1,iel)),psith(omesh.t(3,iel)));
  49         [bp,bn] = Ubern((dpsi2-dvth2)/vthm2);
  50         t2 = omesh.p(:,omesh.t(1,iel))-omesh.p(:,omesh.t(3,iel));
  51         l2 = norm(t2,2);
  52         t2 = t2/l2;
  53         j2 = vthm2*(coeffe(iel)/l2) * ( n(omesh.t(1,iel)) * bp - ...
  54                 n(omesh.t(3,iel)) * bn) * t2;
  55
  56         gg2= omesh.shg(:,1,iel)'*omesh.shg(:,3,iel)*l2^2;
  57
  58         dpsi3 = psi(omesh.t(2,iel))-psi(omesh.t(1,iel));
  59         dvth3 = psith(omesh.t(2,iel))-psith(omesh.t(1,iel));
  60         vthm3 = Utemplogm(psith(omesh.t(2,iel)),psith(omesh.t(1,iel)));
  61         [bp,bn] = Ubern((dpsi3-dvth3)/vthm3);
  62         t3 = omesh.p(:,omesh.t(2,iel))-omesh.p(:,omesh.t(1,iel));
  63         l3 = norm(t3,2);
  64         t3 = t3/l3;
  65         j3 = vthm3*(coeffe(iel)/l3) * ( n(omesh.t(2,iel)) * bp - ...
  66                 n(omesh.t(1,iel)) * bn) * t3;
  67
  68         gg3= omesh.shg(:,2,iel)'*omesh.shg(:,1,iel)*l3^2;
  69
  70         Jk = j1*gg1+j2*gg2+j3*gg3;
  71         jx(iel) = Jk(1);
  72         jy(iel) = Jk(2);
  73
  74 end