% function [x,y] = polygon1_geom(bs,s);
% --------------------
% This routine is a geometry function for the
% polygon, and demonstrates the syntax described
% by "doc pdegeom". We decribe the geometry by
% defining the boundary that encloses it. We describe
% the boundary as a number of non-overlapping
% curves described by this function.
function [x,y] = polygon1_geom(bs,s);

% specify the number of boundary sections
% and store the nodal positions
nbs = 6;  % number of boundary sections
% set (x,y) positions of each vertex
V = [0, 0; ... % node 1 
     1, 1; ... % node 2
     0, 2; ... % node 3
     1.5, 3; ... % node 4
     3, 3; ... % node 5
     3, 0;];  % node 6

% Now, construct a matrix such that each row
% corresponds to a section of the boundary.
% The first column contains the identifier
% of start-point vertex, and the second
% contains the identifier of the end-point
% vertex.
B = zeros(nbs,2);
for k=1:nbs  % for each boundary section
    n1 = k; % vertex at start of section
    % get vertex at end of section
    if(k == nbs)
        n2 = 1;
    else
        n2 = k+1;
    end
    % store values in row k of B
    B(k,:) = [n1, n2];
end

% Now, store a matrix D that has six columns,
% one for each boundary section. We parameterize
% the boundary curve by a contour parameter 0 <= s <= 1,
% such that the section # k starts at vertex BS(k,1)
% at s = 0 and ends at vertex BS(k,2) at s = 1.
%
% In rows 3 and 4, we specify what type of regions are
% found to the left and right of our curve, as we move
% along the boundary in increasing s. Here, since we move
% clockwise around the boundary, the left-hand region is
% outside of the domain and we set a 0 in row 3. The
% right-hand region is inside the computational domain,
% and thus we set row 4 equal to 1.
%
% In this problem, there is only one region of the domain,
% but in other problems, we could specify that the domain
% is comprised of regions 1, 2, 3, ... and we could use
% different PDE coefficients in each region. We then would
% have additional "internal" boundary sections that divide
% the overall domain into these non-overlapping subdomains.
% We would record the appropriate subdomain identifier
% 0, 1, 2, 3, ... for the left-hand side in row 3 and for
% the right hand side in row 4.
D = zeros(4,nbs);
for k=1:nbs  % for each section in the boundary
    % set s = 0 in row 1
    D(1,k) = 0;
    % set s = 0 in row 2
    D(2,k) = 1;
    % as left-hand side of boundary section is
    % outside of the domain, set 0 in row 3
    D(3,k) = 0;
    % as right-hand side of boundary section is
    % inside of the domain, set 1 in row 4
    D(4,k) = 1;
end

% Now, we are ready to return the appropriate
% output arguments

% if there are no input arguments, return the number
% of segments along the boundary
if(nargin == 0)
    x = nbs;
    y = 0;  % dummy value
    return;
end

% if there is only one input argument, bs, return
% the submatrix DS containing the columns of D
% specified in bs
if(nargin == 1)
    x = D(:,bs);
    y = 0;  % dummy value
    return;
end

% Else, if there are two input arguments, bs contains
% some number of boundary sections and s the values
% of s along the section at which x(s) and y(s) are
% to be computed.

% compute size of bs
[bs_numRows,bs_numCols] = size(bs);
% if bs is a scalar, than we are to compute
% the values of x(s) and y(s) at each value
% of s along the same curve. Thus, we expand
% bs to be the same size as s.
if((bs_numRows == 1)&&(bs_numCols == 1))
    bs = bs*(ones(size(s)));
elseif((bs_numRows ~= size(s,1))|(bs_numCols ~= size(s,2)))
     error('pdegeom: size(bs) ~= size(s)');
end
% Now, for each of these points along the
% specified boundary, we compute the corresponding
% positions x(s) and y(s).  Here s parameterizes
% the curve proportional to the arc length, For
% boundary sections that are not line segments, this
% is not true, and the pdearc1 function should
% be used to correct for this.
x = zeros(size(s)); y = zeros(size(s)); 
for m=1:size(s,1)
    for n=1:size(s,2)
        % record label of boundary section
        k = bs(m,n);
        % find the start and end vertex labels
        n1 = B(k,1);  n2 = B(k,2);
        % find the start and end positions
        x1 = V(n1,1);  x2 = V(n2,1);
        y1 = V(n1,2);  y2 = V(n2,2);
        % use linear interpolation to find the
        % value of x(s) and y(s)
        x(m,n) = (1-s(m,n))*x1 + s(m,n)*x2;
        y(m,n) = (1-s(m,n))*y1 + s(m,n)*y2;
    end
end

return;