function initialConditionsTest2D()
% initialConditionsTest2D: test initial condition routines in 2 dimensions.
%
%   initialConditionsTest2D (no arguments)
%
% This function (basically a script file) generates a sequence of 
%   shapes built by constructive solid geometry methods.
%
% We show both the 2D implicit surface (a contour) and the implicit
%   surface function (a 3D surface plot).

% Copyright 2004 Ian M. Mitchell (mitchell@cs.ubc.ca).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Ian Mitchell, 6/23/04

run('../addPathToKernel');

%---------------------------------------------------------------------------
% Standard output file id.
fid = 1;

%---------------------------------------------------------------------------
fprintf(fid, '\n%s %s\n%s', ...
             'A sequence of implicit surface functions in 2D', ...
             'along with their descriptions.', ...
             'Press any key after each display to move to the next');

%---------------------------------------------------------------------------
% A 2D grid.
g2.dim = 2;
g2.min = -2;
g2.max = +2;
g2.N = 101;
g2 = processGrid(g2);

%---------------------------------------------------------------------------
% Some basic shapes.
sphere1 = shapeSphere(g2);
show2D(g2, sphere1);
fprintf(fid, '\nUnit circle at the origin, generated by shapeSphere');
pause;

sphere2 = shapeSphere(g2, [ -0.5; +0.5 ]);
show2D(g2, sphere2);
fprintf(fid, '\nUnit circle in upper left quadrant, generated by shapeSphere');
pause;

cylinder1 = shapeCylinder(g2, 2, [ 0.5; 0.5 ], 0.5);
show2D(g2, cylinder1);
fprintf(fid, '\nVertical slab width 0.5 centered to right at 0.5, %s', ...
             'generated by shapeCylinder');
pause;

rectangle1 = shapeRectangleByCorners(g2);
show2D(g2, rectangle1);
fprintf(fid, '\nUnit square in upper right quadrant, %s', ...
             'generated by shapeRectangleByCorners');
pause;

rectangle2 = shapeRectangleByCorners(g2, [ -1.0; -Inf ], [ +Inf, -1.0 ]);
show2D(g2, rectangle2);
fprintf(fid, '\nUnbounded down and left of [ -1, -1 ], %s', ...
             'generated by shapeRectangleByCorners');
pause;

rectangle3 = shapeRectangleByCenter(g2, [ 1.0; -0.5 ], [ 0.5; 1.0 ]);
show2D(g2, rectangle3);
fprintf(fid, '\nCentered at [ 1, -0.5 ], width [ 0.5, 1.0 ], %s', ...
             'generated by shapeRectangleByCenter');
pause;

hyperplane1 = shapeHyperplane(g2, [ 3.0; 4.0 ], [ 1; 1 ]);
show2D(g2, hyperplane1);
fprintf(fid, '\nNormal [ 3, 4 ] through [ 1, 1 ], %s', ...
             'generated by shapeHyperplane');
pause;

hyperplane2 = shapeHyperplane(g2, [ -1.0; +1.0 ]);
show2D(g2, hyperplane2);
fprintf(fid, '\nNormal [ -1, +1 ] through origin, %s', ...
             'generated by shapeHyperplane');
pause;

%---------------------------------------------------------------------------
% Constructive geometry operations.
union1 = shapeUnion(sphere1, sphere2);
show2D(g2, union1);
fprintf(fid, '\nUnion of first two spheres, %s', ...
             'generated by shapeUnion (NOTE: not signed distance)');
pause;

intersection1 = shapeIntersection(cylinder1, hyperplane1);
show2D(g2, intersection1);
fprintf(fid, '\nIntersection of slab and first hyperplane, %s', ...
             'generated by shapeIntersection');
pause;

complement1 = shapeComplement(rectangle2);
show2D(g2, complement1);
fprintf(fid, '\nEverything but unbounded down and left of [ -1, -1 ], %s', ...
             'generated by shapeComplement');
pause;

difference1 = shapeDifference(sphere1, rectangle1);
show2D(g2, difference1);
fprintf(fid, '\nFirst sphere subtract first unit rectangle, %s', ...
             'generated by shapeDifference');

%---------------------------------------------------------------------------
% Now let's build a pac-man shape (non-convex 2D polygon)
%   by intersections and unions of hyperplanes.

% Build the top first, which is a convex polygon.
ns = [ -1, 0; -1, +1; +1, +2; +1, -2; 0, -1 ];
ps = [ -1, 0; -1, 0.5; 0, +1; 0, 0; 0, -0.5 ];
pacTop = -inf * ones(g2.shape);
for i = 1 : size(ns, 1)
  pacTop = shapeIntersection(pacTop, shapeHyperplane(g2, ns(i,:)', ps(i,:)'));
end

% Bottom is just the top flipped vertically.
pacBottom = flipdim(pacTop, 2);

% Combine them with union to form pac-man shape.
pac = shapeUnion(pacTop, pacBottom);

show2D(g2, pac);
fprintf(fid, '\npac-man shape, generated by combining %s', ...
             'shapeHyperplane, shapeIntersection and shapeUnion');


%---------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------------------------------------
function h = show2D(g, data)
% show2D: display a 2D implicit surface function.
%
%   show2D(g, data)
%
% 2D implicit surface functions can be displayed either by contour plot
%   or by surface plot.  This routine does both.
%
% parameters:
%   g   	Grid structure (see processGrid.m for details).
%   data        Array containing the implicit surface function.


%---------------------------------------------------------------------------
% Set up two figures in which to place the results.
%   Use the same figures every time.
persistent f1 f2
if(isempty(f1))
  f1 = figure;
end
if(isempty(f2))
  f2 = figure;
end

%---------------------------------------------------------------------------
% Contour plot of implicitly defined set.
figure(f1);
contour(g.xs{1}, g.xs{2}, data, [ 0 0 ], 'b-');
grid on;
axis equal;
axis(g.axis);
xlabel('x');  ylabel('y');

%---------------------------------------------------------------------------
% Surface plot of implicit surface function.
figure(f2);
surf(g.xs{1}, g.xs{2}, data);
xlabel('x');  ylabel('y');  zlabel('\phi(x,y)');
grid on;
