   function [A, XY] = hdisk2(nx)
%%  function [A, XY] = hdisk2(nx)
%% generates a half disk domain with 
%% a corresponding FEM matrix. 
%% calls  assmbl1
%%---------- virtual x-y grid
if (nargin == 0) 
  nx     = 21;  
end 
hx = 1/(nx-1);
theta  = [0:hx:1]*pi;
%%---------- generate points 
r = [0.4:0.2:1];
nr = length(r);
X   = kron(r,cos(theta)) ;
Y   = kron(r,sin(theta));
%%---------- to count # triangles 
nt=0;
%%---------- generate triangles 
for j=1:nr-1
    for i=1:nx-1
        nt = nt+1;
        l = (j-1)*nx+ i;
%%---------- first triangle in box: 
        tri(nt,:) = [l, l+1, l+nx+1] ;
%%---------- Second triangle in box: 
        nt = nt+1;
        tri(nt,:) = [l,l+nx+1,l+nx];
    end
end
XY  = [X', Y'];
A = assmbl1(tri, XY); 
figure(1)
gplot(A,XY,'b-');
hold on
axis([-1.1 1.1 0.0 1.1]);
axis equal;
plot([-1 1], [0 0], 'r-', 'linewidth',2);
plot([0 0], [0 1 ], 'r-', 'linewidth',2);

