meshref

 adapt (i.e., refine locally) mesh

0001 function p=meshref(p,varargin) 
0002 % adapt (i.e., refine locally) mesh
0003 maxt=p.maxt; ngen=p.ngen; errbound=p.errbound; noa=nargin-1;k=1; 
0004 while k<=noa; 
0005     switch lower(varargin{k})
0006     case 'maxt'; maxt=varargin{k+1}; k=k+2;
0007     case 'ngen'; ngen=varargin{k+1}; k=k+2;  
0008     case 'eb'; errbound=varargin{k+1}; k=k+2;
0009     otherwise; break; 
0010     end
0011 end
0012 fprintf('current nt=%g, aiming at nt=%g and err=%g, starting error= ', ...
0013     p.nt,maxt,errbound);
0014 g=p.geo; t=p.tria; po=p.points; e=p.edges; 
0015 q=p; lamd=p.tau(p.neq*p.np+1); ps=po; nps=p.np; nt=size(t,2);
0016 alfa=0.15;beta=0.15;mexp=1;Par=0.5;  % Default values
0017 Tripick='pdeadworst';Rmethod='longest'; gen=0; u=p.u; np=p.np; 
0018 while 1 %(size(t,2)<maxt && mpdejmps>errbound && gen<=p.ngen)
0019   [c,a,f,b]=p.f(p,p.u,p.lam); if any(b) f=bgradu2f(p,f,b,p.u); end
0020   errv=pdejmps(p.points,p.tria,c,a,f,p.u,alfa,beta,mexp);
0021   merr=max(max(errv));fprintf(' %g\n', merr); 
0022   if(merr<errbound); fprintf('\nerr-est<p.errbound\n'); break; end; 
0023   if gen>ngen,fprintf('\nnumber of refinements > max\n'); break; end
0024   if nt>maxt, fprintf('number of triangles > max\n'); break; end 
0025   i=feval(Tripick,po,t,c,a,f,u,errv,Par); % size(i), pause
0026   if isempty(i); fprintf('No triangles for refinement found!\n'); break; end 
0027   tl=i'; if size(tl,1)==1; tl=[tl;tl]; end
0028   u=reshape(u,np,length(u)/np); 
0029   [po,e,t,u]=refinemesh(g,po,e,t,u,tl,Rmethod); np=size(po,2); nt=size(t,2);
0030   u=u(:);p.points=po;p.edges=e;p.tria=t;p.np=np;p.nt=nt; gen=gen+1;
0031   fprintf('number of triangles=%g, ',nt); ims=p.imax; p.imax=6*ims; 
0032   [p.u,p.res,p.iter]=nlooppde(p,u,p.lam); p.imax=ims; 
0033   fprintf('res=%g, error-est=',p.res);   
0034 end
0035 if(p.vsw>0) 
0036  figure(p.ifig);pdemesh(p.points,p.edges,p.tria);axis tight; 
0037  plotsol(p,p.pfig,p.pcmp,p.pstyle); 
0038 end
0039 xi=1/p.np; % since refinement, rather decrease xi
0040 if (p.isw>1) fprintf('old xi=%g, new xi=%g\n',p.xi,xi); p.xi=asknu('xi:',xi); 
0041 else p.xi=xi; % if no interaction then use new standard xi
0042 end 
0043 tau=zeros(p.np*p.neq,1); x=ps(1,:); y=ps(2,:);
0044 xn=p.points(1,:); yn=p.points(2,:);
0045 for i=1:p.neq % interpolate tau
0046     z=p.tau((i-1)*nps+1:i*nps);taun=p2interpol(xn,yn,z,x,y); 
0047     tau((i-1)*p.np+1:i*p.np)=taun; 
0048 end
0049 tau(p.np*p.neq+1)=lamd;  % append ds coordinate to tau and normalize
0050  p.tau=tau/xinorm(tau,p.xi); p.headfu(p); 
0051 
0052