octave_packages/secs2d-0.0.8/Utilities/Uscharfettergummel2.m

   1 function SG=Uscharfettergummel2(mesh,v,acoeff,bcoeff)
   2
   3 %
   4 % SG=Ufastscharfettergummel2(mesh,v,acoeff,bcoeff)
   5 %
   6 %
   7 % Builds the Scharfetter-Gummel  matrix for the
   8 % the discretization of the LHS
   9 % of the Drift-Diffusion equation:
  10 %
  11 % $ -\div (a(x) (\grad (b(x) u) -  b(x) u \grad v'(x) ))= f $
  12 %
  13 % where a(x) is piecewise constant
  14 % and v(x),b(x) is piecewise linear, so that
  15 % v'(x) is still piecewise constant
  16 % and u is the unknown
  17 %
  18
  19
  20 % This file is part of
  21 %
  22 %            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
  23 %         -------------------------------------------------------------------
  24 %            Copyright (C) 2004-2006  Carlo de Falco
  25 %
  26 %
  27 %
  28 %  SECS2D is free software; you can redistribute it and/or modify
  29 %  it under the terms of the GNU General Public License as published by
  30 %  the Free Software Foundation; either version 2 of the License, or
  31 %  (at your option) any later version.
  32 %
  33 %  SECS2D is distributed in the hope that it will be useful,
  34 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
  35 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  36 %  GNU General Public License for more details.
  37 %
  38 %  You should have received a copy of the GNU General Public License
  39 %  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
  40
  41
  42 Nnodes = length(mesh.p);
  43 Nelements = length(mesh.t);
  44
  45 areak   = reshape (sum( mesh.wjacdet,1),1,1,Nelements);
  46 shg     = mesh.shg(:,:,:);
  47 M       = reshape (acoeff,1,1,Nelements);
  48
  49
  50 % build local Laplacian matrix
  51
  52 Lloc=zeros(3,3,Nelements);
  53
  54 for inode=1:3
  55   for jnode=1:3
  56
  57     ginode(inode,jnode,:)=mesh.t(inode,:);
  58     gjnode(inode,jnode,:)=mesh.t(jnode,:);
  59     Lloc(inode,jnode,:)  = M .* sum( shg(:,inode,:) .* shg(:,jnode,:),1) .* areak;
  60
  61   end
  62 end
  63
  64 vloc    = v(mesh.t(1:3,:));
  65 bloc    = bcoeff(mesh.t(1:3,:));
  66
  67 blocm1 = Utemplogm(bloc(3,:),bloc(2,:));
  68 blocm2 = Utemplogm(bloc(1,:),bloc(3,:));
  69 blocm3 = Utemplogm(bloc(1,:),bloc(2,:));
  70
  71 [bp12,bm12] = Ubern(((vloc(2,:)-vloc(1,:))-(bloc(2,:)-bloc(1,:)))./blocm3);
  72 [bp13,bm13] = Ubern(((vloc(3,:)-vloc(1,:))-(bloc(3,:)-bloc(1,:)))./blocm2);
  73 [bp23,bm23] = Ubern(((vloc(3,:)-vloc(2,:))-(bloc(3,:)-bloc(2,:)))./blocm1);
  74
  75 bp12 = reshape(blocm3.*bp12,1,1,Nelements).*Lloc(1,2,:);
  76 bm12 = reshape(blocm3.*bm12,1,1,Nelements).*Lloc(1,2,:);
  77 bp13 = reshape(blocm2.*bp13,1,1,Nelements).*Lloc(1,3,:);
  78 bm13 = reshape(blocm2.*bm13,1,1,Nelements).*Lloc(1,3,:);
  79 bp23 = reshape(blocm1.*bp23,1,1,Nelements).*Lloc(2,3,:);
  80 bm23 = reshape(blocm1.*bm23,1,1,Nelements).*Lloc(2,3,:);
  81
  82 SGloc(1,1,:) = -bm12-bm13;
  83 SGloc(1,2,:) = bp12;
  84 SGloc(1,3,:) = bp13;
  85
  86 SGloc(2,1,:) = bm12;
  87 SGloc(2,2,:) = -bp12-bm23;
  88 SGloc(2,3,:) = bp23;
  89
  90 SGloc(3,1,:) = bm13;
  91 SGloc(3,2,:) = bm23;
  92 SGloc(3,3,:) = -bp13-bp23;
  93
  94 ##SGloc=[-bm12-bm13,   bp12     ,   bp13
  95 ##       bm12     ,  -bp12-bm23,   bp23
  96 ##        bm13     ,   bm23     ,  -bp13-bp23];
  97
  98 SG = sparse(ginode(:),gjnode(:),SGloc(:));
  99
 100