function NXQ = CreateNXQ(N,Ngp,Nxi3)

%
% NXQ = CreateNXQ(N,Ngp,Nxi3)
%
% This function creates the NXQ array within a closure element - 
% note that in CMISS the NXQ indices in each direction 
% are -1,1,-2,2,-3,3 - the corresponding 
% indices for the MATLAB implementation are MXI1,PXI1,MXI2,PXI2,MXI3,PXI3.
%
% N - number of perimeter nodes on closure element
% Ngp - number of grid points between each pair of perimeter 
%       nodes
% Nxi3 - number of grid points in the Xi3 direction
%
% Mark Trew - June 2001
%

% set up constants for xi direction indices
  MXI1 = 1; PXI1 = 2; MXI2 = 3; PXI2 = 4; MXI3 = 5; PXI3 = 6;

% set up lookup tables
% Type: 1  2  3  4  5
  D = [ 0 -1 -1 -2 -2  % -xi1 direction
       -N -3  0  0 -2  % +xi1 direction - zeros placed for pts outside closure element
        0  1  1  1  1  % -xi2 direction
        0  1  1  1  1  % +xi2 direction
      ];
  % Mapping array for gp type 1
  M(1,1:N,1:5) = ...
       [
        PXI1*ones(N,1),[1:1:N]',ones(N,1),zeros(N,1),[0:1:N-1]'
       ]; 
  % Mapping array for gp type 2
  M(2,1:6,1:5) = ...
       [
        MXI1 1 -1 -1  0
        PXI1 1  1  1 -1
        PXI1 2  1  1  0
        PXI1 3  1  1  1
        MXI2 1  0  0 -1
        PXI2 1  0  0  1
       ];
  % Mapping array for gp type 3 
  M(3,1:4,1:5) = ...
       [
        MXI1 1 -1 -1  0
        PXI1 1  0  0  0  % zeros placed for grid points outside closure element
        MXI2 1  0  0 -1
        PXI2 1  0  0  1 
       ]; 
  % Mapping array for gp type 4
  M(4,1:5,1:5) = ...
       [
        MXI1 1 -1 -1 -1
        MXI1 2 -1 -1  0
        PXI1 1  0  0  0  % zeros placed for grid points outside closure element
        MXI2 1  0  0 -1
        PXI2 1  0  0  1
       ]; 
  % Mapping array for gp type 5
  M(5,1:6,1:5) = ...
       [
        MXI1 1 -1 -1 -1
        MXI1 2 -1 -1  0
        PXI1 1  1  1  0
        PXI1 2  1  1  1
        MXI2 1  0  0 -1
        PXI2 1  0  0  1
       ];
  % Number of adjacent grind points for each gp type 
  NADJ=[N,6,4,5,6];

% loop over layers in xi3
  for k=1:Nxi3, % layers
    
  % Calculate layer offset
    Lgp = 1 + (Ngp*(Ngp-1)*N)/2;

    % fill the closure plane NXQ
    for c=1:Ngp,
      for i=1:1-rem(1,c)+(c-1)*N,
        
        % determine type of grid point
        type = (1-rem(1,c))*1+...                                                   % type 1
               (rem(1,Ngp+1-c))*rem(1,c)*(1-min(rem(i-1,1-rem(1,c)+c-1),1))*2+... % type 2 
               (1-rem(1,Ngp+1-c))*(1-min(rem(i-1,1-rem(1,c)+c-1),1))*3+...        % type 3
               (1-rem(1,Ngp+1-c))*(min(rem(i-1,1-rem(1,c)+c-1),1))*4+...          % type 4 
               (rem(1,Ngp+1-c))*(min(rem(i-1,1-rem(1,c)+c-1),1))*5;               % type 5
        
        % determine current grid point number
        gpci = GP(c,i,N);
        gpci = gpci + (k-1)*Lgp; % account for layer
        
        % determine sector number for c and i
        s = fix(((i-1)*N)/(1-rem(1,c)+(c-1)*N)) + N*(max(1-i,0)+min((c-1)*N-i,0));

        % fill "zero" NXQ index
        NXQ(MXI1,1,gpci) = D(1,type);
        NXQ(PXI1,1,gpci) = D(2,type);
        NXQ(MXI2,1,gpci) = D(3,type);
        NXQ(PXI2,1,gpci) = D(4,type);
        NXQ(MXI3,1,gpci) = rem(1,k); % zero if k==1
        NXQ(PXI3,1,gpci) = rem(1,Nxi3+1-k); % zero if k==Nxi3

        % fill adjacent grid points
        for g = 1:NADJ(type),
          C = c + M(type,g,3);
          I = i + M(type,g,4)*s + M(type,g,5);
          gpCI = GP(C,I,N);
          gpCI = gpCI + (k-1)*Lgp; % account for layer
          NXQ(M(type,g,1),1+M(type,g,2),gpci) = gpCI;
        end;
        
        % fill top and bottom grid points
        NXQ(MXI3,2,gpci) = NXQ(MXI3,1,gpci)*(gpci - Lgp); % zero if no gp below       
        NXQ(PXI3,2,gpci) = NXQ(PXI3,1,gpci)*(gpci + Lgp); % zero if no gp above        

      end;
    end;
  end;

return;
