octave-secs2d 0.0.8-5 / usr / share / octave / packages / secs2d-0.0.8 / QDDGOX / QDDGOXcompdens.m

function w = QDDGOXcompdens(mesh,Dsides,win,vin,fermiin,d2,toll,maxit,verbose);

%  w = QDDGOXcompdens(mesh,Dsides,win,vin,fermiin,d2,toll,maxit,verbose);

global QDDGOXCOMPDENS_LAP QDDGOXCOMPDENS_MASS QDDGOXCOMPDENS_RHS 
%% Set some usefull constants
VErank = 4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% convert input vectors to columns
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if Ucolumns(win)>Urows(win)
    win=win';
end
if Ucolumns(vin)>Urows(vin)
    vin=vin';
end 
if Ucolumns(fermiin)>Urows(fermiin)
    fermiin=fermiin';
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% convert grid info to FEM form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nodes 		    = mesh.p;
Nnodes		    = size(nodes,2);

elements	    = mesh.t(1:3,:);
Nelements           = size(elements,2);

Dedges    =[];

for ii = 1:length(Dsides)
	Dedges=[Dedges,find(mesh.e(5,:)==Dsides(ii))];
end

% Set list of nodes with Dirichelet BCs
Dnodes = mesh.e(1:2,Dedges);
Dnodes = [Dnodes(1,:) Dnodes(2,:)];
Dnodes = unique(Dnodes);

Dvals  = win(Dnodes);

Varnodes = setdiff([1:Nnodes],Dnodes);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% 		initialization:
%% 		we're going to solve
%% 		$$ -\delta^2 \Lap w_{k+1} + B'(w_k) \delta w_{k+1} =  2 * w_k$$
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% set $$ w_1 = win $$ 
w = win;
wnew = win;

%% let's compute  FEM approximation of $$ L = - \aleph \frac{d^2}{x^2} $$
if (isempty(QDDGOXCOMPDENS_LAP))
    QDDGOXCOMPDENS_LAP = Ucomplap (mesh,ones(Nelements,1));
end
L = d2*QDDGOXCOMPDENS_LAP;

%% now compute $$ G_k = F - V  + 2 V_{th} log(w) $$
if (isempty(QDDGOXCOMPDENS_MASS))
    QDDGOXCOMPDENS_MASS = Ucompmass2 (mesh,ones(Nnodes,1),ones(Nelements,1));
end
G	    = fermiin - vin  + 2*log(w);
Bmat	= QDDGOXCOMPDENS_MASS*sparse(diag(G));
nrm     = 1;
%%%%%%%%%%%%%%%%%%%%%%%%
%%% NEWTON ITERATION START
%%%%%%%%%%%%%%%%%%%%%%%%
converged = 0;
for jnewt =1:ceil(maxit/VErank)
  for k=1:VErank
    [w(:,k+1),converged,G,L,Bmat]=onenewtit(w(:,k),G,fermiin,vin,L,Bmat,jnewt,mesh,Dnodes,Varnodes,Dvals,Nnodes,Nelements,toll);        
    if converged
      break
    end
  end
  if converged
    break
  end
  w  = Urrextrapolation(w);
end	
%%%%%%%%%%%%%%%%%%%%%%%%
%%% NEWTON ITERATION END
%%%%%%%%%%%%%%%%%%%%%%%%
w = w(:,end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% ONE NEWTON ITERATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [w,converged,G,L,Bmat]=onenewtit(w,G,fermiin,vin,L,Bmat,jnewt,mesh,Dnodes,Varnodes,Dvals,Nnodes,Nelements,toll);
  
  global QDDGOXCOMPDENS_LAP QDDGOXCOMPDENS_MASS QDDGOXCOMPDENS_RHS 
  dampit 		= 5;
  dampcoeff	        = 2;
  converged             = 0;
  wnew                  =  w;
  
  res0      = norm((L + Bmat) * w,inf);
  
  
  %% chose $ t_k $ to ensure positivity of  $w$
  mm  = -min(G);
  pause(1)

  if (mm>2)
    tk = max( 1/(mm));
  else
    tk = 1;
  end

  tmpmat = QDDGOXCOMPDENS_MASS*2;
  if (isempty(QDDGOXCOMPDENS_RHS))
    QDDGOXCOMPDENS_RHS = Ucompconst (mesh,ones(Nnodes,1),ones(Nelements,1));
  end
  tmpvect= 2*QDDGOXCOMPDENS_RHS.*w;
  
  %%%%%%%%%%%%%%%%%%%%%%%%
  %%% DAMPING ITERATION START
  %%%%%%%%%%%%%%%%%%%%%%%%
  for idamp = 1:dampit
    %% Compute $ B1mat = \frac{2}{t_k} $
    %% and the (lumped) mass matrix B1mat(w_k)
    B1mat	= tmpmat/tk;
    
    %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
    A 		= L + B1mat + Bmat;
    b 		= tmpvect/tk;
    
    %% Apply boundary conditions
    A (Dnodes,:) = 0;
    b (Dnodes)   = 0;
    b = b - A (:,Dnodes) * Dvals;
    
    A(Dnodes,:)= [];
    A(:,Dnodes)= [];
    
    b(Dnodes)	= [];
    

    wnew(Varnodes) =  A\b;	
    
    
    %% compute $$ G_{k+1} = F - V  + 2 V_{th} log(w) $$
    G	    = fermiin - vin  + 2*log(wnew);
    Bmat	= QDDGOXCOMPDENS_MASS*sparse(diag(G));
    
    res    = norm((L + Bmat) * wnew,inf);
    
    if (res<res0)
      break
    else
      tk = tk/dampcoeff;
    end
  end	
  %%%%%%%%%%%%%%%%%%%%%%%%
  %%% DAMPING ITERATION END
  %%%%%%%%%%%%%%%%%%%%%%%%
  nrm = norm(wnew-w);
  
  if (nrm < toll)
    w= wnew;
    converged = 1;
  else
    w=wnew;
  end