octave_packages/secs2d-0.0.8/QDDGOX/QDDGOXcompdens.m

   1 function w = QDDGOXcompdens(mesh,Dsides,win,vin,fermiin,d2,toll,maxit,verbose);
   2
   3 %  w = QDDGOXcompdens(mesh,Dsides,win,vin,fermiin,d2,toll,maxit,verbose);
   4
   5 global QDDGOXCOMPDENS_LAP QDDGOXCOMPDENS_MASS QDDGOXCOMPDENS_RHS
   6 %% Set some usefull constants
   7 VErank = 4;
   8
   9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  10 %% convert input vectors to columns
  11 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  12
  13 if Ucolumns(win)>Urows(win)
  14     win=win';
  15 end
  16 if Ucolumns(vin)>Urows(vin)
  17     vin=vin';
  18 end
  19 if Ucolumns(fermiin)>Urows(fermiin)
  20     fermiin=fermiin';
  21 end
  22
  23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  24 %% convert grid info to FEM form
  25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  26
  27 nodes               = mesh.p;
  28 Nnodes              = size(nodes,2);
  29
  30 elements            = mesh.t(1:3,:);
  31 Nelements           = size(elements,2);
  32
  33 Dedges    =[];
  34
  35 for ii = 1:length(Dsides)
  36         Dedges=[Dedges,find(mesh.e(5,:)==Dsides(ii))];
  37 end
  38
  39 % Set list of nodes with Dirichelet BCs
  40 Dnodes = mesh.e(1:2,Dedges);
  41 Dnodes = [Dnodes(1,:) Dnodes(2,:)];
  42 Dnodes = unique(Dnodes);
  43
  44 Dvals  = win(Dnodes);
  45
  46 Varnodes = setdiff([1:Nnodes],Dnodes);
  47 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  48 %%
  49 %%              initialization:
  50 %%              we're going to solve
  51 %%              $$ -\delta^2 \Lap w_{k+1} + B'(w_k) \delta w_{k+1} =  2 * w_k$$
  52 %%
  53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  54
  55 %% set $$ w_1 = win $$
  56 w = win;
  57 wnew = win;
  58
  59 %% let's compute  FEM approximation of $$ L = - \aleph \frac{d^2}{x^2} $$
  60 if (isempty(QDDGOXCOMPDENS_LAP))
  61     QDDGOXCOMPDENS_LAP = Ucomplap (mesh,ones(Nelements,1));
  62 end
  63 L = d2*QDDGOXCOMPDENS_LAP;
  64
  65 %% now compute $$ G_k = F - V  + 2 V_{th} log(w) $$
  66 if (isempty(QDDGOXCOMPDENS_MASS))
  67     QDDGOXCOMPDENS_MASS = Ucompmass2 (mesh,ones(Nnodes,1),ones(Nelements,1));
  68 end
  69 G           = fermiin - vin  + 2*log(w);
  70 Bmat    = QDDGOXCOMPDENS_MASS*sparse(diag(G));
  71 nrm     = 1;
  72 %%%%%%%%%%%%%%%%%%%%%%%%
  73 %%% NEWTON ITERATION START
  74 %%%%%%%%%%%%%%%%%%%%%%%%
  75 converged = 0;
  76 for jnewt =1:ceil(maxit/VErank)
  77   for k=1:VErank
  78     [w(:,k+1),converged,G,L,Bmat]=onenewtit(w(:,k),G,fermiin,vin,L,Bmat,jnewt,mesh,Dnodes,Varnodes,Dvals,Nnodes,Nelements,toll);
  79     if converged
  80       break
  81     end
  82   end
  83   if converged
  84     break
  85   end
  86   w  = Urrextrapolation(w);
  87 end
  88 %%%%%%%%%%%%%%%%%%%%%%%%
  89 %%% NEWTON ITERATION END
  90 %%%%%%%%%%%%%%%%%%%%%%%%
  91 w = w(:,end);
  92
  93 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  94 %%%%%%%%% ONE NEWTON ITERATION
  95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  96 function [w,converged,G,L,Bmat]=onenewtit(w,G,fermiin,vin,L,Bmat,jnewt,mesh,Dnodes,Varnodes,Dvals,Nnodes,Nelements,toll);
  97
  98   global QDDGOXCOMPDENS_LAP QDDGOXCOMPDENS_MASS QDDGOXCOMPDENS_RHS
  99   dampit                = 5;
 100   dampcoeff             = 2;
 101   converged             = 0;
 102   wnew                  =  w;
 103
 104   res0      = norm((L + Bmat) * w,inf);
 105
 106
 107   %% chose $ t_k $ to ensure positivity of  $w$
 108   mm  = -min(G);
 109   pause(1)
 110
 111   if (mm>2)
 112     tk = max( 1/(mm));
 113   else
 114     tk = 1;
 115   end
 116
 117   tmpmat = QDDGOXCOMPDENS_MASS*2;
 118   if (isempty(QDDGOXCOMPDENS_RHS))
 119     QDDGOXCOMPDENS_RHS = Ucompconst (mesh,ones(Nnodes,1),ones(Nelements,1));
 120   end
 121   tmpvect= 2*QDDGOXCOMPDENS_RHS.*w;
 122
 123   %%%%%%%%%%%%%%%%%%%%%%%%
 124   %%% DAMPING ITERATION START
 125   %%%%%%%%%%%%%%%%%%%%%%%%
 126   for idamp = 1:dampit
 127     %% Compute $ B1mat = \frac{2}{t_k} $
 128     %% and the (lumped) mass matrix B1mat(w_k)
 129     B1mat       = tmpmat/tk;
 130
 131     %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
 132     A           = L + B1mat + Bmat;
 133     b           = tmpvect/tk;
 134
 135     %% Apply boundary conditions
 136     A (Dnodes,:) = 0;
 137     b (Dnodes)   = 0;
 138     b = b - A (:,Dnodes) * Dvals;
 139
 140     A(Dnodes,:)= [];
 141     A(:,Dnodes)= [];
 142
 143     b(Dnodes)   = [];
 144
 145
 146     wnew(Varnodes) =  A\b;
 147
 148
 149     %% compute $$ G_{k+1} = F - V  + 2 V_{th} log(w) $$
 150     G       = fermiin - vin  + 2*log(wnew);
 151     Bmat        = QDDGOXCOMPDENS_MASS*sparse(diag(G));
 152
 153     res    = norm((L + Bmat) * wnew,inf);
 154
 155     if (res<res0)
 156       break
 157     else
 158       tk = tk/dampcoeff;
 159     end
 160   end
 161   %%%%%%%%%%%%%%%%%%%%%%%%
 162   %%% DAMPING ITERATION END
 163   %%%%%%%%%%%%%%%%%%%%%%%%
 164   nrm = norm(wnew-w);
 165
 166   if (nrm < toll)
 167     w= wnew;
 168     converged = 1;
 169   else
 170     w=wnew;
 171   end
 172
 173