function [Nk,colNk,GradNk,colGradNk,w] = FEMGauss3on27
% [Nk,GradNk,w] = gauss3on27
% Computes the interpolation functions, their gradients,
% and weight values at 3x3x3 Gauss points in a 27 node master element.
% The field Nk (27*27,1) contains the values of the 27 interpolation
% functions on the 27 Gauss points. 
% Nk(1:27,1) contains the values of the 27 Shape functions at the first
% Gauss point, Nk(28:54,1) contains the 27 Shape functions at the second
% Gauss point, etc. (the blocking: Shape functions within Gausspoints)
% GradNk(27*27,3): Gradients of the shape functions. E.g.,
% GradNk(1:27,1:3): the gradients of the 27 shape functions at the 1st GaussPt
% GradNk(28:54,1:3): the gradients of the 27 shape functions at the 2nd GaussPt

ix = Q27index;
a = sqrt(3/5);     % three point rule

xp = a*ix;        % array of Gauss point coordinates in master element.
ax = [5/9,8/9,5/9]'; 
go = ix+2;         % turned into an index matrix acting on ax.
w = ax(go(:,1)).*ax(go(:,2)).*ax(go(:,3));   % weights attached to gausspt. 
cx =[ 0, -0.5, 0.5;
      1,   0 ,  -1;
      0,  0.5, 0.5];     % coefficient matrix.
cxd = [ -0.5, 1;
	0  ,-2;
	0.5, 1];

pol1 = cx*[ones(1,27);xp(:,1)';(xp(:,1)').^2];  % [1 x x*x];
pol2 = cx*[ones(1,27);xp(:,2)';(xp(:,2)').^2];  % [1 y y*y];
pol3 = cx*[ones(1,27);xp(:,3)';(xp(:,3)').^2];  % [1 z z*z];
pol1d = cxd*[ones(1,27);xp(:,1)'];  % [1 x x*x];
pol2d = cxd*[ones(1,27);xp(:,2)'];  % [1 y y*y];
pol3d = cxd*[ones(1,27);xp(:,3)'];  % [1 z z*z];

Nk=pol1(go(:,1),:).*pol2(go(:,2),:).*pol3(go(:,3),:);

x = pol1d(go(:,1),:).*pol2(go(:,2),:).*pol3(go(:,3),:);
y = pol1(go(:,1),:).*pol2d(go(:,2),:).*pol3(go(:,3),:);
z = pol1(go(:,1),:).*pol2(go(:,2),:).*pol3d(go(:,3),:);

colGradNk = [x(:),y(:),z(:)];
bloc = [x;y;z];
GradNk = zeros(27,81);
GradNk(:) = bloc(:);

colNk = Nk(:);

return;