function [x,y]=geome(bs,s)
%
%Geometry M-file for a particular problem. This Geometry M-file can be used
%by mesh generators within the matlab contex. 
%

%Number of boundary segments for this particular geometry is 7:
nbs=7;

%If the function is called with no arguments, return only the number
%of boundary edges:
if nargin == 0
  x=nbs;
  return
end

%Parameters used to create the boundary geometry:
d1=1;
d2=d1/3;
d3=1;
d4=d3/2;
r=d1;

%The parametrization parameter s will run along the boundary 
%(i.e. s is the arc-length). Specify the values of s at every edge end points 
%(edge e(i) has end points parametrized by s(i)=si and s(i+1).
s_1=d1;  
s_2=s_1+1/3/r*d2^2*((1+2*r/d2)^(3/2)-1);
s_3=s_2+d1+d2;
s_4=s_3+pi/2*r;
s_5=s_4+d3+d4;
s_6=s_5+1/3/r*d4^2*((1+2*r/d4)^(3/2)-1);
s_7=s_6+d3;

%The total length of the boundary is s7 or:
L =2*d1+2*d3+d2+d4+pi/2*r+...
      1/3/r*(d2^2*((1+2*r/d2)^(3/2)-1)+d4^2*((1+2*r/d4)^(3/2)-1));

%The following matrix describes the boundary:
%1st row: values of the parametrization parameter where each edge starts.
%2nd row: values of the parametrization parameter where each edge ends.
%3rd row: Label for left hand region.
%4rd row: Label for right hand region. (the region outside the problem region
%                                       has label 0)
%Each column describes a boundary segment.
dl=[ 0   s_1  s_2  s_3  s_4  s_5  s_6;
    s_1  s_2  s_3  s_4  s_5  s_6  s_7;
     1    1    1    1    1    1    1 ;
     0    0    0    0    0    0    0 ];

%If the function is called using one argument (i.e. bs), return the columns
%of matrix dl especified by the elements of bs. The elements of bs 
%must be positive integers. Otherwise, matlab rounds to the next bigger 
%integer.
if nargin == 1
  x=dl(:,bs);
  return
end

%If the function call uses 2 arguments, we must check that bs has the 
%same size than s, or bs is a scalar. 
x=zeros(size(s));
y=zeros(size(s));
[m,n]=size(bs);
if m==1 & n==1,
  bs=bs*ones(size(s)); %expand bs to the size of s
elseif m~=size(s,1) | n~=size(s,2),
  error('bs must be scalar or of the same size as s')
end

%The following defines the boundary geometry using 7 'pieces' and then 
%assembling all of them together. The parameter is ss.
ss_1=d1;  
ss_2=ss_1+d2;
ss_3=ss_2+d1+d2;
ss_4=ss_3+pi/2;
ss_5=ss_4+d3+d4;
ss_6=ss_5+d4;
ss_7=ss_6+d3;

ds=1/400;
ss1=[0:ds:ss_1];
ss2=[ss_1+ds:ds:ss_2];
ss3=[ss_2+ds:ds:ss_3];
ss4=[ss_3+ds:ds:ss_4];
ss5=[ss_4+ds:ds:ss_5];
ss6=[ss_5+ds:ds:ss_6];
ss7=[ss_6+ds:ds:ss_7];
ss=[ss1 ss2 ss3 ss4 ss5 ss6 ss7];

xx1=ss1;
yy1=zeros(size(ss1));
xx2=ss2;
yy2=r/d2*(ss2-ss_1).*(2-(ss2-ss_1)/d2);
xx3=ss_3-ss3;
yy3=r*ones(size(ss3));
xx4=-r*sin(ss4-ss_3);
yy4=r*cos(ss4-ss_3);
xx5=-r*ones(size(ss5));
yy5=-(ss5-ss_4);
xx6=r/d4*(ss6-ss_5).*(2-(ss6-ss_5)/d4)-r;
yy6=d4/(ss_6-ss_5)*(ss6-ss_5)-(d3+d4);
xx7=zeros(size(ss7));
yy7=ss7-ss_6-d3;
xx=[xx1 xx2 xx3 xx4 xx5 xx6 xx7];
yy=[yy1 yy2 yy3 yy4 yy5 yy6 yy7];

%Restrict the parameter ss to the values of parameter s (especified as input).
%The parameter s takes (a) value/s from 0 to L.
ss=pdearcl(ss,[xx;yy],s,0,L);

for i=1:length(ss)
if ss(i) >= 0 & ss(i) <= ss_1
  x(i)=ss(i);
  y(i)=0;
elseif ss(i) >= ss_1 & ss(i) <= ss_2
  x(i)=ss(i);
  y(i)=r/d2*(ss(i)-ss_1).*(2-(ss(i)-ss_1)/d2);
elseif ss(i) >= ss_2 & ss(i) <= ss_3
  x(i)=ss_3-ss(i);
  y(i)=r;
elseif ss(i) >= ss_3 & ss(i) <= ss_4
  x(i)=-r*sin(ss(i)-ss_3);
  y(i)=r*cos(ss(i)-ss_3);
elseif ss(i) >= ss_4 & ss(i) <= ss_5
  x(i)=-r;
  y(i)=-(ss(i)-ss_4);
elseif ss(i) >= ss_5 & ss(i) <= ss_6
  x(i)=r/d4*(ss(i)-ss_5).*(2-(ss(i)-ss_5)/d4)-r;
  y(i)=d4/(ss_6-ss_5)*(ss(i)-ss_5)-(d3+d4);
elseif ss(i) >= ss_6 & ss(i) <= ss_7
  x(i)=0;
  y(i)=ss(i)-ss_6-d3;
end
end

%figure, plot(xx,yy), axis('equal'), grid on
%hold on, plot(x,y), axis('equal'), grid on

