octave_packages/secs2d-0.0.8/Utilities/Udriftdiffusion2.m

   1 function c=Udriftdiffusion2(mesh,Dsides,guess,M,U,V,Vth,u)
   2
   3 %%
   4 %  c=Udriftdiffusion(mesh,Dsides,guess,M,U,V,Vth,u)
   5 % solves the drift diffusion equation
   6 % $ -div ( u ( \nabla (n Vth) - n \nabla V)) + M = U $
   7 %%
   8
   9 % This file is part of
  10 %
  11 %            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
  12 %         -------------------------------------------------------------------
  13 %            Copyright (C) 2004-2006  Carlo de Falco
  14 %
  15 %
  16 %
  17 %  SECS2D is free software; you can redistribute it and/or modify
  18 %  it under the terms of the GNU General Public License as published by
  19 %  the Free Software Foundation; either version 2 of the License, or
  20 %  (at your option) any later version.
  21 %
  22 %  SECS2D is distributed in the hope that it will be useful,
  23 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25 %  GNU General Public License for more details.
  26 %
  27 %  You should have received a copy of the GNU General Public License
  28 %  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
  29
  30 global DDG_RHS DDG_MASS %DEBUG_SGM
  31
  32 Nnodes    = max(size(mesh.p));
  33 Nelements = max(size(mesh.t));
  34
  35
  36 % Set list of nodes with Dirichelet BCs
  37 Dnodes=Unodesonside(mesh,Dsides);
  38
  39 % Set values of Dirichelet BCs
  40 Bc     = guess(Dnodes);
  41
  42 % Set list of nodes without Dirichelet BCs
  43 Varnodes = setdiff([1:Nnodes],Dnodes);
  44
  45 % Build LHS matrix and RHS
  46 A = Uscharfettergummel2(mesh,V,u,Vth);
  47 if (isempty(DDG_MASS))
  48   DDG_MASS=Ucompmass2(mesh,ones(Nnodes,1),ones(Nelements,1));
  49 end
  50 A = A + DDG_MASS.*spdiag(M,0);
  51
  52 if (isempty(DDG_RHS))
  53         DDG_RHS=Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1));
  54 end
  55 b = DDG_RHS.*U;
  56
  57
  58
  59 %%%%%%%%DEBUG%%%%%%%%%%%
  60 %DEBUG_SGM = A;
  61
  62 %% Apply boundary conditions
  63 A (Dnodes,:) = 0;
  64 b (Dnodes)   = 0;
  65 b = b - A (:,Dnodes) * Bc;
  66
  67 A(Dnodes,:)= [];
  68 A(:,Dnodes)= [];
  69
  70 b(Dnodes)       = [];
  71
  72
  73 c = guess;
  74 c(Varnodes) = A \ b;
  75