function elements=SElem(type, gp, NXQ)

%
% elements=SElem(type, gp, NXQ)
%
% This function returns an array of structured element grid points for
% grid point gp of type type using the NXQ array.
%
% type - the type of grid point
% gp - the grid point number
% NXQ - the NXQ array
%
% Mark Trew - June 2001
% 

% Indexing arrays
  MXI1 = 1; PXI1 = 2; MXI2 = 3; PXI2 = 4; MXI3 = 5; PXI3 = 6;
  N = abs(NXQ(2,1,gp)); % this will only be a "correct" value if type=1
  NELEM = [N,6,4,5,6];
  S(1,1:N+1,1:2)=[2*ones(N+1,1),[[1:1:N]';1]];
  S(2,1:6+1,1:2)=[PXI1,1;PXI1,2;PXI1,3;PXI2,1;MXI1,1;MXI2,1;PXI1,1];
  S(3,1:4+1,1:2)=[PXI1,1;PXI2,1;MXI1,1;MXI2,1;PXI1,1];
  S(4,1:5+1,1:2)=[PXI1,1;PXI2,1;MXI1,2;MXI1,1;MXI2,1;PXI1,1];
  S(5,1:6+1,1:2)=[PXI1,1;PXI1,2;PXI2,1;MXI1,2;MXI1,1;MXI2,1;PXI1,1];

  % determine if a layer of elements below the grid pt can be constructed
  off = 0;
  if NXQ(5,1,gp) ~= 0
    for i=1:NELEM(type),
      elements(i,1) = gp;
      elements(i,2) = NXQ(S(type,i,1),1+S(type,i,2),gp);
      elements(i,3) = gp;
      elements(i,4) = NXQ(S(type,i+1,1),1+S(type,i+1,2),gp);

      elements(i,5) = NXQ(5,2,gp);
      elements(i,6) = NXQ(S(type,i,1),1+S(type,i,2),NXQ(5,2,gp));
      elements(i,7) = NXQ(5,2,gp);
      elements(i,8) = NXQ(S(type,i+1,1),1+S(type,i+1,2),NXQ(5,2,gp));
    end;
    off = NELEM(type);
  end;
  
  % determine if a layer of elements above the grid pt can be constructed
  if NXQ(6,1,gp) ~= 0
    for i=1:NELEM(type),
      elements(off+i,1) = gp;
      elements(off+i,2) = NXQ(S(type,i,1),1+S(type,i,2),gp);
      elements(off+i,3) = gp;
      elements(off+i,4) = NXQ(S(type,i+1,1),1+S(type,i+1,2),gp);

      elements(off+i,5) = NXQ(6,2,gp);
      elements(off+i,6) = NXQ(S(type,i,1),1+S(type,i,2),NXQ(6,2,gp));
      elements(off+i,7) = NXQ(6,2,gp);
      elements(off+i,8) = NXQ(S(type,i+1,1),1+S(type,i+1,2),NXQ(6,2,gp));
    end;
  end;

return;




