function [sgrid] = generate_ijk_grid(flag,nsubx,nsuby,lx,ly)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% M. Peszynska, 02/2005
%% © Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function [grid] = generate_ijk_grid(flag,nx,ny,lx,ly)
%% generates rectangular grid representing the shape <flag> 
%% with size parameters nx,ny,lx,ly 
%% grid is padded by zeros which denote boundary nodes
%% or those outside the region
%% nsubx,nsuby are numbers of intervals in each direction
%% resulting grid is node based: nodes are 1..nsubx+1
%% internal nodes are 2..nsubx 
%% grid(i,j) = 1 mena sthe node is in the grid
%% numbering of nodes can be obtained with myindex.m
%%
%% on some simple regions generate_ijk_grid followed by mynumindex is
%% equivalent to MATLAB's numgrid
%% the grids can be later merged using mergegrids.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% only few shapes are implemented 
%% flag = 'R' generates a rectangle
%% flag = 'F' generates a rectangle without the center
%% flag = 'L' generates an 'L' shape within nx x ny rectangle, lx, ly
%% flag = 'C' generates a circular shape within nx x ny rectangle, lx, ly
%% flag = 'S' generates a 'step' shape within nx x ny rectangle, lx, ly
%% flag = 'Holes' generates a box with lx holes each of size 2*ly-1, centered in random locations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% examples:
%% fiss = generate_ijk_grid('F',10,10,2,2);
%% block = generate_ijk_grid('R',6,6)
%%%% (these generate an interesting compatible structure)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin <=2 , nsuby =nsubx; end;
nx = nsubx + 1; ny = nsuby + 1;

sgrid = zeros(nx,ny);
switch flag
    case 'R'
        for i=2:nx-1 for j=2:ny-1 sgrid(i,j)= 1;end;end;
    case 'F'
        if nargin <=3 lx = 1;ly = 1;end;
        for j=2: ly+1 for i=2:nx-1 sgrid(i,j)=1;end; end
        for j=ly+2: ny-ly -1 
            for i=2:lx+1 sgrid(i,j) = 1; end;
            for i=nx-lx:nx-1 sgrid(i,j) = 1; end;
        end;
        for j=ny-ly: ny-1 for i=2:nx-1 sgrid(i,j)=1; end; end
    case 'L'
        if nargin <= 3 lx=1;ly=1;end;
        for j=2:ly+1 for i=2:nx-1 sgrid(i,j) = 1;end; end
        for j=ly+2: ny-1 
            for i=2:lx+1 sgrid(i,j) = 1; end;
        end;
    case 'C'
        dx = 2/nx; dy = 2/ny; 
        for j=2:nx-1 for k=2:ny-1 
                if norm([-1+dx*j,-1+dy*k],2) < 1 sgrid(j,k)=1;end;
            end;end;
    case 'S'
        for j=ly+1:ny-1 for i=2:nx-1 sgrid(i,j) = 1; end; end;
        for i= lx+1: nx-1 for j=2:ly sgrid(i,j) =1; end; end;

    case 'Holes'
        %% generate box
        for i=2:nx-1 for j=2:ny-1 sgrid(i,j)= 1;end;end;
        %% add holes: lx  = number of inclusions , ly - random state from which we start
        if lx >= (nx-4) lx=nx-4; end; 
        rand('state',ly);
        ir = randperm(nx-4); jr = randperm(ny-4);
        for k = 1:lx
              ci = ir(k); cj = jr(k); lp = randperm(ly); l = ly;
              for m=-l:1:l  for n=-l:1:l 
                    sgrid(max(min(ci + m, nx-2),3),max(min(cj + n,ny-2),3) ) = 0; 
              end;end;
        end;
        sgrid(2,2:ny-1)= 1;
        sgrid(2:nx-1,2)= 1;
    case 'Multi'
        nx = nsubx * lx + 1; ny = nsuby *ly + 1;

        sgrid = zeros(nx,ny);
        for i=2:nx-1 for j=2:ny-1 sgrid(i,j)= 1;end;end;
        
end;
    


