function gplocs = CreateGPLocs(N,Ngp,Nxi3)

%
% gplocs = CreateGPLocs(N,Ngp,Nxi3)
% 
% This function generates the grid point (xi1, xi2, xi3) locations in a 
% closure element.
%
% N - number of perimeter nodes
% Ngp - number of grid points on each face
% Nxi3 - number of grid points in the xi3 direction
%
% Mark Trew - June 2001
%

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

    % loop over rows and columns
    for c=1:Ngp,
      for i=1:1-rem(1,c)+(c-1)*N,
        gp = GP(c,i,N);
        gp = gp + (k-1)*Lgp;

        gplocs(gp,:) = [(c-1)/(Ngp-1),(i-1)/(1-rem(1,c)+(c-1)*N),(k-1)/(1-rem(1,k)+Nxi3-1)];
      end;
    end;
  end;

return;