cont

 main cont. routine for -div(c grad u)+au-f-b grad u=0, 
 c,a,f, b encoded in p.f, BC encoded in p.bcf 
 do p.nsteps steps of (initial) size p.ds, starting from p.u, p.lam, 
 direction tau, if it's already there

0001 function p=cont(p) 
0002 % main cont. routine for -div(c grad u)+au-f-b grad u=0,
0003 % c,a,f, b encoded in p.f, BC encoded in p.bcf
0004 % do p.nsteps steps of (initial) size p.ds, starting from p.u, p.lam,
0005 % direction tau, if it's already there
0006 if(p.pnamesw==1) % set filenames to current input variable name
0007     dum1=p.count; dum2=p.branch; dum3=p.bcount; % remeber counters
0008     [p,ok]=setfn(p,inputname(1)); 
0009     p.count=dum1; p.branch=dum2; p.bcount=dum3;
0010     if ok~=1; return; end
0011 end
0012 p.ntime=0; p.tstime=0; p.sptime=0; totime=cputime; % reset timing
0013 is=0; % stepcounter
0014 ineg0=-1; % to save last #neg EVals
0015 p.headfu(p);            % output user defined header
0016 while is<p.nsteps      % continuation loop (p.nsteps steps)
0017    if (length(p.tau)~=p.neq*p.np+1)||(p.restart>0) % initial step
0018       [p,iok]=inistep(p); if iok==0; return; end 
0019       plotsol(p,1,1,p.pstyle); is=is+2;
0020    end
0021    stime=cputime; iter=0; res=10*p.tol; stepok=0; % stepok=1 after successful step
0022    dss=p.ds; % save current ds for restart at current p.lam after mesh-ref-in cfail
0023    while stepok==0     % loop to find next cont point
0024       lamd=p.tau(p.neq*p.np+1); 
0025       u1=p.u+p.ds*p.tau(1:p.neq*p.np); lam1=p.lam+p.ds*lamd;  % predictor
0026       ntime=cputime;
0027       if(p.parasw==0 || (p.parasw==1 && abs(lamd)>p.lamdtol)) % fixed lam corrector
0028           [u1,res,iter,Gu,Glam]=nlooppde(p,u1,lam1); p.meth='nat'; 
0029       else % arclength-corrector
0030          [u1,lam1,res,iter,Gu,Glam]=nloopext(p,u1,lam1,p.ds); p.meth='arc'; 
0031       end
0032       p.ntime=p.ntime+cputime-ntime; % newton-loop time, accumulated
0033       ds=p.ds; % for output below, now stepsize control (re convergence)
0034       [p,stepok,u1,lam1,res,iter,Gu,Glam]=sscontrol(p,u1,lam1,res,iter,Gu,Glam,dss); 
0035       if(stepok==-1); return; end; % ABORT cont
0036    end              % stepok==0
0037    if p.spcalcsw>0 % calculate EVals
0038        sptime=cputime; ineg0=p.ineg; 
0039        [p.ineg,p.muv]=spcalc(Gu,p); 
0040        p.sptime=p.sptime+cputime-sptime; % spectral-time, accumulated
0041    end
0042    % form extended matrix and compute new tangent
0043    amat=[[Gu Glam]; [p.xi*p.tau(1:p.neq*p.np)' (1-p.xi)*p.tau(p.neq*p.np+1)]];  
0044    tau1=p.blss(amat,[zeros(p.neq*p.np,1);1],p,lam1); tau1=tau1/xinorm(tau1,p.xi); 
0045    if(p.bifchecksw>0)     % check for bifurcation via reduced |det(A)|
0046      newds=p.ds; p.ds=ds; % use ds from BEFORE last sscontrol!
0047      biftime=cputime; p=bifdetec(p,u1,lam1,tau1,Gu,Glam,ineg0); 
0048      p.biftime=p.biftime+cputime-biftime; p.ds=newds; 
0049    end
0050    p.restart=0; p.u=u1; p.lam=lam1; p.tau=tau1;      % step accepted!
0051    p.res=res; p.iter=iter; p.lamd=lamd; % store stuff (if e.g. p.ufu needs iter)
0052    if(p.errchecksw>0); p.err=errcheck(p);   % ERRCHECK and possibly mesh-refinement
0053      if(p.err>p.errbound && p.errbound>0)                            
0054        if(p.errchecksw==1) % just give warning!
0055          fprintf('err.est.=%g>errbound=%g. Consider mesh-refinement.\n', p.err, p.errbound);          
0056        end 
0057        if(p.errchecksw==2); % adapt mesh
0058            fprintf('err.est.=%g>errbound=%g. Adapting mesh\n', p.err, p.errbound); 
0059            p=meshadac(p,'eb',p.errbound); 
0060        end 
0061        if(p.errchecksw>2); % refine mesh
0062           p=meshref(p,'eb',p.errbound); 
0063        end          
0064      end %p.err>p.errbound
0065    end %p.errchecksw>0
0066    if mod(p.count,p.amod)==0; % adapt mesh ignoring errbound
0067        fprintf('adapting mesh\n'); p=meshadac(p,'maxt',p.maxt,'eb',0); end 
0068    if mod(p.count,p.brmod)==0 % put step branch
0069        brout=p.outfu(p,p.u,p.lam); % userfu to append to bif-branches
0070        p.branch=[p.branch [p.count; p.ineg; p.lam; p.err; brout]]; % put on branch
0071        figure(p.brfig); hold on; % plot point
0072        if p.ineg<=0; plot(lam1,real(brout(p.bpcmp)),'*'); 
0073        else plot(lam1,real(brout(p.bpcmp)),'+'); end 
0074    end 
0075    if(p.count>0 && mod(p.count,p.smod)==0) % save to file
0076        fname=[p.pname,sprintf('%i',p.count),'.mat']; save(fname,'p'); end
0077    p.tstime=p.tstime+cputime-stime; % total step time (accumulated)
0078    if(mod(p.count,p.pmod)==0); plotsol(p,p.pfig,p.pcmp,p.pstyle); end % plot sol
0079    cstop=p.ufu(p,brout,ds); % user function, typically printout
0080    p.count=p.count+1; is=is+1; 
0081    if(cstop==1) break; end % p.ufu returned stop!
0082 end % while is<p.nsteps
0083 % some postprocessing, i.e., save the last point
0084 if(p.smod~=0); fname=[p.pname,sprintf('%i',p.count-1),'.mat']; save(fname,'p'); end 
0085 p.totime=cputime-totime; 
0086 if(p.timesw>0) 
0087   fprintf('Timing: total=%g, av.step=%g, av.newton=%g, av.spcalc=%g\n',...
0088       p.totime,p.tstime/is,p.ntime/is,p.sptime/is);
0089 end