% 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;