1function S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta) 2 3 ## -*- texinfo -*- 4 ## 5 ## @deftypefn {Function File} @ 6 ## {@var{S}} = Uscharfettergummel3 (@var{mesh}, @var{alpha}, @ 7 ## @var{gamma}, @var{eta}, @var{beta}) 8 ## 9 ## 10 ## Builds the Scharfetter-Gummel matrix for the 11 ## discretization of the LHS 12 ## of the equation: 13 ## 14 ## @iftex 15 ## @tex 16 ## $ -div ( \alpha \gamma ( \eta \vect{\nabla} u - \vect{beta} u )) = f $ 17 ## @end tex 18 ## @end iftex 19 ## @ifinfo 20 ## -div (@var{alpha} * @var{gamma} (@var{eta} grad u - @var{beta} u )) = f 21 ## @end ifinfo 22 ## 23 ## where: 24 ## @itemize @minus 25 ## @item @var{alpha} is an element-wise constant scalar function 26 ## @item @var{eta}, @var{gamma} are piecewise linear conforming 27 ## scalar functions 28 ## @item @var{beta} is an element-wise constant vector function 29 ## @end itemize 30 ## 31 ## Instead of passing the vector field @var{beta} directly 32 ## one can pass a piecewise linear conforming scalar function 33 ## @var{phi} as the last input. In such case @var{beta} = grad @var{phi} 34 ## is assumed. If @var{phi} is a single scalar value @var{beta} 35 ## is assumed to be 0 in the whole domain. 36 ## 37 ## Example: 38 ## @example 39 ## [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:1/3:1],[0:1/3:1],1,1:4); 40 ## mesh = Umeshproperties(mesh); 41 ## x = mesh.p(1,:)'; 42 ## Dnodes = Unodesonside(mesh,[2,4]); 43 ## Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 44 ## Varnodes = setdiff(1:Nnodes,Dnodes); 45 ## alpha = ones(Nelements,1); eta = .1*ones(Nnodes,1); 46 ## beta = [ones(1,Nelements);zeros(1,Nelements)]; 47 ## gamma = ones(Nnodes,1); 48 ## f = Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1)); 49 ## S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 50 ## u = zeros(Nnodes,1); 51 ## u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 52 ## uex = x - (exp(10*x)-1)/(exp(10)-1); 53 ## assert(u,uex,1e-7) 54 ## @end example 55 ## 56 ## @seealso{Ucomplap, Ucompconst, Ucompmass2, Uscharfettergummel} 57 ## @end deftypefn 58 59 60 %% This file is part of 61 %% 62 %% SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator 63 %% ------------------------------------------------------------------- 64 %% Copyright (C) 2004-2006 Carlo de Falco 65 %% 66 %% 67 %% 68 %% SECS2D is free software; you can redistribute it and/or modify 69 %% it under the terms of the GNU General Public License as published by 70 %% the Free Software Foundation; either version 2 of the License, or 71 %% (at your option) any later version. 72 %% 73 %% SECS2D is distributed in the hope that it will be useful, 74 %% but WITHOUT ANY WARRANTY; without even the implied warranty of 75 %% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 76 %% GNU General Public License for more details. 77 %% 78 %% You should have received a copy of the GNU General Public License 79 %% along with SECS2D; If not, see <http://www.gnu.org/licenses/>. 80 81 82 Nnodes = columns(mesh.p); 83 Nelements = columns(mesh.t); 84 85 alphaareak = reshape (alpha.*sum( mesh.wjacdet,1)',1,1,Nelements); 86 shg = mesh.shg(:,:,:); 87 88 89 % build local Laplacian matrix 90 91 Lloc=zeros(3,3,Nelements); 92 93 for inode=1:3 94 for jnode=1:3 95 ginode(inode,jnode,:)=mesh.t(inode,:); 96 gjnode(inode,jnode,:)=mesh.t(jnode,:); 97 Lloc(inode,jnode,:) = sum( shg(:,inode,:) .* shg(:,jnode,:),1)... 98 .* alphaareak; 99 end 100 end 101 102 103 x = mesh.p(1,:); 104 x = x(mesh.t(1:3,:)); 105 y = mesh.p(2,:); 106 y = y(mesh.t(1:3,:)); 107 108 if all(size(beta)==1) 109 v12=0;v23=0;v31=0; 110 elseif all(size(beta)==[2,Nelements]) 111 v12 = beta(1,:) .* (x(2,:)-x(1,:)) + beta(2,:) .* (y(2,:)-y(1,:)); 112 v23 = beta(1,:) .* (x(3,:)-x(2,:)) + beta(2,:) .* (y(3,:)-y(2,:)); 113 v31 = beta(1,:) .* (x(1,:)-x(3,:)) + beta(2,:) .* (y(1,:)-y(3,:)); 114 elseif all(size(beta)==[Nnodes,1]) 115 betaloc = beta(mesh.t(1:3,:)); 116 v12 = betaloc(2,:)-betaloc(1,:); 117 v23 = betaloc(3,:)-betaloc(2,:); 118 v31 = betaloc(1,:)-betaloc(3,:); 119 end 120 121 etaloc = eta(mesh.t(1:3,:)); 122 123 eta12 = etaloc(2,:)-etaloc(1,:); 124 eta23 = etaloc(3,:)-etaloc(2,:); 125 eta31 = etaloc(1,:)-etaloc(3,:); 126 127 etalocm1 = Utemplogm(etaloc(2,:),etaloc(3,:)); 128 etalocm2 = Utemplogm(etaloc(3,:),etaloc(1,:)); 129 etalocm3 = Utemplogm(etaloc(1,:),etaloc(2,:)); 130 131 gammaloc = gamma(mesh.t(1:3,:)); 132 geloc = gammaloc.*etaloc; 133 134 gelocm1 = Utemplogm(geloc(2,:),geloc(3,:)); 135 gelocm2 = Utemplogm(geloc(3,:),geloc(1,:)); 136 gelocm3 = Utemplogm(geloc(1,:),geloc(2,:)); 137 138 [bp12,bm12] = Ubern( (v12 - eta12)./etalocm3); 139 [bp23,bm23] = Ubern( (v23 - eta23)./etalocm1); 140 [bp31,bm31] = Ubern( (v31 - eta31)./etalocm2); 141 142 bp12 = reshape(gelocm3.*etalocm3.*bp12,1,1,Nelements).*Lloc(1,2,:); 143 bm12 = reshape(gelocm3.*etalocm3.*bm12,1,1,Nelements).*Lloc(1,2,:); 144 bp23 = reshape(gelocm1.*etalocm1.*bp23,1,1,Nelements).*Lloc(2,3,:); 145 bm23 = reshape(gelocm1.*etalocm1.*bm23,1,1,Nelements).*Lloc(2,3,:); 146 bp31 = reshape(gelocm2.*etalocm2.*bp31,1,1,Nelements).*Lloc(3,1,:); 147 bm31 = reshape(gelocm2.*etalocm2.*bm31,1,1,Nelements).*Lloc(3,1,:); 148 149 Sloc(1,1,:) = (-bm12-bp31)./reshape(etaloc(1,:),1,1,Nelements); 150 Sloc(1,2,:) = bp12./reshape(etaloc(2,:),1,1,Nelements); 151 Sloc(1,3,:) = bm31./reshape(etaloc(3,:),1,1,Nelements); 152 153 Sloc(2,1,:) = bm12./reshape(etaloc(1,:),1,1,Nelements); 154 Sloc(2,2,:) = (-bp12-bm23)./reshape(etaloc(2,:),1,1,Nelements); 155 Sloc(2,3,:) = bp23./reshape(etaloc(3,:),1,1,Nelements); 156 157 Sloc(3,1,:) = bp31./reshape(etaloc(1,:),1,1,Nelements); 158 Sloc(3,2,:) = bm23./reshape(etaloc(2,:),1,1,Nelements); 159 Sloc(3,3,:) = (-bm31-bp23)./reshape(etaloc(3,:),1,1,Nelements); 160 161 S = sparse(ginode(:),gjnode(:),Sloc(:)); 162 163endfunction 164 165%!test 166%! [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:1/3:1],[0:1/3:1],1,1:4); 167%! mesh = Umeshproperties(mesh); 168%! x = mesh.p(1,:)'; 169%! Dnodes = Unodesonside(mesh,[2,4]); 170%! Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 171%! Varnodes = setdiff(1:Nnodes,Dnodes); 172%! alpha = ones(Nelements,1); eta = .1*ones(Nnodes,1); 173%! beta = [ones(1,Nelements);zeros(1,Nelements)]; 174%! gamma = ones(Nnodes,1); 175%! f = Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1)); 176%! S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 177%! u = zeros(Nnodes,1); 178%! u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 179%! uex = x - (exp(10*x)-1)/(exp(10)-1); 180%! assert(u,uex,1e-7) 181 182%!test 183%! [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:1/3:1],[0:1/3:1],1,1:4); 184%! mesh = Umeshproperties(mesh); 185%! x = mesh.p(1,:)'; 186%! Dnodes = Unodesonside(mesh,[2,4]); 187%! Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 188%! Varnodes = setdiff(1:Nnodes,Dnodes); 189%! alpha = ones(Nelements,1); eta = .1*ones(Nnodes,1); 190%! beta = x; 191%! gamma = ones(Nnodes,1); 192%! f = Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1)); 193%! S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 194%! u = zeros(Nnodes,1); 195%! u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 196%! uex = x - (exp(10*x)-1)/(exp(10)-1); 197%! assert(u,uex,1e-7) 198 199%!test 200%! [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:1/3:1],[0:1/3:1],1,1:4); 201%! mesh = Umeshproperties(mesh); 202%! x = mesh.p(1,:)'; 203%! Dnodes = Unodesonside(mesh,[2,4]); 204%! Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 205%! Varnodes = setdiff(1:Nnodes,Dnodes); 206%! alpha = 10*ones(Nelements,1); eta = .01*ones(Nnodes,1); 207%! beta = x/10; 208%! gamma = ones(Nnodes,1); 209%! f = Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1)); 210%! S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 211%! u = zeros(Nnodes,1); 212%! u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 213%! uex = x - (exp(10*x)-1)/(exp(10)-1); 214%! assert(u,uex,1e-7) 215 216%!test 217%! [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:1/3:1],[0:1/3:1],1,1:4); 218%! mesh = Umeshproperties(mesh); 219%! x = mesh.p(1,:)'; 220%! Dnodes = Unodesonside(mesh,[2,4]); 221%! Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 222%! Varnodes = setdiff(1:Nnodes,Dnodes); 223%! alpha = 10*ones(Nelements,1); eta = .001*ones(Nnodes,1); 224%! beta = x/100; 225%! gamma = 10*ones(Nnodes,1); 226%! f = Ucompconst(mesh,ones(Nnodes,1),ones(Nelements,1)); 227%! S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 228%! u = zeros(Nnodes,1); 229%! u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 230%! uex = x - (exp(10*x)-1)/(exp(10)-1); 231%! assert(u,uex,1e-7) 232 233%!test 234%! [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:1/1e3:1],[0:1/2:1],1,1:4); 235%! mesh = Umeshproperties(mesh); 236%! x = mesh.p(1,:)'; 237%! Dnodes = Unodesonside(mesh,[2,4]); 238%! Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 239%! Varnodes = setdiff(1:Nnodes,Dnodes); 240%! alpha = 3*ones(Nelements,1); eta = x+1; 241%! beta = [ones(1,Nelements);zeros(1,Nelements)]; 242%! gamma = 2*x; 243%! ff = 2*(6*x.^2+6*x) - (6*x+6).*(1-2*x)+6*(x-x.^2); 244%! f = Ucompconst(mesh,ff,ones(Nelements,1)); 245%! S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 246%! u = zeros(Nnodes,1); 247%! u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 248%! uex = x - x.^2; 249%! assert(u,uex,5e-3) 250 251%!test 252%! [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:1/1e3:1],[0:1/2:1],1,1:4); 253%! mesh = Umeshproperties(mesh); 254%! x = mesh.p(1,:)'; 255%! Dnodes = Unodesonside(mesh,[2,4]); 256%! Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 257%! Varnodes = setdiff(1:Nnodes,Dnodes); 258%! alpha = ones(Nelements,1); eta = ones(Nnodes,1); 259%! beta = 0; 260%! gamma = x+1; 261%! ff = 4*x+1; 262%! f = Ucompconst(mesh,ff,ones(Nelements,1)); 263%! S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 264%! u = zeros(Nnodes,1); 265%! u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 266%! uex = x - x.^2; 267%! assert(u,uex,1e-7) 268 269%!test 270%! [mesh.p,mesh.e,mesh.t] = Ustructmesh([0:.1:1],[0:.1:1],1,1:4); 271%! mesh = Umeshproperties(mesh); 272%! x = mesh.p(1,:)';y = mesh.p(2,:)'; 273%! Dnodes = Unodesonside(mesh,[1:4]); 274%! Nnodes = columns(mesh.p); Nelements = columns(mesh.t); 275%! Varnodes = setdiff(1:Nnodes,Dnodes); 276%! alpha = ones(Nelements,1); diff = 1e-2; eta=diff*ones(Nnodes,1); 277%! beta =[ones(1,Nelements);ones(1,Nelements)]; 278%! gamma = x*0+1; 279%! ux = y.*(1-exp((y-1)/diff)) .* (1-exp((x-1)/diff)-x.*exp((x-1)/diff)/diff); 280%! uy = x.*(1-exp((x-1)/diff)) .* (1-exp((y-1)/diff)-y.*exp((y-1)/diff)/diff); 281%! uxx = y.*(1-exp((y-1)/diff)) .* (-2*exp((x-1)/diff)/diff-x.*exp((x-1)/diff)/(diff^2)); 282%! uyy = x.*(1-exp((x-1)/diff)) .* (-2*exp((y-1)/diff)/diff-y.*exp((y-1)/diff)/(diff^2)); 283%! ff = -diff*(uxx+uyy)+ux+uy; 284%! f = Ucompconst(mesh,ff,ones(Nelements,1)); 285%! S = Uscharfettergummel3(mesh,alpha,gamma,eta,beta); 286%! u = zeros(Nnodes,1); 287%! u(Varnodes) = S(Varnodes,Varnodes)\f(Varnodes); 288%! uex = x.*y.*(1-exp((x-1)/diff)).*(1-exp((y-1)/diff)); 289%! assert(u,uex,1e-7) 290 291