Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / secs2d-0.0.8 / Utilities / Uscharfettergummel2.m

diff --git a/octave_packages/secs2d-0.0.8/Utilities/Uscharfettergummel2.m b/octave_packages/secs2d-0.0.8/Utilities/Uscharfettergummel2.m
new file mode 100644 (file)
index 0000000..fada6cf
--- /dev/null
+++ b/octave_packages/secs2d-0.0.8/Utilities/Uscharfettergummel2.m
@@ -0,0 +1,100 @@
+function SG=Uscharfettergummel2(mesh,v,acoeff,bcoeff)
+
+%
+% SG=Ufastscharfettergummel2(mesh,v,acoeff,bcoeff)
+% 
+%
+% Builds the Scharfetter-Gummel  matrix for the 
+% the discretization of the LHS 
+% of the Drift-Diffusion equation:
+%
+% $ -\div (a(x) (\grad (b(x) u) -  b(x) u \grad v'(x) ))= f $
+%
+% where a(x) is piecewise constant
+% and v(x),b(x) is piecewise linear, so that 
+% v'(x) is still piecewise constant
+% and u is the unknown
+%
+
+
+% This file is part of 
+%
+%            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
+%         -------------------------------------------------------------------
+%            Copyright (C) 2004-2006  Carlo de Falco
+%
+%
+%
+%  SECS2D is free software; you can redistribute it and/or modify
+%  it under the terms of the GNU General Public License as published by
+%  the Free Software Foundation; either version 2 of the License, or
+%  (at your option) any later version.
+%
+%  SECS2D is distributed in the hope that it will be useful,
+%  but WITHOUT ANY WARRANTY; without even the implied warranty of
+%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%  GNU General Public License for more details.
+%
+%  You should have received a copy of the GNU General Public License
+%  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
+
+
+Nnodes = length(mesh.p);
+Nelements = length(mesh.t);
+
+areak   = reshape (sum( mesh.wjacdet,1),1,1,Nelements);
+shg     = mesh.shg(:,:,:);
+M       = reshape (acoeff,1,1,Nelements);      
+
+
+% build local Laplacian matrix 
+
+Lloc=zeros(3,3,Nelements);     
+
+for inode=1:3
+  for jnode=1:3
+
+    ginode(inode,jnode,:)=mesh.t(inode,:);
+    gjnode(inode,jnode,:)=mesh.t(jnode,:);
+    Lloc(inode,jnode,:)  = M .* sum( shg(:,inode,:) .* shg(:,jnode,:),1) .* areak;
+
+  end
+end            
+
+vloc    = v(mesh.t(1:3,:));
+bloc    = bcoeff(mesh.t(1:3,:));
+
+blocm1 = Utemplogm(bloc(3,:),bloc(2,:));
+blocm2 = Utemplogm(bloc(1,:),bloc(3,:));
+blocm3 = Utemplogm(bloc(1,:),bloc(2,:));
+
+[bp12,bm12] = Ubern(((vloc(2,:)-vloc(1,:))-(bloc(2,:)-bloc(1,:)))./blocm3);
+[bp13,bm13] = Ubern(((vloc(3,:)-vloc(1,:))-(bloc(3,:)-bloc(1,:)))./blocm2);
+[bp23,bm23] = Ubern(((vloc(3,:)-vloc(2,:))-(bloc(3,:)-bloc(2,:)))./blocm1);
+
+bp12 = reshape(blocm3.*bp12,1,1,Nelements).*Lloc(1,2,:);
+bm12 = reshape(blocm3.*bm12,1,1,Nelements).*Lloc(1,2,:);
+bp13 = reshape(blocm2.*bp13,1,1,Nelements).*Lloc(1,3,:);
+bm13 = reshape(blocm2.*bm13,1,1,Nelements).*Lloc(1,3,:);
+bp23 = reshape(blocm1.*bp23,1,1,Nelements).*Lloc(2,3,:);
+bm23 = reshape(blocm1.*bm23,1,1,Nelements).*Lloc(2,3,:);
+
+SGloc(1,1,:) = -bm12-bm13;
+SGloc(1,2,:) = bp12;
+SGloc(1,3,:) = bp13;
+
+SGloc(2,1,:) = bm12;
+SGloc(2,2,:) = -bp12-bm23; 
+SGloc(2,3,:) = bp23;
+
+SGloc(3,1,:) = bm13;
+SGloc(3,2,:) = bm23;
+SGloc(3,3,:) = -bp13-bp23;
+
+##SGloc=[-bm12-bm13,   bp12     ,   bp13
+##       bm12     ,  -bp12-bm23,   bp23
+##        bm13     ,   bm23     ,  -bp13-bp23];
+
+SG = sparse(ginode(:),gjnode(:),SGloc(:));
+
+