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

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

%xi locations 
xi1 = 0:0.05:1;
xi2 = 0:0.0625:1;

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

% 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);
         u(i,j) = EvalClosureElement(mxi1(i,j),mxi2(i,j),1,z,N);
   end;
end;

mesh(xsurf2,ysurf2,u);
axis equal;
xlabel('x');
ylabel('y');
zlabel('u');
view(2);







