function StElem(GridPoint,Ngp,Nxi3,dataset)

% Display the closure element

% Nodal lists
if dataset=='sphere'
  SphereDat
elseif dataset=='cylind'
  CylinderDat
elseif (dataset == 'invert')
  InvertDat
elseif (dataset == 'apexda')
  ApexProcData
end

% Nodal locations
nloc = [x(:,1)';y(:,1)';z(:,1)'];    

% Number of faces in xi2 direction
[nc,nn]=size(nloc);
N = (nn-2)/2;

%%%%%%% Construct NXQ and use to determine structured elements %%%%%%%%
% Real element vertices
[GP,GPxi1,GPxi2,GPxi3]=TestNXQ(GridPoint,N,Ngp,Nxi3);
[nv,nl]=size(GPxi1);
for l=1:nl,
  for v=1:nv,
    GPx(v,l) = EvalClosureElement(GPxi1(v,l),GPxi2(v,l),GPxi3(v,l),x,N);
    GPy(v,l) = EvalClosureElement(GPxi1(v,l),GPxi2(v,l),GPxi3(v,l),y,N);
    GPz(v,l) = EvalClosureElement(GPxi1(v,l),GPxi2(v,l),GPxi3(v,l),z,N);
  end;
end;

% grid points
[ngp,nc]=size(GP);
for gp=1:ngp,
  gpx(gp) = EvalClosureElement(GP(gp,1),GP(gp,2),GP(gp,3),x,N); 
  gpy(gp) = EvalClosureElement(GP(gp,1),GP(gp,2),GP(gp,3),y,N); 
  gpz(gp) = EvalClosureElement(GP(gp,1),GP(gp,2),GP(gp,3),z,N); 
end;


% Plot the elements and grid points
clf;
hold on;
plot3(GPx,GPy,GPz,'b-','linewidth',1);
plot3(gpx,gpy,gpz,'ko','linewidth',2);
view(3);
hold off;
axis equal;
grid on;
xlabel('x');ylabel('y');zlabel('z');

%%%%%%% Calculate and draw surfaces of closure element %%%%%%%%

%xi locations 
xi1 = 0:0.1:1;
xi2 = 0:0.05:1;
xi3 = 0:0.1:1;

% grid of xi1 and xi2 points
[mxi1,mxi2]=meshgrid(xi1,xi2);

% bottom surface
[nxi1,nxi2]=size(mxi1);
for i=1:nxi1,
   for j=1:nxi2,
         xsurf1(i,j) = EvalClosureElement(mxi1(i,j),mxi2(i,j),0.0,x,N);
         ysurf1(i,j) = EvalClosureElement(mxi1(i,j),mxi2(i,j),0.0,y,N);
         zsurf1(i,j) = EvalClosureElement(mxi1(i,j),mxi2(i,j),0.0,z,N);
   end;
end;

% top surface
for i=1:nxi1,
   for j=1:nxi2,
         xsurf2(i,j) = EvalClosureElement(mxi1(i,j),mxi2(i,j),1,x,N);
         ysurf2(i,j) = EvalClosureElement(mxi1(i,j),mxi2(i,j),1,y,N);
         zsurf2(i,j) = EvalClosureElement(mxi1(i,j),mxi2(i,j),1,z,N);
   end;
end;

% outside boundary - xi1 = 1
for i=1:length(xi3),
   for j=1:length(xi2),
         xsurf3(j,i) = EvalClosureElement(1,xi2(j),xi3(i),x,N);
         ysurf3(j,i) = EvalClosureElement(1,xi2(j),xi3(i),y,N);
         zsurf3(j,i) = EvalClosureElement(1,xi2(j),xi3(i),z,N);
         csurf3(j,i) = xi3(i);
   end;
end;

xsurf = [xsurf1,xsurf2,xsurf3];
ysurf = [ysurf1,ysurf2,ysurf3];
zsurf = [zsurf1,zsurf2,zsurf3];

csurf = [mxi1,mxi2,csurf3];

figure(2);
cval = 0.*zsurf;
surf(xsurf,ysurf,zsurf,csurf);
shading interp;
axis equal;
hold on;
mesh(xsurf,ysurf,zsurf,cval);
hidden off;
plot3(x(:,1),y(:,1),z(:,1),'kp');
hold off;
xlabel('x');
ylabel('y');
zlabel('z');

figure(3);
plot3(xsurf,ysurf,zsurf,'ko','linewidth',2);
view(3);
