octave_packages/secs2d-0.0.8/Utilities/Udriftdiffusion.m

   1 function c=Udriftdiffusion(mesh,Dsides,guess,M,U,V,u)
   2
   3 %%
   4 %  c=Udriftdiffusion(mesh,Dsides,guess,M,U,V,u)
   5 % solves the drift diffusion equation
   6 % $ -div ( u ( \nabla n - 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
  33 if (Ucolumns(guess)>Urows(guess))
  34     guess=guess';
  35 end
  36
  37 if (Ucolumns(V)>Urows(V))
  38     V=V';
  39 end
  40
  41 if (Ucolumns(U)>Urows(U))
  42     U=U';
  43 end
  44 Nnodes    = max(size(mesh.p));
  45 Nelements = max(size(mesh.t));
  46
  47
  48 % Set list of nodes with Dirichelet BCs
  49 Dnodes=Unodesonside(mesh,Dsides);
  50
  51 % Set values of Dirichelet BCs
  52 Bc     = guess(Dnodes);
  53 % Set list of nodes without Dirichelet BCs
  54 Varnodes = setdiff([1:Nnodes],Dnodes);
  55
  56
  57 % Build LHS matrix and RHS
  58 A = Uscharfettergummel(mesh,V,u);
  59 if (isempty(DDG_MASS))
  60   DDG_MASS=Ucompmass2(mesh,ones(Nnodes,1),ones(Nelements,1));
  61 end
  62 A = A + DDG_MASS.*spdiag(M,0);
  63
  64 if (isempty(DDG_RHS))
  65         DDG_RHS=Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1));
  66 end
  67 b = DDG_RHS.*U;
  68
  69
  70
  71 %%%%%%%%DEBUG%%%%%%%%%%%
  72 %DEBUG_SGM = A;
  73
  74 %% Apply boundary conditions
  75 A (Dnodes,:) = 0;
  76 b (Dnodes)   = 0;
  77 b = b - A (:,Dnodes) * Bc;
  78
  79 A(Dnodes,:)= [];
  80 A(:,Dnodes)= [];
  81
  82 b(Dnodes)       = [];
  83
  84 % Enforce positivity
  85 % dA = abs(diag(A));
  86 % edA= abs(sum(triu(A,1)+tril(A,-1),1));
  87 % nondominant = find(dA < edA');
  88 % for ii=1:length(nondominant)
  89 %     A(nondominant(ii),nondominant(ii)) = edA(ii);
  90 % end
  91
  92 c = guess;
  93
  94 c(Varnodes) = A \ b;
  95