% pde_ex1.m
% This program demonstrates the use of the
% command line interface of the pde toolkit
% for solving BVPs in 2-D.
function iflag = pde_ex1();
iflag = 0;

show_all_figs = 0;

% using geometry function polygon1_geom below,
% show geometry of the computational domain
figure;
pdegplot('polygon1_geom');

% initialize mesh
[P,E,T] = initmesh('polygon1_geom');
% jiggle mesh
P = jigglemesh(P,E,T);
% plot mesh
if(show_all_figs)
    figure;
end
pdemesh(P,E,T);
% do global refinement
[P2,E2,T2] = refinemesh('polygon1_geom',P,E,T);
if(show_all_figs)
    figure;
end
pdemesh(P2,E2,T2);

% demonstrate Delaunay tesselation directly
% for points in initial mesh
x = P(1,:)';  y = P(2,:)';
tri = delaunay(x,y);
if(show_all_figs)
    figure;
end
triplot(tri,x,y);

% form voronoi polyhedra
[vx,vy] = voronoi(x,y,tri);
if(show_all_figs)
    figure;
end
voronoi(x,y,tri);

% retain fine mesh and discard initial one
P = P2; E = E2;  T = T2;
clear P2 E2 T2;

% compute values of f = x^2 + y^2 at node points
x = P(1,:)';  y = P(2,:)';
f = zeros(size(x));
f = (x.^2) + (y.^2);

% plot f(x,y) in various formats
% 3-D surface plot showing mesh
if(show_all_figs)
    figure;
end
pdeplot(P,E,T,'zdata',f);
% produce smoother plot by interpolating to
% standard (x,y) positions and using regular
% plotting routine
if(show_all_figs)
    figure;
end
pdeplot(P,E,T,'xydata',f,'zdata',f);
% switch to more common color scheme and add title
if(show_all_figs)
    figure;
end
pdeplot(P,E,T,'xydata',f,'zdata',f,...
    'colormap','jet','title','f = x^2 + y^2');
xlabel('x');  ylabel('y');
% make 2-D filled contour plot
if(show_all_figs)
    figure;
end
pdeplot(P,E,T,'xydata',f,'zdata',f,...
    'colormap','jet','title','f = x^2 + y^2', ...
    'zstyle','off');
xlabel('x');  ylabel('y');


% Now, we solve Poisson equation with source term f
% using the adaptive mesh solver for elliptic PDEs.
% type doc adaptmesh to see syntax, and doc assempde
% for information on how we specify coefficients of PDE
geom_mfile = 'polygon1_geom';
bound_mfile = 'pde_ex1_bound';
c = 1;  a = 0;  f = 'pde_ex1_f';
[u,P2,E2,T2] = adaptmesh(...
    'polygon1_geom', 'pde_ex1_bound', ... % names of m-files
    c, a, f);  % coefficients in PDE
% plot the solution
figure;
pdeplot(P2,E2,T2,'xydata',u,'zdata',u,...
    'colormap','jet','title','u(x,y)', ...
    'zstyle','off');
xlabel('x');  ylabel('y');

% Finally, make a master plot
figure;
% mesh returned from adaptmesh
subplot(2,2,1);  pdemesh(P2,E2,T2);
title('Mesh with adaptive refinement during solution');
axis([0 3 0 3]);
% source function
subplot(2,2,2);
x = P(1,:)';  y = P(2,:)';
f = (x.^2) + (y.^2);
pdeplot(P,E,T,'xydata',f,'zdata',f,...
    'colormap','jet','title','f = x^2 + y^2', ...
    'zstyle','off');
axis([0 3 0 3]);
% solution, filled contour plot
subplot(2,2,3);
pdeplot(P2,E2,T2,'xydata',u,'zdata',u,...
    'colormap','jet','title','u(x,y)', ...
    'zstyle','off');
axis([0 3 0 3]);
% compute gradient of solution
[ux,uy] = pdegrad(P2,T2,u);
subplot(2,2,4);
pdeplot(P2,E2,T2,'xydata',u,'zdata',u, ...
    'colormap','jet', 'zstyle','off', ...
    'flowdata', [ux; uy]);
axis tight;

% save results to .mat file
save pde_ex1.mat;

iflag = 1;
return;