meshcheck

 check if FEM-solution is converged to PDE-solution by refining mesh 
 and comparing refined solution u1=p1.u with u0=p.u; 
 for this calc difference ud=u1-un with un=interpol of u0 to new mesh 
 Print max(abs(ud)). If sw>0 also plot.

0001 function [p1,ud]=meshcheck(p,varargin)
0002 % check if FEM-solution is converged to PDE-solution by refining mesh
0003 % and comparing refined solution u1=p1.u with u0=p.u;
0004 % for this calc difference ud=u1-un with un=interpol of u0 to new mesh
0005 % Print max(abs(ud)). If sw>0 also plot.
0006 err=errcheck(p); fprintf('error-estimate: %g\n', err); 
0007 if (nargin-1==0) sw=0; else sw=varargin{1}; end 
0008 p0=p;p1=p;neq=p.neq; p1.errbound=0;
0009 p1.isw=0; % switch off interaction during meshref
0010 p1=meshref(p1,'maxt',2*p1.nt);
0011 un=zeros(neq*p1.np,1); % now interpolate p0.u to new grid
0012 x=p0.points(1,:); y=p0.points(2,:); 
0013 xn=p1.points(1,:);yn=p1.points(2,:); onp=size(p1.points,2);
0014 for i=1:neq 
0015    z=p0.u((i-1)*p0.np+1:i*p0.np);
0016    ui=p2interpol(xn,yn,z,x,y); 
0017    un((i-1)*onp+1:i*onp)=ui; 
0018 end
0019 ud=p1.u-un; m1=max(abs(p1.u)); m2=max(abs(ud)); m3=m2/m1;
0020 fprintf('infty-norm(un-uo)=%g,  rel-err=%g\n',m2,m3);
0021 p0=p1; p0.u=ud; % put difference into p0 for plotting
0022 if(sw>0) plotsol(p0,p0.ifig,sw,p0.pstyle); end