meshada

 adapt mesh for pde in struct. p, called by meshadac, 
 no interpol of tau here since done in meshadac

0001 function p=meshada(p,varargin) 
0002 % adapt mesh for pde in struct. p, called by meshadac,
0003 % no interpol of tau here since done in meshadac
0004 maxt=p.maxt; ngen=p.ngen; errbound=p.errbound; noa=nargin-1; k=1; 
0005 while k<=noa 
0006    switch lower(varargin{k})
0007        case 'maxt'; maxt=varargin{k+1}; k=k+2;
0008        case 'ngen'; ngen=varargin{k+1}; k=k+2;
0009        case 'eb'; errbound=varargin{k+1}; k=k+2;
0010        otherwise; break
0011    end
0012 end
0013 g=p.geo; t=p.tria; po=p.points; e=p.edges; alfa=0.15;beta=0.15;mexp=1;Par=0.5; 
0014 Tripick='pdeadworst';Rmethod='longest'; % Default values
0015 np=p.np;nt=size(t,2);gen=0; u=p.u; 
0016 while 1 
0017   [c,a,f,b]=p.f(p,p.u,p.lam); if any(b) f=bgradu2f(p,f,b,p.u); end
0018   errv=pdejmps(p.points,p.tria,c,a,f,p.u,alfa,beta,mexp);
0019   p.err=max(max(errv));fprintf(' %g\n', p.err); 
0020   if(p.err<errbound/2); fprintf('\nerr-est<p.errbound/2\n'); break; end; 
0021   if gen>=ngen,fprintf('\nnumber of refinements >= max\n'); break; end
0022   if nt>maxt, fprintf('number of triangles >= max\n'); break; end 
0023   i=feval(Tripick,po,t,c,a,f,u,errv,Par); % size(i), pause
0024   if isempty(i); fprintf('No triangles for refinement found!\n'); break; end 
0025   tl=i'; if size(tl,1)==1; tl=[tl;tl]; end
0026   u=reshape(u,np,length(u)/np); 
0027   [po,e,t,u]=refinemesh(g,po,e,t,u,tl,Rmethod); np=size(po,2); nt=size(t,2);
0028   u=u(:);p.points=po;p.edges=e;p.tria=t;p.np=np;p.nt=nt; gen=gen+1;
0029   fprintf('number of triangles=%g, ',nt); ims=p.imax; p.imax=6*ims; 
0030   [p.u,p.res,p.iter]=nlooppde(p,u,p.lam); p.imax=ims; 
0031   fprintf('res=%g,  error-est=',p.res); 
0032 end