1 function w = QDDGOXcompdens(mesh,Dsides,win,vin,fermiin,d2,toll,maxit,verbose);
3 % w = QDDGOXcompdens(mesh,Dsides,win,vin,fermiin,d2,toll,maxit,verbose);
5 global QDDGOXCOMPDENS_LAP QDDGOXCOMPDENS_MASS QDDGOXCOMPDENS_RHS
6 %% Set some usefull constants
9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10 %% convert input vectors to columns
11 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13 if Ucolumns(win)>Urows(win)
16 if Ucolumns(vin)>Urows(vin)
19 if Ucolumns(fermiin)>Urows(fermiin)
23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24 %% convert grid info to FEM form
25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 Nnodes = size(nodes,2);
30 elements = mesh.t(1:3,:);
31 Nelements = size(elements,2);
35 for ii = 1:length(Dsides)
36 Dedges=[Dedges,find(mesh.e(5,:)==Dsides(ii))];
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);
46 Varnodes = setdiff([1:Nnodes],Dnodes);
47 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50 %% we're going to solve
51 %% $$ -\delta^2 \Lap w_{k+1} + B'(w_k) \delta w_{k+1} = 2 * w_k$$
53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55 %% set $$ w_1 = win $$
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));
63 L = d2*QDDGOXCOMPDENS_LAP;
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));
69 G = fermiin - vin + 2*log(w);
70 Bmat = QDDGOXCOMPDENS_MASS*sparse(diag(G));
72 %%%%%%%%%%%%%%%%%%%%%%%%
73 %%% NEWTON ITERATION START
74 %%%%%%%%%%%%%%%%%%%%%%%%
76 for jnewt =1:ceil(maxit/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);
86 w = Urrextrapolation(w);
88 %%%%%%%%%%%%%%%%%%%%%%%%
89 %%% NEWTON ITERATION END
90 %%%%%%%%%%%%%%%%%%%%%%%%
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);
98 global QDDGOXCOMPDENS_LAP QDDGOXCOMPDENS_MASS QDDGOXCOMPDENS_RHS
104 res0 = norm((L + Bmat) * w,inf);
107 %% chose $ t_k $ to ensure positivity of $w$
117 tmpmat = QDDGOXCOMPDENS_MASS*2;
118 if (isempty(QDDGOXCOMPDENS_RHS))
119 QDDGOXCOMPDENS_RHS = Ucompconst (mesh,ones(Nnodes,1),ones(Nelements,1));
121 tmpvect= 2*QDDGOXCOMPDENS_RHS.*w;
123 %%%%%%%%%%%%%%%%%%%%%%%%
124 %%% DAMPING ITERATION START
125 %%%%%%%%%%%%%%%%%%%%%%%%
127 %% Compute $ B1mat = \frac{2}{t_k} $
128 %% and the (lumped) mass matrix B1mat(w_k)
131 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
132 A = L + B1mat + Bmat;
135 %% Apply boundary conditions
138 b = b - A (:,Dnodes) * Dvals;
146 wnew(Varnodes) = A\b;
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));
153 res = norm((L + Bmat) * wnew,inf);
161 %%%%%%%%%%%%%%%%%%%%%%%%
162 %%% DAMPING ITERATION END
163 %%%%%%%%%%%%%%%%%%%%%%%%