%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to CMSC/AMSC 661 Spring 2005 Homework 2      
%  hmwk2.m Dianne P. O'Leary   03/2005                   
%
%  In this homework, we use Matlab's pdetool to solve
%  6 problems:
%          - grad . (a grad u) = f in Omega
%  with Dirichlet boundary conditions and true solution
%          u(x,y) = -(x-.5).^3+cos(x.*y)
%  We use three domains:
%     -- the unit circle
%     -- a 'pacman' shape
%     -- a convex polygon.
%  We use 2 coefficient definitions:
%     --  a = 1 (Poisson's equation)
%     --  a discontinuous at x = -.5
%  We solve each problem over 3 meshes (2 refinements)
%  and measure the infinity norm of the error at the
%  mesh points, and the ratio of errors over successive
%  meshes.
%  This code is derived from that produced by Matlab's
%  GUI pdetool.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over the domains, and initialize the pdetool.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for idomain=1:3,
  [pde_fig,ax]=pdeinit;
  pdetool('appl_cb',1);
  set(ax,'DataAspectRatio',[1 1 1]);
  set(ax,'PlotBoxAspectRatio',[1.5 1 1]);
  set(ax,'XLim',[-1.5 1.5]);
  set(ax,'YLim',[-1 1]);
  set(ax,'XTickMode','auto');
  set(ax,'YTickMode','auto');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the geometry description for the three domains.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if (idomain==1)
% Geometry description:
    sgeom = 'Circle';
    pdeellip(0,0,1,1,0,'E1');
%   pdecirc(0,0,1,'E1');  should work as well, but error bigger
    set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','E1')
    nseg = 4;
    area = pi;
  elseif (idomain==2)
    sgeom = 'Pacman';
    pdeellip(0,0,1,1,0,'E1');
    pdepoly([0,5,5*cos(pi/3)],[0,0,5*sin(pi/3)],'P1');
    set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','E1-P1')
    nseg = 6;
    area = pi*240/360;
  elseif(idomain==3)
    sgeom = 'polygon';
    pdepoly([1,1,-.8,-1,.2],[0,1,.5,-.5,-.4],'E1');
    set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','E1')
    nseg = 5;
    area = 2; % a gross approximation
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the boundary conditions.  
% The domain has nseg segments.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  pdetool('changemode',0)
  for i=nseg:-1:1,
    pdesetbd(i,'dir',1,'1','-(x-.5).^3+cos(x.*y)')
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over the PDE coefficients and set  a.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for icoef = 1:2,

    if (icoef == 1)
      scoef = 'the Poisson equation';
      pdeseteq(1,'1.0','0.0','(6*(x-.5) + cos(x.*y).*(y.^2+x.^2))',...
              '1.0','0:10','0.0','0.0','[0 100]')
      setuprop(pde_fig,'currparam',...
          ['1.0                                ';...
          '0.0                                ';...
          '(6*(x-.5) + cos(x.*y).*(y.^2+x.^2))';...
          '1.0                                '])
     else
      scoef = 'with discontinuous coefficients';
      pdeseteq(1,'(x<=(-.5)).*(x+5)+4.5*(x>(-.5))','0.0',...
      '(x<=-.5).*((x+5).*(cos(x.*y).*(x.^2+y.^2)+6*(x-.5))+3*(x-.5).^2+y.*sin(x.*y))+4.5*(x>-.5).*(6*(x-.5)+cos(x.*y).*(y.^2+x.^2))',...
      '1.0','0:10','0.0','0.0','[0 100]')
      setuprop(pde_fig,'currparam',...
      ['(x<=(-.5)).*(x+5)+4.5*(x>(-.5))                                                                                             ';...
       '0.0                                                                                                                         ';...
       '(x<=-.5).*((x+5).*(cos(x.*y).*(x.^2+y.^2)+6*(x-.5))+3*(x-.5).^2+y.*sin(x.*y))+4.5*(x>-.5).*(6*(x-.5)+cos(x.*y).*(y.^2+x.^2))';...
       '1.0                                                                                                                         '])
     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate the initial mesh. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    setuprop(pde_fig,'Hgrad',1.3);
    setuprop(pde_fig,'refinemethod','regular');
    pdetool('initmesh')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the solver and plotting parameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    setuprop(pde_fig,'solveparam',...
        str2mat('0','1000','10','pdeadworst',...
        '0.5','longest','0','1E-4','','fixed','Inf'))

% Plotflags and user data strings:
    setuprop(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]);
    setuprop(pde_fig,'colstring','');
    setuprop(pde_fig,'arrowstring','');
    setuprop(pde_fig,'deformstring','');
    setuprop(pde_fig,'heightstring','');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over the meshes.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ngrids = 5;
    for i=1:ngrids

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate the initial mesh (if i=1) or refine the mesh.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       if (i>1)
         pdetool('refine')
       end
       pdetool('solve')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Export information from PDETool to workspace and 
% compute the error at the mesh points.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

       [pmesh,e,t,ucomp,l,c,a,f,d,b,g] = getpetuc;
       utrue =-(pmesh(1,:)-.5).^3 + cos(pmesh(1,:).*pmesh(2,:));
       errpde(i) = norm(ucomp-utrue','inf');
       ntri(i) =max(size(t));
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Print the results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    disp(sprintf('Solving %s on the %s domain.',scoef,sgeom))
    disp(sprintf('Errors on the the five meshes:'))
    disp(sprintf(' %8.3e ',errpde))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We estimate the mesh size h by assuming that we have
% divided the domain into equilateral right triangles.
%
% We calculate pred_h, the value of h that matches the
% error ratio, assuming that convergence is O(h^2 log(1/h).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hest  = sqrt(2)*sqrt(area./ntri);

% Ratio of errors; asymptotically, should be 4
est = [0, errpde(1:ngrids-1)./errpde(2:ngrids)];
pred_h = [0,exp(-est(2:ngrids)*log(2)./(4-est(2:ngrids)))];

disp('grid  estimated_h predicted_h  error ratio   ')
for i=1:ngrids,
  disp(sprintf('%3d      %7.2e   %7.2e   %7.2e   ', ...
           i,hest (i),pred_h(i),est(i)))
end
  end % for icoef

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare for the next problem by allowing the pdetool
% to forget the current data. (Lamphier/Harvey/Mathworks)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  pde_fig=findobj(allchild(0),'flat','Tag','PDETool');
  if ~isempty(pde_fig),
    h = findobj(get(pde_fig,'Children'),'flat','Tag','PDEFileMenu');
    flags = get(h,'UserData');
    flags(1) = 0;
    set(h,'UserData',flags)
  end

end % for idomain